! 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 motion of a charged particle in a magnetic ! neutral line field using DLSODE solver ! ------------------------------------------------------------------------------ PROGRAM DrivenPendulum ! IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT NONE INTEGER, PARAMETER :: neq = 4, maxord = 12 INTEGER, PARAMETER :: liw = 20 + neq INTEGER, PARAMETER :: lrw = 22 + neq*(maxord+1) + 3*neq + neq**2 INTEGER :: meth, miter, mf REAL*4 :: E = 0.5d0, Bn = 0.275d0, Pz = -0.9215d0 REAL*8 :: rwork(lrw), Ysos, Energy, Veff INTEGER :: iwork(liw), ErrorStat, I, NumSteps, Flag = -1 REAL*8 :: Y(neq), Yold(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 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! ***" OPEN(UNIT = 30, FILE = "Neutral2.txt", STATUS = "UNKNOWN", IOSTAT = ErrorStat) IF (ErrorStat /= 0) STOP "*** Error in opening output file! ***" ! ! ...... INPUT: parameters and initial values ! NumSteps = 100000 Pi = 4.0D0*DATAN(1.0D0) DO J = -10, 10 ! DO K = -24, 24 Flag = 1 NFrac = 10000 X0 = DBLE(J) / 5.0d0 ! Vx0 = DBLE(K) / 25.0d0 ! X0 = 3.5d0 Vx0 = 0.0d0 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 Ysos = DSQRT(-2.0d0 * Pz) yold(2) = Ysos yold(1) = X0 yold(3) = Vx0 yold(4) = Vy0 ! ! ! ...... 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 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.05) THEN PRINT *, "ERROR Erergy not conserved!", Energy EXIT END IF Time = T2 ! WRITE(30,'(5(E12.5,3X))') Time, y(1), y(2), y(3), y(4) IF (Yold(2) < YSos .AND. Y(2) > YSos) THEN DO L = 1, neq yt(L) = yold(L) END DO CALL Henon(yt, YSos) END IF DO M = 1, Neq Yold(M) = Y(M) END DO IF (K == 0 .AND. I == 1) THEN PRINT *, J , K ! print where program is END IF ! IF (MOD(I,1000) == 0) PRINT *, I END DO ! End trajectory loop END DO ! END DO CLOSE(20) STOP ! ------------------------------------------------------------------------------ CONTAINS ! ------------------------------------------------------------------------------ SUBROUTINE Henon(Y, sos) INTEGER, PARAMETER :: neqs = 4 REAL*8 :: y(neqs), t, h, Veff, Energy, sos EXTERNAL f2,jac2 ! itol = 1 ! rtol = 1.d-6 ! atol = 1.d-10 ! itask = 1 ! istate = 1 ! iopt = 0 ! meth = 1 ! miter = 1 ! mf = 10*meth + miter t = y(2) t2 = sos CALL dlsode(f2,neqs,y,T,T2,itol,rtol,atol,itask,istate, & & iopt,rwork,lrw,iwork,liw,jac2,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 *, "Henon" PAUSE END IF ! 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! Henon", Energy ! PAUSE ! END IF WRITE(30,*) Y(1), Y(2), Y(3), Y(4) RETURN END SUBROUTINE END ! ! ! ...... THIS subroutine computes the equations of motion for the ! simple harmonic oscillator (mass on spring) ! 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 SUBROUTINE f2(neqs,t,y,ydot) real*8 :: y(neqs),ydot(neqs),t COMMON/DERIV/ E, Bn, Pz ydot(1) = y(3) / y(4) ! xdot ydot(2) = 1.0d0 ! ydot ydot(3) = (Bn ** 2 * y(1) * ((1.0d0/2.0d0 * (y(2) ** 2 - Bn ** 2 * y(1) ** 2)) + Pz)) / y(4) !vxdot ydot(4) = (-1.0d0 * y(2) * ((1.0d0/2.0d0 * (y(2) ** 2 - Bn ** 2 * y(1) ** 2)) + Pz)) / y(4) ! vydot RETURN END SUBROUTINE SUBROUTINE jac2(neq,t,y,ml,mu,pd,nrowpd) real*8 :: y(neq),pd(nrowpd,*),t COMMON/DERIV/ E, Bn, Pz pd(1,1) = 0.0d0 pd(1,2) = 0.0d0 pd(1,3) = 1.0d0 / Y(4) pd(1,4) = -1.0d0 / Y(4) ** 2 pd(2,1) = 0.0d0 pd(2,2) = 0.0d0 pd(2,3) = 0.0d0 pd(2,4) = 0.0d0 PD(3,1) = -((Bn ** 4 * Y(1) ** 2) / Y(4)) + (Bn ** 2 * (Pz + (-(Bn ** 2 *Y(1) ** 2) + Y(2) ** 2) / 2.0d0)) / Y(4) pd(3,2) = (Bn ** 2 * Y(1) * Y(2)) / Y(4) pd(3,3) = 0.0d0 pd(3,4) = -((Bn ** 2 * Y(1) * (Pz + (-(Bn ** 2 * Y(1) ** 2) + Y(2) ** 2) / 2.0d0)) / Y(4) ** 2) pd(4,1) = (Bn ** 2 * Y(1) * Y(2)) / Y(4) pd(4,2) = -(Y(2) ** 2 / Y(4)) - (Pz + (-(Bn ** 2 * Y(1) ** 2) + Y(2) ** 2) / 2.0d0) / Y(4) pd(4,3) = 0.0d0 pd(4,4) = (Y(2) * (Pz + (-(Bn ** 2 * Y(1) ** 2) + Y(2) ** 2) / 2.0d0)) / Y(4) ** 2 RETURN END