wiki:fortran

Version 12 (modified by Leon Kos, 7 years ago) (diff)

Move comments out

FORTRAN 95 examples

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 input
!values 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 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

!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