| Version 16 (modified by , 8 years ago) (diff) |
|---|
-
The following FORTRAN examples are showing various computational and …
- Suppose that a charge particle move in magnetic filed. The force from …
- 1. Integral function calculation using trapeze, rectangle and Simpson …
- 2. Integral function calculation using Simpson 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 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 layer
- 24. Electron motion in uniform electric filed
The following FORTRAN examples are showing various computational and numerical methods for mathematical modeling in the field of plasma physics.
Suppose that a charge particle move in magnetic filed. The force from the magnetic filed is always perpendicular to its direction of motion, for that reason the motion is circular. That force is called Lorentz force and depends on the particle velocity and its charge. But if on the charge particle motion is affect plus electric filed, the motion will be treated as the superposition of a relatively fast circular motion around a point called the guiding center and a relatively slow drift of this point.
1. Integral function calculation using trapeze, rectangle and Simpson method
There are three methods which are using for calculation definite integrals in Fortran. First is rectangle method, where the function is integrated using segment integration and like a resultant we obtain rectangle formula. Second is trapeze method and the same like in the first method after integration we obtain the trapeze formula. And the third method is Simpson method where the function which is integrated, is parabolic and after integration we obtain Simpson's formula. In all methods are used segment integration, which means that interval of integration a and b is divided into N segments. The splitting points xn are called nodes. In this program we calculate the same integral but using different methods the get the difference between each 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 Simpson method
Calculate an integral using only Simpson method. Look at program 1 for introduction
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
To compare the resultants which are obtained in the program 1 with a graphical integration and to get the precision of integration, in this program first is obtained the graph of the integral dependence in Origin,and after that it is compared with the resultants in program 1
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
The introduction is the same like in program 3. But here the resultant of integration is compare with program 2.
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
To understand the algorithm of how is working Simpson method, in this program is obtain how with Simpson formula and using the functions f1, f2 and f3 can be integrated a definite integral
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
To main characteristics if the electron motion are the coordinate and the velocity. This program it is a solution of a task which is written below. All units for the inputvalues are in CGS system. 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
The main idea in plasma physics is to understand the motion of the charge particles. For that reason first we will make a solution of the electron motion but only in magnetic filed. Solving this problem by using the equation of electron motion can be obtain the phase trajectory. Electron filed E0=0 Part 1. Using only dependence of B0 or initial magnetic filed
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 magnetic filed.
The introduction is like in program 7. The difference is to solve the problem how will be the electron motion if the velocities vx and vy depends of the magnetic filed Bz which is a function of the coordinate y. E0=0. Part 2
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
To make a conclusion how is the electric motion looks like using plus electric filed. And to get a picture of the motion in the plasma in this program we are using different functions of dependencies for vx and vy. And the electric filed is E0=3. The particle begins to drift.
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
The introduction is like in program 8, but the difference is that here we have trigonometrical dependencies of the initial velocities. Part 3
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
Also for one of the main characteristics in the plasma is the energy and the power of the charge particles. In this program are obtained the energy and the power of the protons in the plasma. Using the basics equations for it.
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
To solve a partial differential equations need to satisfy the condition y(x0)=y0 which is called Cauchy task. The most effective and the commonly used method of the solution of Cauchy task is Runge-Kuta method. It is based on a approximation of the function y, which is obtained by expanding the function in Taylor series. Using method Runge-Kuta can be obtain the basic differential equation for 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
Using method Runge-Kuta obtain the differential equations for oscillator, the same like in program 12 but the difference is that here we are using onother function for du(1) which is velocity. Part 2
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
The same like in program 12 and 12 using Runge Kutta obtain the differential equations but with another example and using !another function for du(1).Part 3
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
The method Runge-Kuta is not used only for obtaining the differential equation for oscilloscope. The same method can be used in different tasks where the differential equations need be solve. Using method Runge-Kuta, it is easy to find the electron acceleration in electric filed.
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
Also the method Runge Kuta is used for solving more complex problems in plasma physics, like SGA (synchrotron gyro-magnetic auto resonant) which is a self sustaining ECR plasma in magnetic field where the time is rising up. First the method Runge-Kuta is used for auto resonant condition and after that using the equations of plasma physics for the SGA mode it are determined the conditions.
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. Using method Runge-Kuta for determination the conditions of the electron motion.
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
The same introductions like in programs 16 and 17 but the difference is that here we have dependence of speed light. To get the relative conditions for the motion. 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
One of the most important thing in plasma physics is to get the motion equations of the particles. To get that equations it is used Boris method. This method is divided in for sections: using the momentum of the particles, their rotations, using the electrical field and at the end to determinate the coordinates. Using this method in this program is analyzed 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
To find a solutions for solving nonlinear equation first need to find the roots in the initial interval. There are a lot of methods for how can these roots be obtained. In this program we are using method Newton-Raphson. This method is more rapidly convergent. To find the roots in the interval we are choosing a number which will serve as an initial approximation.
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,*) "enter the fineness of the partition x" read(*,*) n ! do i=1,n write(0,*) "enter the fineness of the partition 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
To obtain the one dimensional Poisson equation is used the swap method and after the solving the Poisson equation can be determinate 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 propagation
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 layer
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
To integrate the ordinary differential equations is used method Euler. This method is applicable to a system of first order differential equations. In this program this method will help to find the velocity of the electron in uniform electric filed.
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
