! Copyright 2008 Matthew Norton ! ! This program is free software; you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation; either version 2 of the License, or ! (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with this program; if not, write to the Free Software ! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA ! ! NOTE: This is in FORTRAN 90 and you will need a FORTRAN 90 complier to run this code ! ------------------------------------------------------------------------------ ! This program apporoximates the Lyapunov Characteristic Exponent of the neutral ! line magnetic field ! ------------------------------------------------------------------------------ PROGRAM DrivenPendulum ! IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT NONE INTEGER, PARAMETER :: neq = 4, maxord = 12 INTEGER, PARAMETER :: liw = 20 + neq, NumSteps = 10000000 INTEGER, PARAMETER :: lrw = 22 + neq*(maxord+1) + 3*neq + neq**2 REAL*8, PARAMETER :: D0 = 0.01 INTEGER :: meth, miter, mf REAL*4 :: E = 0.5d0, Bn = 0.275d0, Pz = -0.9215d0 REAL*8 :: rwork(lrw), Energy, Veff, Dk, LambdaK = 0.0d0 INTEGER :: iwork(liw), ErrorStat, I REAL*8 :: Y(neq), X0, Y0, Vx0, Vy0, TwoPi, H, Td, Time, T2, Wg, H0 REAL*8 :: Lx, Ly, LVx, xs, ys, Vxs, atol, rtol, Pi, YT(neq) INTEGER :: NFrac, M, itol, itask, istate, iopt, J, K, L, Tau, TCount = 1 REAL(8), DIMENSION(NumSteps) :: Xin, Yin, Vxin, Vyin EXTERNAL f,jac COMMON/DERIV/ E, Bn, Pz OPEN(UNIT = 20, FILE = "Neutral.txt", STATUS = "UNKNOWN", IOSTAT = ErrorStat) IF (ErrorStat /= 0) STOP "*** Error in opening output file 1! ***" OPEN(UNIT = 30, FILE = "NeutralIn.txt", STATUS = "UNKNOWN", IOSTAT = ErrorStat) IF (ErrorStat /= 0) STOP "*** Error in opening input file! ***" OPEN(UNIT = 40, FILE = "LCE.txt", STATUS = "UNKNOWN", IOSTAT = ErrorStat) IF (ErrorStat /= 0) STOP "*** Error in opening output file 2! ***" ! ! ...... INPUT: parameters and initial values ! CALL GetOld() Tau = NumSteps / 1000000 Pi = 4.0D0*DATAN(1.0D0) NFrac = 10000 X0 = 2.00d0 + D0 Vx0 = 0.5d0 ! PRINT *, X0, Vx0 Y0 = DSQRT(-2.0d0 * Pz) Vy0 = DSQRT(8.0d0 * E - 4.0d0 * Pz**2 - 4.0d0 * Vx0 ** 2 + 4.0d0 * Bn ** 2 * Pz* & & X0 ** 2 - Bn ** 4 * X0 ** 4 - 4.0d0 * Pz * Y0 ** 2 + 2.0 * Bn ** 2 * X0& & ** 2 * Y0 ** 2 - Y0 ** 4) / 2.0d0 ! ! ! ...... COMPUTE timestep h to be a fraction of the natural period, and ! set up initial arrays ! TwoPi = 8.0D0*DATAN(1.0D0) H0 = TwoPi / NFrac !H = H0 Time = 0.0d0 y(1) = x0 ! X y(2) = y0 ! Y y(3) = Vx0 ! X dot y(4) = Vy0 ! Y dot ! ...... DLSODE parameters (see users manual) ! itol = 1 rtol = 1.d-6 atol = 1.d-10 itask = 1 istate = 1 iopt = 0 meth = 1 miter = 1 mf = 10*meth + miter DO I = 1, NumSteps Wg = DSQRT(y(2) ** 2 + Bn ** 4 * y(1) ** 2) ! gyro frequency H = TwoPi / (NFrac * Wg) ! calculate step size IF(H > H0) THEN H = H0 ! If step size is greater than max value, use max value END IF !Print *, H T2 = Time + H ! update time CALL dlsode(f,neq,y,Time,T2,itol,rtol,atol,itask,istate, & & iopt,rwork,lrw,iwork,liw,jac,mf) ! ! ...... CHECK dlsode error conditions ! 140 FORMAT(d15.5,d16.5,d14.3,i5,d14.3) IF (istate .lt. 0) THEN PRINT*, 'istate = ',istate,' stopping run' PRINT *, "time", I ! PAUSE EXIT END IF ! ! ...... OUTPUT ! ! test for energy conservation Veff = 1.0d0/2.0d0 * (Pz + 1.0d0/2.0d0 * (y(2) ** 2 - Bn ** 2 * y(1) ** 2)) ** 2 Energy = 1.0d0/2.0d0 * (y(3) ** 2 +y(4) ** 2) + Veff !PRINT *, Energy IF (ABS(Energy - E) >= 0.005) THEN PRINT *, "ERROR Erergy not conserved!", Energy, I EXIT END IF Time = T2 IF (MOD(I,Tau) == 0) THEN Dk = DSQRT((Y(1) - Xin(I)) ** 2 + (Y(2) - Yin(I)) ** 2 + (Y(3) - & &Vxin(I)) ** 2 + (Y(4) - Vyin(I)) ** 2) PRINT *, Dk Y(1) = Xin(I) + D0 / Dk * Y(1) Y(2) = Yin(I) + D0 / Dk * Y(2) Y(3) = Vxin(I) + D0 / Dk * Y(3) Y(4) = Vyin(I) + D0 / Dk * Y(4) LambdaK = LambdaK + DLOG(Dk / D0) TCount = TCount + 1 WRITE(40,*) I, 1.0d0 / DBLE(Tcount * Tau) * LambdaK END IF WRITE(20,*) Y(1), Y(2), Y(3), Y(4) IF (MOD(I,1000000) == 0) PRINT *, I END DO ! End trajectory loop PRINT *, "LOOP EXITED" LambdaK = LambdaK / DBLE(TCount * Tau) PRINT *, "Lyapunov Characteristic Exponent is: ", LambdaK CLOSE(20) !200 FORMAT(4(1X, F12.9)) STOP ! ------------------------------------------------------------------------------ CONTAINS ! ------------------------------------------------------------------------------ SUBROUTINE GetOld() REAL*8 :: Dummy INTEGER :: ErrorStat, I DO I = 1, NumSteps READ(30,*) Xin(I), Yin(I), Vxin(I), Vyin(I) END DO CLOSE(30) RETURN !200 FORMAT(4(1X, F12.9)) END SUBROUTINE END PROGRAM ! ! ...... THIS subroutine computes the equations of motion for the ! neutral line field ! SUBROUTINE f(neq,t,y,ydot) real*8 :: y(neq),ydot(neq),t COMMON/DERIV/ E, Bn, Pz ydot(1) = y(3) ! xdot ydot(2) = y(4) ! ydot ydot(3) = Bn ** 2 * y(1) * ((1.0d0/2.0d0 * (y(2) ** 2 - Bn ** 2 * y(1) ** 2)) + Pz) !vxdot ydot(4) = -1.0d0 * y(2) * ((1.0d0/2.0d0 * (y(2) ** 2 - Bn ** 2 * y(1) ** 2)) + Pz) ! vydot RETURN END ! ! ! ...... THIS subroutine computes the jacobian of the equations of ! motion: partial derivative pd(i,j) = d(ydot(i))/dy(j) ! SUBROUTINE jac(neq,t,y,ml,mu,pd,nrowpd) real*8 :: y(neq),pd(nrowpd,*),t COMMON/DERIV/ E, Bn, Pz pd(1,1) = 0.0 pd(1,2) = 0.0 pd(1,3) = 1.0 pd(1,4) = 0.0 pd(2,1) = 0.0 pd(2,2) = 0.0 pd(2,3) = 0.0 pd(2,4) = 1.0 PD(3,1) = -(Bn ** 4 * Y(1) ** 2) + Bn ** 2 * (Pz + (-(Bn ** 2 * Y(1) ** 2) + Y(2) ** 2)/2.0d0) pd(3,2) = Bn ** 2 * y(2) * y(1) pd(3,3) = 0.0 pd(3,4) = 0.0 pd(4,1) = Bn ** 2 * y(1) * y(2) pd(4,2) = -Pz - Y(2) ** 2 + (Bn ** 2 * Y(1) ** 2 - Y(2) ** 2) / 2.0d0 pd(4,3) = 0.0 pd(4,4) = 0.0 RETURN END