| Version 7 (modified by , 8 years ago) (diff) |
|---|
-
FORTRAN 95 examples
- 1. Integral function calculation using trapeze, rectangle and Simpson …
- 2. Integral function calculation using Simspon method
- 3. Calculation dependence between integral function and coordinates
- 4. Calculation dependence between integral function and coordinates
- 5. Simpson formula for Simpson method of integration
- 6. Characteristics of the electron motion, which is located between …
- 7. Electron motion only in magnetic field
- 8. Electron motion in electric and magnetic filed.
- 9. Electron motion in electric and magnetic filed
- 10. Electron motion only in magnetic field
- 11. Determination the proton energy and power
- 12. Runge-Kutta method for harmonic oscillations
- 13. Runge-Kutta method for harmonic oscillations
- 14. Runge-Kutta method for harmonic oscillations
- 15. Electron acceleration in uniform electric field
- 16. Determination the conditions of electron capture in SGA mode and …
- 17. Determination the dependencies between the electron momentum and …
- 18. Electron acceleration in a uniform electric field using speed light
- 19. Charged particle motion in magnetic mirror trap under ECR
- 20. Solving nonlinear equation
- 21. Poisson equations solver
- 22. Distribution of Langmuir waves in a uniform plasma
- 23. Proton acceleration by the filed of electron leyer
- 24. Electron motion in uniform electric filed
FORTRAN 95 examples
1. Integral function calculation using trapeze, rectangle and Simpson method
real x,dx,s a=-0.5 b=0.5 n=100 dx=(b-a)/n S1=0 S2=0 S3=0 do i=0,n S1=S1+(f(i*dx+a)+f((i+1)*dx+a))/2*dx !integral which is calculate by trapeze method S2=S2+ f(i*dx+a)*dx !integral which is calculate with rectangle method if(mod(i,2)==0) then S3=S3+dx/3*2*f(i*dx+a) !integral which is calculate by Simpson method else S3=S3+dx/3*4*f(i*dx+a) end if end do print *, 'Trap = ', S1 , 'Rec=', S2 , 'Sim=', S3 pause contains function f(x) f=1/sqrt(1-x**2) !integral function end function f end
2. Integral function calculation using Simspon method
real x,dx,s a=0 b=1.57 n=100 dx=(b-a)/n S=0 do i=0,n if(mod(i,2)==0) then S=S+dx/3*2*f(i*dx+a) !integral which is calculate by Simpson method else S=S+dx/3*4*f(i*dx+a) end if end do print *, 'Sim = ', S pause contains function f(x) f=sin(x)*sin(2*x)*sin(3*x) !integral function end function f end
3. Calculation dependence between integral function and coordinates
! Get a graph of the dependence in Origin, differentiate and integrate the graph and !compare the resultants with the integrals which are calculated in the first program real x,dx, a, b open(1,file='y2.dat', status='unknown')!the document y2 is in the project folder a=-0.5 b=0.5 n=100 dx=(b-a)/n do i=1,n write(1,*) i*dx + a, 1/sqrt(1-(i*dx+a)**2)!integral function end do end
4. Calculation dependence between integral function and coordinates
! Get a graph of the dependence in Origin, differentiate and integrate the graph and !compare the resultant with the integral which is calculated in the second program real x,dx, a, b open(1,file='y2.dat', status='unknown')!the document y2 is in the project folder a=0 b=1.14 n=100 dx=(b-a)/n do i=1,n write(1,*) i*dx + a, sin(i*dx+a)*sin(2*i*dx+a)*sin(3*i*dx+a)!integral function end domethod of end
5. Simpson formula for Simpson method of integration
program Simpson integer n real y1, y2, y3, x, h, a, s n=100 h=2/n do i=1,n x(i)=-1+h(i-1) y1(i)=f1(x(i)) y2(i)=f2(x(i),h) y3(i)=f3(x(i),h) a=h/6*(y1+4*y2+y3)!Simpson's formula s=sum(a) end do contains real function f1(x) real x f1=x*(5-4*x)**(-0.5) !first function end function real function f2(x,h) real x, h f2=(x+h/2)*(5-4*(x+h/2))**(-0.5) !second function end function real function f3(x,h) real x, h f3=(x+h)*(5-4*(x+h))**(-0.5) !third function end function print *, smethodof !Simpson method end program
6. Characteristics of the electron motion, which is located between charged coaxial rings
!Between two coaxial rings, charged with q1=q2=1nC and their radius are equal to r1=1cm !and r2=2cm, are located at distance 5 cm between each other. At halfway between the rings is !located electron with initial energy equal to 0. Find the dependencies between the !velocity(energy)and the time,and electron coordinate with time. Also plot the phase !trajectory V(x). And explain the resultants. program rings integer :: n=10000 real(8) :: r1=1.0, r2=2.0, q=3.0, d=5.0, eq=4.8e-10, m=9.1e-28 !the unis are in CGS system real(8) :: dt, t, x, vx t=1.0e-6 dt=t/float(n) vx=0 x=2.5 open (1,file='rings',status='unknown') do i=1,n E=q*x/sqrt((x**2+r1**2)**3)- q*(5-x)/sqrt(((5-x)**2+r2**2)**3) !energy vx=vx+(eq*E/m)*dt !phase trajectory x=x+vx*dt !coordinate write(1,*) x, i*dt, vx end do end program
7. Electron motion only in magnetic field
!Using the equation of electron motion find the phase trajectory. Electron filed E=0 !Part 1. Using only dependence of B0 integer :: n=1000 real(8) :: dt, t, vm, b, rl, om, x, y, vx, vy, B0 real(8) :: eq=-4.8e-10, me=9.1e-28, c = 3.0e10 !B0=me*c*om/ez B0 = 200.0 !E0 = 5000 !A = 1e-30 x=1.0 y=1.0 vx=1000.0 vy=1000.0 a=100 !b = B0/B0 !om = eq*B0/me*c !vm = A*om !print *, vm !rl = c/om !vx = vm/vm !vy = vm/vm !x = x/rl !y = y/rl !print *, rl t = 1.0e-8 dt = t/float(n) open (1,file='1',status='unknown') do i=1,n vx = vx+(B0*eq)/(c*me)*vy*dt vy = vy-(B0*eq)/(c*me)*vx*dt x = x+vx*dt y = y+vy*dt write(1,*) x, vx end do pause end program
8. Electron motion in electric and magnetic filed.
!Like in program 7.Part 1 integer :: n=1000 real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0 x = 1.0 y = 1.0 a = 0.5 vx = 1.0e09 vy = 1.0e09 E0=0.0 B0 = 200.0 b=(2*3.14*6)/(1.0e09) !B0 = 1e-25 !E0 = 5000 !A = 1e-30 !b = B0/B0 !om = eq*B0/me*c !vm = A*om !print *, vm !rl = c/om !vx = vm/vm !vy = vm/vm !x = x/rl !y = y/rl !print *, rl t = 1.00e-08 dt = t/float(n) open (1,file='1',status='unknown') do i=1,n Bz = B0*(1+a*y) !Bz = B0 x = x + vx*dt y = y + vy*dt !vx = vx + (vy*eq*Bz)/(c*me)*dt !vy = vy - (vx*eq*Bz)/(c*me)*dt if(i*dt.gt.b) then E0 = 5.0 end if vx = vx +(eq*E0)/(me)*dt+(vy*eq*Bz)/(c*me)*dt vy = vy + (eq*E0)/(me)*dt-(vx*eq*Bz)/(c*me)*dt ! x = x + vx*dt !y = y + vy*dt !print *, vx, vy write (1,*) x, y end do pause end program
9. Electron motion in electric and magnetic filed
!Like in program 8. Part 2. the differnece is that here we are using different functions !for vx and vy integer :: n=1000 real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0 x = 1.0 y = 1.0 a = 0.5 vx = 1.0e09 vy = 1.0e09 E0 = 0.0 B0 = 200.0 b = (2*3.14)/(1.0e09) t = 1.00e-08 dt = t/float(n) open (1,file='1',status='unknown') do i=1,n Bz = B0*(1+a*y) !Bz = B0 x = x + vx*dt y = y + vy*dt vx = vx + (vy*eq*Bz)/(c*me)*dt vy = vy - (vx*eq*Bz)/(c*me)*dt if(i*dt.gt.b) then E0 = 3.0 end if vx = vx + (vy*eq*Bz)/(c*me)*dt vy = vy + (eq*E0)/(me)*dt - (vx*eq*Bz)/(c*me)*dt x = x - vx*dt y = y + vy*dt !print *, vx, vy write (1,*) x, y end do pause end program
10. Electron motion only in magnetic field
! Like in program 7. Part 2 here we have dependence of magnetic filed in z coordinate integer :: N = 1500 real(8) :: x = 2.5, y = 0.5, B0 = 500.0, We = 4.0054e-08, b = 2.0, pi = 3.14159, m = 9.1e-28, e = -4.8e-10, c = 2.99e10 real(8) vx, vy, v, dt, t v = sqrt(2*We/m) vx = sin(pi/4.0)*v vy = cos(pi/4.0)*v t = 5e-9 dt=t/float(N) !Norma !tau = t*10e8 !me = m*1e28 !qe = e*1e10 !vl = c*1e-10 !ve = v/c !vx = sin(pi/4.0)*v !vy = cos(pi/4.0)*v !dt = tau/float(N) open (1,file='mov',status='replace') do i=1,N Bz = B0*(1-b**2*(x**2+y**2)) x = x + vx*dt y = y + vy*dt vx = vx + vy/m*Bz*e/c*dt vy = vy - vx/m*Bz*e/c*dt print *, x, y pause write(1,*) y, x end do end program
11. Determination the proton energy and power
integer :: n=1000 real(8) :: dt, t, x, y, vx, vy, v, E, W real(8) :: mp = 1.67e-24, q = 4.8e-10 x = 1.0 y = 1.0 d = 1.0 a = 30 W0 = 1000 v = sqrt(2*W0/mp) t = 1.00e-05 dt = t/float(n) open (1,file='1',status='unknown') do i=1,n vx = v*sin(a) vy = v*cos(a) x = x + vx*dt y = y + vy*dt E = (4*3.14*mp*y)/d W = (E*q)/ x !print *, vx, vy write (1,*) W, t end do pause end program
12. Runge-Kutta method for harmonic oscillations
! metod Runge-Kuta. Oscillator. Part 1 external fct, out real aux(8,2) common /a/ b, lam, f0, dt, k, n real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/) real :: u(2) = (/0., 0.5/) real :: du(2) = (/0.5, 0.5/) real :: b=0.003 integer :: nd = 2 real :: dt = 2*3.14159265 ! real :: lam = 0.1 real :: f0 = 0.00 open (unit=2, file = 'osc_rk.dat', status='new') write(*,*) 'x=', u(1), 'v=', u(2) pause call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=', ih stop end subroutine fct (t, u, du) real u(2), du(2), lam common /a/ b, lam, f0, dt, k, n du(1) = -u(2)-2.*b*u(1) du(2) = u(1) return end subroutine out (t, u, du, ih, nd, pr) real u(2), du(2), pr(5), lam common /a/ b, lam, f0, dt, k, n if (t.ge.k*(dt/20.)) then write (2,1) t,u(1),u(2) ! k = k + 1 else end if 1 format (f8.3, f8.3, f8.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
13. Runge-Kutta method for harmonic oscillations
! metod Runge-Kuta. Oscillator. Part 2 (the differnece between part 1 is that here we are ! using onother function for du(1)) external fct, out real aux(8,2) common /a/ b, lam, f0, dt, k, n, c real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/) real :: u(2) = (/0., 0.5/) real :: du(2) = (/0.5, 0.5/) real :: b=0.003 integer :: nd = 2 real:: c = 2.0 real :: dt = 2*3.14159265 ! real :: lam = 0.1 real :: f0 = 0.00 n = 1 k = 1 open (unit=2, file = 'osc_rk.dat', status='new') write(*,*) 'x=', u(1), 'v=', u(2) pause call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=', ih stop end subroutine fct (t, u, du) real u(2), du(2), lam common /a/ b, lam, f0, dt, k, n, c du(1) = c*cos(0.5*t)-u(2)-2.*b*u(1) du(2) = u(1) return end subroutine out (t, u, du, ih, nd, pr) real u(2), du(2), pr(5), lam common /a/ b, lam, f0, dt, k, n, c if (t.ge.k*(dt/20.)) then write (2,1) t,u(1),u(2) ! k = k + 1 else end if 1 format (f8.3, f8.3, f8.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
14. Runge-Kutta method for harmonic oscillations
!Part 3. Another example and using another function for du(1) external fct, out real aux(8,2) common /a/ dt, k, n, A, Vm, om real :: pr(5) = (/0.,2*3.14159*10., 2*3.14159/50, .0001, .0/) real :: du(2) = (/0.5, 0.5/) integer :: nd = 2 real :: dt=2*3.14159265 real :: t A = 2.0 om = 1.0 Vm = A*om u(1) = Vm/Vm u(2) = 0 k = 1 n=1 open (unit=2, file = 'osc_rk0.dat' ,status='unknown') write(*,*) 'x=' , u(1), 'v=', u(2) pause call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=', ih stop end subroutine fct (t, u, du) real u(2), du(2) common /a/ dt, k, n, A, Vm, om du(1) = -u(2)/(1+0.1*sin(om*dt)) du(2) = u(1) return end subroutine out (t, u, du, ih, nd, pr) real u(2), du(2), pr(5) common /a/ dt, k, n, A, Vm, om if(t.ge.k*(dt/50.)) then write (2,1) t/om, u(1)*Vm, u(2)*A k=k+1 else end if 1 format (f8.3, f8.3, f8.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
15. Electron acceleration in uniform electric field
! Using metod Runge-Kuta. . external fct, out real aux(8,4) common /a/ dt, k, n, vm, A, om, rl, x, y, r, u(1), u(2), u(3), u(4) real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/) integer :: nd = 2 real :: dt = 2*3.14159265 real :: t real :: c = 3.0e10 !B0=me*c*om/ez B0 = 200.0 E0 = 5000 eq = 4.8e-10 me = 9.1e-28 om = eq*B0/me*c A = 2.0 vm = A*om n = 1 k = 1 b = B0/B0 rl = c/om u(1) = vm/vm u(2) = vm/vm u(3) = x/rl u(4) = y/rl open (unit=2, file = '1_rk.dat', status='unknown') write(*,*) 'x=', u(3), 'vx=', u(1), 'y=', u(4), 'vy=', u(2) call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=', ih end subroutine fct (t, u, du) real u(3), du(3), u(4), du(4) common /a/ dt, k, n, vm, A, om, rl, x, y, r, u(1), u(2), u(3), u(4) du(1) = -u(3) du(2) = u(4) du(3) = u(1) du(4) = u(2) return end subroutine out (t, u, du, ih, nd, pr) common /a/ dt, k, n, om, vm, A, rl, x, y, r, u(1), u(2), u(3), u(4) if (t.ge.k*(dt/20.)) then write (4,1) t, u(1), u(2), u(3), u(4) k = k + 1 else end if 1 format (f8.3, f8.3, f8.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
16. Determination the conditions of electron capture in SGA mode and the bunch time
! Using metod Runge-Kuta. Resonant and auto resonant. ! SGA (synchrotron gyro-magnetic auto resonant) is a self sustaining ECR plasma in !magnetic field which is rising in time. external fct, out real aux(8,2) common /a/ dt, k, n, g0, om, B0, B, a real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/) real :: u(2) = (/0., 0.5/) real :: du(2) = (/0.5, 0.5/) integer :: nd = 2 real :: dt = 2*3.14159265 real :: t, m0, B, B0, a real :: f0 = 0.00 f = 2.45e09 pi = 3.14159265 ez = 4.8e-10 m0 = 9.1e-28 c = 3.0e10 om = 2*pi*f E = 2. w = 1000. a = 0.09 g0 = ez*E/(m0*c*om) B0 = m0*c*om/ez B = B0*(1+a*t) b = B/B0 u(1) = pi u(2) = w/511000.+1 write (*,*) g0 pause n = 1 k = 1 open (unit=2, file = 'res_autorer.dat', status='unknown') call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=' ,ih end subroutine fct (t, u, du) real u(2), du(2), lam common /a/ dt, k, n, g0, om, B, B0, a du(1) = (a*t+1-u(2))/u(2)+g0*sqrt(1./(u(2)**2-1))*sin(u(1)) !equation for electron phase !in SGA du(2) = -g0*sqrt(1.-1./u(2)**2)*cos(u(1)) !equation for electron total energy in SGA !write (*,*) du(1), du(2) !pause return end subroutine out (t, u, du, ih, nd, pr) real u(2), du(2), pr(5), lam common /a/ dt, k, n, g0, om, B, B0, a if (t.ge.k*(dt/2.)) then write (2,1) t,u(1),(u(2)-1)*511. k = k + 1 else end if 1 format (3e10.3) return end !plot the dependencies of the phase and the time, between the energy and the time SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
17. Determination the dependencies between the electron momentum and time, electron momentum and coordinate
! Using the same condition like in program 16 ! metod Runge-Kuta. external fct, out real aux(8,4) common /a/ dt, k, n, g0, om, B0, B, a, y, rl real :: pr(5)=(/0., 2*3.14159*200., 2*3.14159/50, .0001, .0/) real :: u(4) = (/0., 0., 0., 0./) real :: du(4) = (/0.25, 0.25, 0.25, 0.25/) integer :: nd = 4 real :: dt = 2*3.14159265 real :: t, m0, B, B0, a real :: f0 = 0.00 f = 2.45e09 pi = 3.14159265 ez = 4.8e-10 m0 = 9.1e-28 c = 3.0e10 om = 2*pi*f E0 = 5. E = E0 a = 0.00034 B0 = m0*c*om/ez g0 = E0/B0 B = B0 b = B/B0 rl = c/om u(1) = 0. u(2) = 0. u(3) = 0. u(4) = 0. write (*,*) g0, B0, b, E, rl pause n = 1 k = 1 open (unit=2, file = 'res_autores.dat', status='unknown') call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=' ,ih end subroutine fct (t, u, du) real u(4), du(4), lam common /a/ dt, k, n, g0, om, B, B0, a, y, rl y=sqrt(1+u(1)**2+u(2)**2) du(1) = -g0*cos(t)- u(2)*(1+a*t)/y du(2) = -g0*sin(t)+u(1)*(1+a*t)/y du(3) = u(1)/y du(4) = u(2)/y !write (*,*) du(1), du(2), du(3), du(4) !pause return end subroutine out (t, u, du, ih, nd, pr) real u(4), du(4), pr(5), lam common /a/ dt, k, n, g0, om, B, B0, a, y, rl if (t.ge.k*(dt)) then write (2,1) t,(y-1)*511. ,u(3)*rl,u(4)*rl k = k + 1 else end if 1 format (5e12.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(4),DERY(4),AUX(8,4),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
18. Electron acceleration in a uniform electric field using speed light
! Using metod Runge-Kuta external fct, out real aux(8,4) common /a/ dt, k, n, vm, A, om, rl real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/) real :: u(3) = (/0., 0.5/) real :: du(3) = (/0.5, 0.5/) real :: u(4) = (/0., 0.5/) real :: du(4) = (/0.5, 0.5/) integer :: nd = 2 real :: dt = 2*3.14159265 ! real :: lam = 0.1 real :: t real :: c = 3.0e10 B0 = 200.0 E0 = 5000 eq = 4.8e-10 me = 9.1e-28 om = (eq*B0)/(me*c) A = 2.0 vm = A*om n = 1 k = 1 b = B0/B0 rl = c/om u(1) = vm/c u(2) = vm/c u(3) = x/rl u(4) = y/rl open (unit=2, file = '1_rk.dat', status='unknown') call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=', ih end subroutine fct (t, u, du) real u(3), du(3), u(4), du(4) common /a/ dt, k, n, vm, A, om, rl du(1) = -u(2) du(2) = u(1) du(3) = u(1) du(4) = u(2) return end subroutine out (t, u, du, ih, nd, pr) real u(3), du(3), u(4), du(4), pr(5), lam common /a/ dt, k, n, om, vm, A if (t.ge.k*(dt/20.)) then write (2,1) t, u(1), u(2), u(3), u(4) ! k = k + 1 else end if 1 format (f8.3, f8.3, f8.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
19. Charged particle motion in magnetic mirror trap under ECR
! Using Boris method get the equations of the charge particle motion. And analyze the ECR !(electron cyclotron resonance) phenomenon in parabolic approximation of the magnetic !field (mirror trap) for different input values real :: ez = 4.8e-10, c=3.0e10, me=9.1e-28, f=2.4e09, pi=3.14159265 common/a/ gx(250),gy(250),x,y,z,ux,uy,uz,rl,gm,tp,pi,dt,cl2,al,dl,d open (unit=1, file = 'input_t.dat') open (unit=2, file = 'res.dat') ! ! call cpu_time(start) ! read(1,*) kt,kd,E,al,R,B00,x,y,z,W0,fi write(2,1) kt,kd,E,al,R,B00,x,y,z,W0,fi 1 format('kt='I7,2x,'kd=',I3,2x,'E(kV/cm)=',f4.1,2x,'al=',e10.3,2x, & 'R=',f5.3,2x,'B00(%)=',f8.5/'x(cm) =',f4.1,2x,'y(cm) =',f4.1,2x, & 'z(cm)=',f4.1,2x,'W0(eV)=',e10.3,2x,'fi=',f5.1) 2 format('g0=',e10.3,2x,'B0=',f6.1,2x,'cl=',e10.3,2x,'rl =',f6.2/ & 'wmax=',e10.3,2x,'v0=',e10.3,2x,'w00=',e10.3,2x,'rc=',e10.3) pi2=2.0*pi DL=8.0 D=12.0 fi=fi/180.*pi v0=sqrt(1.-1./(1.+w0/511000.)**2) write(*,*) v0*3.0e8 U0=v0/sqrt(1.0-v0**2) W00=511000.*(sqrt(1.0+u0**2)-1.0) om=2.0*pi*f rl=c/om zm=dl/2.0 cl=(zm/rl)/sqrt(R-1.0) cl2=cl*cl zm=zm/rl E=E*1.0e05/3.0e04 g0=ez*E/(me*c*om) B0=me*c*om/ez wmax=511.0*2.0*g0**(2.0/3.0) rc=v0*3.0e10/om ! B0=875 ! B=875 write(2,2) g0, B0, cl, rl, wmax, v0, w00, rc B=B/B0 x=x/rl y=y/rl z=z/rl ux=u0*sin(fi) uy=u0*cos(fi) uz=0. km=0. kp=0 m=250 dt=2.*pi/m do i=1,m gx(i)=-g0*cos(pi2*(i-1)/m)*dt/2.0 gy(i)=-g0*sin(pi2*(i-1)/m)*dt/2.0 end do write(2,*) 'wmax=', 2.*g0**(2./3.)*511. wm=0 do k=1,kt l=k-kp tm=k*dt tp=-(1.+b00)*dt/2. call motion(l) if (k.eq.km+kd) then w=(gm-1.0)*511.0 write(2,3) tm/(2.*pi), x*rl,y*rl,z*rl, w km=km+kd else end if if (k.eq.kp+m) then kp=kp+m else end if end do 3 format(7f12.5) call cpu_time(finish) write(2,*) 'cpu-time =', finish-start end subroutine motion(l) common/a/ gx(250),gy(250),x,y,z,ux,uy,uz,rl,gm,tp,pi,dt,cl2,al,dl,d r = sqrt(x*x + y*y) bxp=-x*z/cl2 byp=-y*z/cl2 bzp=1.0+z*z/cl2-r*r/(2.*cl2) ! bxp=0.0 ! byp=0.0 ! bzp=1.0 u=sqrt(ux**2+uy**2) !Boris's scheme of partical motion uxm = ux + gx(l) uym = uy + gy(l) uzm = uz gmn = sqrt(1 + uxm*uxm + uym*uym + uzm*uzm) tg = tp/gmn tx = tg * bxp ty = tg * byp tz = tg * bzp uxr = uxm + uym*tz - uzm*ty uyr = uym + uzm*tx - uxm*tz uzr = uzm + uxm*ty - uym*tx gmn = sqrt(1 + uxm*uxm + uym*uym + uzm*uzm) txyzr= 1.+ tx*tx + ty*ty + tz*tz sx = 2*tx/txyzr sy = 2*ty/txyzr sz = 2*tz/txyzr uxp = uxm + uyr*sz - uzr*sy uyp = uym + uzr*sx - uxr*sz uzp = uzm + uxr*sy - uyr*sx ux = uxp + gx(l) uy = uyp + gy(l) uz = uzp gm = sqrt (1. + ux**2 + uy**2 + uz**2) gt=dt/gm x = x + ux*gt y = y + uy*gt z = z + uz*gt if (abs(z*rl).ge.dl/2.0.or.r*rl.ge.d/2.0) then write(*,*) 'out',x*rl,y*rl,z*rl stop else endif return end
20. Solving nonlinear equation
!using method Newton-Raphson program heat1 implicit none ! Variables real a,b,c,d,dx,dt,x_i,t_j integer n,m,i,j real, allocatable:: U(:,:) open(10,file='result_out1.txt') a=1 b=4 c=1 d=15 write(0,*) "vvedite melkost razbijenija po x" read(*,*) n ! do i=1,n write(0,*) "vvedite melkost razbijenija po t" read(*,*) m dx=(b-a)/n dt=(d-c)/m allocate(U(n,m)) do j=1,m t_j=j*dt U(1,j)=3 U(n,j)=1 end do do i=2,n-1 x_i=i*dx U(i,1)=1 end do do j=1,m-1 do i=2,n-1 U(i,j+1)=U(i,j)+0.005*(dt/(dx)**2)*(U(i+1,j)-2*U(i,j)+U(i-1,j)) print *,i,j+1,U(i,j+1) end do end do !print *, x_i !end do !do j=1,m !print *, t_j do i=1,n write(10,*) i*dx,U(i,m/100),U(i,m/2),U(i,m) enddo end program heat1
21. Poisson equations solver
! Using swap method in one dimensional case, determine the electric filed of the ! electron layer ! Get the plots of the dependencies between the density and the electric filed with ! coordinate z program poisn1 real d,n0 ! j - number of grid points, al - cr - boudary coefficients ! q(1000),fi(1000) - charge and potentials at grid points ! common/psn/j,al,bl,cl,ar,br,cr,q(1000),fi(1000) ! ! Ez(1000) - electric field strength at grid points ! dimension Ez(1000) ! ! open files: p_in.dat - initial data; in-dstr.dat - initial charge distribution ! fi_Ez.dat - values of potential and electric field strength at grid points ! open(5,file='p_in.dat',status='old') open(6,file='in-dstr.dat',status='new') open(7,file='fi_Ez.dat',status='new') ! ! number of grid points, other initial parameters ! read(5,1) j,al,bl,cl,ar,br,cr 1 format(5x,i5/(5x,f20.8)) write(6,3) j,al,bl,cl,ar,br,cr 3 format(5x,'j=',i5,//2x,6e10.3//) pi=3.14159265 n0=10e+08 d=0.4 ! thick of a layer jz=400 ! number of filled cells del=d/jz ! step (in meters) !z=2 ! charge dencity rho=n0*abs(sin(2*pi*k/400)) qz=rho ! ! spatial initial distribution ! do k=1,j q(k)=0 if(k.ge.50.and.k.le.350) q(k)=n0*abs(sin(2*pi*k/400)) write(*,*) qz end do ! ! Poisson equation solver ! call pnsolv ! ! Calculation of electric field strength at grid points ! do k=2,j-1 Ez(k)=(fi(k-1)-fi(k+1))/(2.*del) write(7,7) k,q(k),fi(k),Ez(k) ! (V/m) end do ! ! Test (analytic solution) ! ! Ean= write(7,8) Ean ! (V/m) 7 format(3x,I3,3x,3e15.5) 8 format(3x,e15.5) ! stop end ! ! SUBROUTINE PNSOLV FOR POISSON EQUATION SOLVING ! subroutine pnsolv common/psn/j,al,bl,cl,ar,br,cr,q(1000),fi(1000) dimension z(1000),y(1000) ! ! (calculation z and y coefficients) ! j1=j+1 jm1=j-1 cm=ar+br z(j)=ar/cm y(j)=cr/cm do i=1,jm1 m=j1-i cm=z(m)-2. z(m-1)=-1./cm y(m-1)=(q(m-1)-y(m))/cm end do ! ! (calculation fi(j) at grid points) ! f0=(al*y(1)-cl)/(al-bl-al*z(1)) fi(1)=z(1)*f0+y(1) do i=1,jm1 fi(i+1)=fi(i)*z(i+1)+y(i+1) end do return end
22. Distribution of Langmuir waves in a uniform plasma
! ! 1D plasma simulation (demo version) ! using numerical method created by particle cell method ! calculate the wavelength and the velocity of wave propogation program plasma_1D real n,m0 common/psn/ nze,j,al,bl,cl,ar,br,cr,q(5000),f(5000),ze(100001) common /E/ g0e,dt,rn,Uz(100001),Ez(5000) common/i/ nzi,qi(5000),zi(100001),Uzi(100001) open(4,file='E(z).dat',status='new') open(7,file='E(t).dat',status='new') open(5,file='P-input.dat',status='old') open(6,file='res.dat',status='new') open(unit=2, file='input.dat') call cpu_time(start) read(2,*) kt, dt, n, cfe ! read(5,1) j,al,bl,cl,ar,br,cr 1 format(5x,i5/(5x,f20.8)) write(6,3) j,al,bl,cl,ar,br,cr 3 format(5x,'j=',i5,//2x,6e10.3//) ! ! Input data ! pi=3.14159265 e=1.6E-19 m0=9.1E-31 c=3.0E8 om=sqrt(n*e*e/(m0*8.85e-12)) write(*,*) 'omega', om Tosc=2.0*pi/om rn=c/om dt=dt*2.0*pi d=0.1e-1 ! thick of a layer jz=4000 ! number of filled cells del=d/jz ! step (in meters) rho=-1.0e16 ! charge dencity qz=rho*1.6e-19*(del**2)/8.85e-12 write(*,*) 'qz =', qz ! qz=0.0e-7 do i=1,4000 q(i)=0. ! Charge of electrons at grid points qi(i)=0. ! Charge of ions at grid points f(i)=0. ! Potential at grid points end do ! call ingen_ue ! call subroutine for initial space electrons distribution call ingen_ui ! call subroutine for initial space ions distribution 7 format(3x,I4,3x,e12.5) 8 format(//) ! do i=1,4000 q(i)=(q(i)-qi(i))*qz ! Charge at grid points (electrons+ions) end do ! open(11,file='shift.dat',status='new') ! do i=1,4000 ! write(11,*) i,q(i),qi(i) ! end do ! close(11) ! stop ! write(6,*) om,Tosc ! do i=1,nze uz(i)=0. ! initial velocities of electrons end do do i=1,nzi uzi(i)=0. ! initial velocities of ions end do kp0=50 kp=50 ! ! Cycle in time ! do k=1,kt t=k*dt call pnsolv ! call Poisson's solver ! do i=2,3999 Ez(i)=(f(i-1)-f(i+1)) !/(2.0*del) ! Calculation of electric field strength ! at grid points end do open(11,file='shift.dat',status='new') do i=1,4000 write(11,*) i,f(i),Ez(i) end do close(11) stop write(7,*) t/om/1.0e-9,Ez(1750) ! write down time (nanosec) and E-field do i=1,4000 q(i)=0. ! Charge and potential at grid point = 0 f(i)=0. end do ! call move ! Integrating motion equation for electrons ! ! open(11,file='shift.dat',status='new') ! do i=1,4000 ! write(11,*) i,q(i),E(i) ! end do ! close(11) ! stop ! do i=1,4000 q(i)=(q(i)-qi(i))*qz ! Charge at grid points (electrons+ions) end do end do ! ! ! End of time cycle ! ! Diagnostics: q(i), qi(i), Ez(i) ! open(3,file='El_dstr.dat',status='unknown') open(13,file='Ions_dstr.dat',status='unknown') do i=1,4000 write(3,*) (i-1)*0.0005,q(i) end do close(3) do i=1,4000 write(13,*) (i-1)*0.0005,qi(i) end do close(13) do i=1,3999 write(4,*) (i-1)*0.0005,Ez(i) end do call cpu_time(finish) write(6,*) 'calculaton time ',finish-start,' sec' end ! ! SUBROUTINE PNSOLV FOR THE POISSON EQUATION SOLVING (Poisson's solver) ! subroutine pnsolv common/psn/ nze,j,al,bl,cl,ar,br,cr,q(5000),f(5000),ze(100001) dimension c(5000),s(5000) ! ! DOWNWARDS SCANNING ! ! j1=j+1 j1=4001 ! jm1=j-1 jm1=3999 cm=ar+br c(4000)=ar/cm s(4000)=cr/cm do i=1,jm1 m=j1-i cm=c(m)-2. c(m-1)=-1./cm s(m-1)=(q(m-1)-s(m))/cm end do ! ! UPWARDS SCANNING ! f0=(al*s(1)-cl)/(al-bl-al*c(1)) f(1)=c(1)*f0+s(1) do i=1,jm1 f(i+1)=f(i)*c(i+1)+s(i+1) end do return end
23. Proton acceleration by the filed of electron leyer
!Using method Runge-Kutta find the max possible acceleration of the proton (maximum !gradient of the magnetic filed). Plot the graphics dependencies between the power and !time, coordinate difference and time external fct, out common /a/ g0, dt, k, pp2b, gme, gmi, dbdz, pi, rl, d real aux(8,4),u(4),du(4),pr(5),n integer :: nd = 4 real :: pi = 3.14159 real :: ez = 4.8e-10, c=3.0e10, me=9.1e-28, f=2.4e09 open (unit=1, file = 'input.txt') open (unit=2, file = 'res.txt') read(1,*) kt,dt,n,d,B0,W0,dbdz write(2,1) kt,dt,n,d,B0,W0,dbdz 1 format('kt=',I5,2x,'dt=',f4.1,2x,'n(cm-3)=',e10.3,2x,'d(cm) =',f4.1,2x, & 'B(Gs)=',f7.1/'W0(keV)=',e10.3,2x,'dBdz(cm-1) =',e10.3) om=ez*B0/(me*c) rl=c/om E0=2.*pi*n*ez*d g0=ez*E0/(me*om*c) d=d/rl dbdz=dbdz*rl b=B0/B0 dt=2*pi/10.0 write(*,*) 'g0, om, b, dbdz', g0,om,b,dbdz pr(1)=0. pr(2)=20.0*pi*kt pr(3)=2.0*pi/100. pr(4)=1.0e-03 pr(5)=0. v0=sqrt(1.-1./(1.+w0/511.)**2) U0=v0/sqrt(1.0-v0**2) W00=511.*(sqrt(1.0+u0**2)-1.0) gme = sqrt(1.0+u0**2) write(*,*) 'v0, u0, w00, gme', v0, u0,w00,gme write(*,*) 'rl', rl u(1)=0. u(2)=0.000/rl u(3)=0.0 u(4)=-0.04/rl dt=50.0*pi k = 1 do i=1,4 du(i)=0.25 end do call rkgs(pr, u, du, nd, ih, fct, out, aux) write(*,*) 'ih=', ih end subroutine fct (t, u, du) real u(4), du(4) common /a/ g0, dt, k, pp2b, gme, gmi, dbdz, pi, rl, d if(u(2)*rl.gt.210.) then write(*,*) 'dbzu =',1.0-dbdz*u(2),dbdz,u(2) stop else end if if(gme**2 - 1.0 - u(1)**2.le.0.) then write(*,*) 'utr =', gme**2 - 1.0 - u(1)**2,gme,u(1),t stop else end if utr=sqrt(gme**2 - 1.0 - u(1)**2) du(1) = 0.5*utr**2/gme*dbdz ! norm du(2) = u(1)/gme du(3) = g0*(u(2)-u(4))/(d/2.)/1836. du(4) = u(3) if(1.0-dbdz*u(2).le.0.) then write(*,*) 'B is zero?=',1.0-dbdz*u(2) stop else end if return end subroutine out (t, u, du, ih, nd, pr) real u(4), du(4), pr(4) common /a/ g0, dt, k, pp2b, gme, gmi, dbdz, pi, rl, d if (t.ge.k*2*pi/5.0) then write (2,1) t/(2.*pi), u(2)*rl, u(4)*rl, (u(2)-u(4))*rl, & 9.1e-28*1836.*(u(3)*3.0e10)**2/2./(1.6e-12)/1.0e06 k = k + 1 if(abs(u(2)-u(4))*rl.gt.d/2.0*rl) then write(*,*) u(2)*rl, u(4)*rl, d/2.*rl write(*,*) 'out of layer' stop else end if if(u(2)*rl.gt.200.) then write(*,*) u(2)*rl, u(4)*rl, d/2.*rl write(2,*) 'end of acceleration length' stop else end if else end if 1 format (f6.1,2x,6e11.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(4),DERY(4),AUX(8,4),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
24. Electron motion in uniform electric filed
!using method Euler real v,x,dt,t,dv,Fm open(2,file='y2.dat', status='unknown') q1=1.0e-9 q2=1.0e-9 r1=0.1 r2=0.2 d=0.5 t0=d/v0 n=100 pi=3.141593 E0=0.0 v0=0.0 dt=t0/n x=x/d write(*,*) dt open(1,file='xv.dat', status='unknown') Fm=q1/me/(vo**2) print*, Fm do k=1,n t=k*dt v=v+Fm*sin(2.0*pi*x)*dt x=x+v*dt write(1,*) t*t0/1.0e-9, v*v0, x*d, sin(2.0*pi*x) if(x.ge.1.0) then write(1,*) t*t0/1.0e-9, v*vo, x*d stop else end if end do end
