wiki:fortran

Version 7 (modified by ivasileska, 4 years ago) (diff)

--

FORTRAN 95 examples

1. Integral function calculation using trapeze, rectangle and Simpson method

real x,dx,s
a=-0.5
b=0.5
n=100
dx=(b-a)/n 
S1=0
S2=0
S3=0
do i=0,n
       S1=S1+(f(i*dx+a)+f((i+1)*dx+a))/2*dx !integral which is calculate by trapeze method
       S2=S2+ f(i*dx+a)*dx !integral which is calculate with rectangle method
if(mod(i,2)==0) then 
       S3=S3+dx/3*2*f(i*dx+a) !integral which is calculate by Simpson method
   else 
       S3=S3+dx/3*4*f(i*dx+a)
   end if       
             
end do 
print *, 'Trap = ', S1 , 'Rec=', S2 , 'Sim=', S3
pause

contains

function f(x) 
f=1/sqrt(1-x**2) !integral function
end function f
end

2. Integral function calculation using Simspon method

real x,dx,s
a=0
b=1.57
n=100
dx=(b-a)/n
S=0
do i=0,n
    if(mod(i,2)==0) then 
       S=S+dx/3*2*f(i*dx+a) !integral which is calculate by Simpson method 
   else 
       S=S+dx/3*4*f(i*dx+a)
   end if
end do
 
print *, 'Sim = ', S 
pause

contains

function f(x) 
f=sin(x)*sin(2*x)*sin(3*x) !integral function
end function f
end

3. Calculation dependence between integral function and coordinates

! Get a graph of the dependence in Origin, differentiate and integrate the graph and 
!compare the resultants with the integrals which are calculated in the first program
real x,dx, a, b
open(1,file='y2.dat', status='unknown')!the document y2 is in the project folder 
a=-0.5
b=0.5
n=100
dx=(b-a)/n
do i=1,n
    write(1,*) i*dx + a, 1/sqrt(1-(i*dx+a)**2)!integral function
end do
end

4. Calculation dependence between integral function and coordinates

! Get a graph of the dependence in Origin, differentiate and integrate the graph and 
!compare the resultant with the integral which is calculated in the second program
real x,dx, a, b
open(1,file='y2.dat', status='unknown')!the document y2 is in the project folder 
a=0
b=1.14
n=100
dx=(b-a)/n
do i=1,n
    write(1,*) i*dx + a, sin(i*dx+a)*sin(2*i*dx+a)*sin(3*i*dx+a)!integral function
end domethod of 
end

5. Simpson formula for Simpson method of integration

program Simpson

integer n
real y1, y2, y3, x, h, a, s
n=100
h=2/n
do i=1,n
   
        x(i)=-1+h(i-1)
        y1(i)=f1(x(i))
        y2(i)=f2(x(i),h)
        y3(i)=f3(x(i),h)
        
         a=h/6*(y1+4*y2+y3)!Simpson's formula
         s=sum(a)
end do

contains
real function f1(x)
    real x
    f1=x*(5-4*x)**(-0.5) !first function
  end function
real function f2(x,h)
    real x, h
    f2=(x+h/2)*(5-4*(x+h/2))**(-0.5) !second function
  end function
real function f3(x,h)
    real x, h
    f3=(x+h)*(5-4*(x+h))**(-0.5) !third function

 end function
print *, smethodof !Simpson method 
end program

6. Characteristics of the electron motion, which is located between charged coaxial rings

!Between two coaxial rings, charged with q1=q2=1nC and their radius are equal to r1=1cm 
!and r2=2cm, are located at distance 5 cm between each other. At halfway between the rings is
!located electron with initial energy equal to 0. Find the dependencies between the 
!velocity(energy)and the time,and electron coordinate with time. Also plot the phase
!trajectory V(x). And explain the resultants.   
program rings

integer :: n=10000
real(8) :: r1=1.0, r2=2.0, q=3.0, d=5.0, eq=4.8e-10, m=9.1e-28 !the unis are in CGS system
real(8) ::  dt, t, x, vx

t=1.0e-6
dt=t/float(n)
vx=0
x=2.5


open (1,file='rings',status='unknown')

do i=1,n 
   E=q*x/sqrt((x**2+r1**2)**3)- q*(5-x)/sqrt(((5-x)**2+r2**2)**3) !energy
   vx=vx+(eq*E/m)*dt !phase trajectory
   x=x+vx*dt !coordinate
   write(1,*) x, i*dt, vx
end do

end program

7. Electron motion only in magnetic field

!Using the equation of electron motion find the phase trajectory. Electron filed E=0
!Part 1. Using only dependence of B0
integer :: n=1000
real(8) :: dt, t, vm, b, rl, om, x, y, vx, vy, B0
real(8) :: eq=-4.8e-10, me=9.1e-28, c = 3.0e10
!B0=me*c*om/ez 
B0 = 200.0
!E0 = 5000
!A = 1e-30
x=1.0
y=1.0
vx=1000.0
vy=1000.0
a=100
!b = B0/B0
!om = eq*B0/me*c
!vm = A*om
!print *, vm
!rl = c/om
!vx = vm/vm
!vy = vm/vm
!x = x/rl
!y = y/rl
!print *, rl
t = 1.0e-8
dt = t/float(n)

open (1,file='1',status='unknown')

do i=1,n 
     vx = vx+(B0*eq)/(c*me)*vy*dt
     vy = vy-(B0*eq)/(c*me)*vx*dt
     x = x+vx*dt
     y = y+vy*dt
   write(1,*) x, vx
end do

pause
end program

8. Electron motion in electric and magnetic filed.

!Like in program 7.Part 1
integer :: n=1000
real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b
real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0 
x = 1.0
y = 1.0
a = 0.5

vx = 1.0e09
vy = 1.0e09
E0=0.0
B0 = 200.0

b=(2*3.14*6)/(1.0e09)

!B0 = 1e-25
!E0 = 5000
!A = 1e-30
!b = B0/B0
!om = eq*B0/me*c
!vm = A*om
!print *, vm
!rl = c/om
!vx = vm/vm
!vy = vm/vm
!x = x/rl
!y = y/rl
!print *, rl
t = 1.00e-08
dt = t/float(n)


open (1,file='1',status='unknown')

do i=1,n
   
   Bz = B0*(1+a*y) 
   !Bz = B0
   x = x + vx*dt
   y = y + vy*dt

   !vx = vx + (vy*eq*Bz)/(c*me)*dt
   !vy = vy - (vx*eq*Bz)/(c*me)*dt

   
   if(i*dt.gt.b) then 
E0 = 5.0
end if
   
  vx = vx +(eq*E0)/(me)*dt+(vy*eq*Bz)/(c*me)*dt
   vy = vy + (eq*E0)/(me)*dt-(vx*eq*Bz)/(c*me)*dt
  ! x = x + vx*dt
   !y = y + vy*dt

   !print *, vx, vy



  

          write (1,*) x, y 

    end do
 pause  
end program

9. Electron motion in electric and magnetic filed

!Like in program 8. Part 2. the differnece is that here we are using different functions 
!for vx and vy
integer :: n=1000
real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b
real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0 
x = 1.0
y = 1.0
a = 0.5

vx = 1.0e09
vy = 1.0e09
E0 = 0.0
B0 = 200.0

b = (2*3.14)/(1.0e09)

t = 1.00e-08
dt = t/float(n)


open (1,file='1',status='unknown')

do i=1,n
   
   Bz = B0*(1+a*y) 
   !Bz = B0
   x = x + vx*dt
   y = y + vy*dt

   vx = vx + (vy*eq*Bz)/(c*me)*dt
   vy = vy - (vx*eq*Bz)/(c*me)*dt

   
   if(i*dt.gt.b) then 
E0 =  3.0
   end if
   
  vx = vx + (vy*eq*Bz)/(c*me)*dt
  vy = vy + (eq*E0)/(me)*dt  - (vx*eq*Bz)/(c*me)*dt
  x = x - vx*dt
  y = y + vy*dt

   !print *, vx, vy



  

          write (1,*) x, y 

    end do
 pause  
end program

10. Electron motion only in magnetic field

! Like in program 7. Part 2 here we have dependence of magnetic filed in z coordinate
integer :: N = 1500
real(8) :: x = 2.5, y = 0.5, B0 = 500.0, We = 4.0054e-08, b = 2.0, pi = 3.14159, m = 9.1e-28, e = -4.8e-10, c = 2.99e10
real(8) vx, vy, v, dt, t

v = sqrt(2*We/m)
vx = sin(pi/4.0)*v
vy = cos(pi/4.0)*v
t = 5e-9
dt=t/float(N)
!Norma
!tau = t*10e8
!me = m*1e28
!qe = e*1e10
!vl = c*1e-10
!ve = v/c
!vx = sin(pi/4.0)*v
!vy = cos(pi/4.0)*v
!dt = tau/float(N)

open (1,file='mov',status='replace')
do i=1,N
   Bz = B0*(1-b**2*(x**2+y**2))
   x = x + vx*dt
   y = y + vy*dt
   vx = vx + vy/m*Bz*e/c*dt
   vy = vy - vx/m*Bz*e/c*dt
   print *, x, y
   pause
   write(1,*) y, x
end do
end program

11. Determination the proton energy and power

integer :: n=1000
real(8) :: dt, t, x, y, vx, vy, v, E, W
real(8) ::  mp = 1.67e-24, q = 4.8e-10
x = 1.0
y = 1.0
d = 1.0
a = 30
W0 = 1000

v = sqrt(2*W0/mp)


t = 1.00e-05
dt = t/float(n)


open (1,file='1',status='unknown')

do i=1,n
   
  vx = v*sin(a)
  vy = v*cos(a)
  x = x + vx*dt
   y = y + vy*dt
  
  E = (4*3.14*mp*y)/d
  W = (E*q)/ x
   

   

   !print *, vx, vy



  
          write (1,*) W, t

    end do
 pause  
end program

12. Runge-Kutta method for harmonic oscillations

!  metod Runge-Kuta. Oscillator. Part 1


external fct, out
real aux(8,2)
common /a/ b, lam, f0, dt, k, n

real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)
real :: u(2) = (/0., 0.5/) 
real :: du(2) = (/0.5, 0.5/)
real :: b=0.003
integer :: nd = 2
          
real :: dt = 2*3.14159265 
   
! real :: lam = 0.1
real :: f0 = 0.00
 
open (unit=2, file = 'osc_rk.dat', status='new') 

write(*,*) 'x=', u(1), 'v=', u(2)
pause
call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=', ih

stop
end

subroutine fct (t, u, du)

real u(2), du(2), lam 
common /a/ b, lam, f0, dt, k, n

du(1) = -u(2)-2.*b*u(1)                  
du(2) = u(1)
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(2), du(2), pr(5), lam
common /a/ b, lam, f0, dt, k, n

if (t.ge.k*(dt/20.)) then

 write (2,1) t,u(1),u(2)  ! 
  k = k + 1
   else
    end if

1 format (f8.3, f8.3, f8.3)

return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               

13. Runge-Kutta method for harmonic oscillations

!  metod Runge-Kuta. Oscillator. Part 2 (the differnece between part 1 is that here we are 
!  using onother function for du(1))

external fct, out
real aux(8,2)
common /a/ b, lam, f0, dt, k, n, c 

real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/)
real :: u(2) = (/0., 0.5/) 
real :: du(2) = (/0.5, 0.5/)
real :: b=0.003
integer :: nd = 2
real::  c = 2.0
real :: dt = 2*3.14159265 

   
! real :: lam = 0.1
real :: f0 = 0.00

n = 1 
k = 1
 
open (unit=2, file = 'osc_rk.dat', status='new') 

write(*,*) 'x=', u(1), 'v=', u(2)
pause
call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=', ih

stop
end

subroutine fct (t, u, du)

real u(2), du(2), lam 
common /a/ b, lam, f0, dt, k, n, c

du(1) = c*cos(0.5*t)-u(2)-2.*b*u(1)                  
du(2) = u(1)
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(2), du(2), pr(5), lam
common /a/ b, lam, f0, dt, k, n, c

if (t.ge.k*(dt/20.)) then

 write (2,1) t,u(1),u(2)  ! 
  k = k + 1
   else
    end if

1 format (f8.3, f8.3, f8.3)

return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               

14. Runge-Kutta method for harmonic oscillations

!Part 3. Another example and using another function for du(1)
external fct, out
real aux(8,2)
common /a/ dt, k, n, A, Vm, om
real :: pr(5) = (/0.,2*3.14159*10., 2*3.14159/50, .0001, .0/)
real :: du(2) = (/0.5, 0.5/)
integer :: nd = 2
real :: dt=2*3.14159265
real :: t
A = 2.0
om = 1.0
Vm = A*om
u(1) = Vm/Vm
u(2) = 0
k = 1
n=1

open (unit=2, file = 'osc_rk0.dat' ,status='unknown')
write(*,*) 'x=' , u(1), 'v=', u(2)
pause
call rkgs(pr, u, du, nd, ih, fct, out, aux)
write(2,*) 'ih=', ih
stop
end
subroutine fct (t, u, du)
real u(2), du(2)
common /a/ dt, k, n, A, Vm, om
du(1) = -u(2)/(1+0.1*sin(om*dt))
du(2) = u(1)
return
end


subroutine out (t, u, du, ih, nd, pr)
real u(2), du(2), pr(5)
common /a/ dt, k, n, A, Vm, om

if(t.ge.k*(dt/50.)) then
write (2,1) t/om, u(1)*Vm, u(2)*A
k=k+1
else
end if

1 format (f8.3, f8.3, f8.3)
return
end

      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               

15. Electron acceleration in uniform electric field

!  Using metod Runge-Kuta. .

external fct, out
real aux(8,4)
common /a/  dt, k, n, vm, A, om, rl, x, y, r, u(1), u(2), u(3), u(4)

real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)


integer :: nd = 2
real :: dt = 2*3.14159265 
real :: t
real :: c = 3.0e10
!B0=me*c*om/ez 
B0 = 200.0
E0 = 5000
eq = 4.8e-10
me = 9.1e-28
om = eq*B0/me*c
A = 2.0
vm = A*om
n = 1 
k = 1
b = B0/B0
rl = c/om
u(1) = vm/vm
u(2) = vm/vm
u(3) = x/rl
u(4) = y/rl
 
open (unit=2, file = '1_rk.dat', status='unknown') 

write(*,*) 'x=', u(3), 'vx=', u(1), 'y=', u(4), 'vy=', u(2)
call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=', ih

end

subroutine fct (t, u, du)

real u(3), du(3), u(4), du(4)
common /a/ dt, k, n, vm, A, om, rl, x, y, r,  u(1), u(2), u(3), u(4)
  
du(1) = -u(3)
du(2) = u(4)
du(3) = u(1)
du(4) = u(2)              
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)

common /a/ dt, k, n, om, vm, A, rl, x, y, r,  u(1), u(2), u(3), u(4)

if (t.ge.k*(dt/20.)) then

 write (4,1) t, u(1), u(2), u(3), u(4)   
  k = k + 1
   else
    end if

1 format (f8.3, f8.3, f8.3)

return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               


16. Determination the conditions of electron capture in SGA mode and the bunch time

!  Using metod Runge-Kuta. Resonant and auto resonant. 
! SGA (synchrotron gyro-magnetic auto resonant) is a self sustaining ECR plasma in 
!magnetic field which is rising in time. 

external fct, out
real aux(8,2)
common /a/  dt, k, n, g0, om, B0, B, a 

real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/)
real :: u(2) = (/0., 0.5/) 
real :: du(2) = (/0.5, 0.5/)

integer :: nd = 2
real :: dt = 2*3.14159265 

real :: t, m0, B, B0, a
real :: f0 = 0.00
f = 2.45e09
pi = 3.14159265
ez = 4.8e-10
m0 = 9.1e-28
c = 3.0e10
om = 2*pi*f
E = 2.
w = 1000.
a = 0.09
g0 = ez*E/(m0*c*om) 
B0 = m0*c*om/ez
B = B0*(1+a*t)
b = B/B0
u(1) = pi
u(2) = w/511000.+1
write (*,*) g0
pause


 
n = 1 
k = 1 
 
open (unit=2, file = 'res_autorer.dat', status='unknown') 


call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=' ,ih 



end

subroutine fct (t, u, du)

real u(2), du(2), lam 
common /a/ dt, k, n, g0, om, B, B0, a

du(1) = (a*t+1-u(2))/u(2)+g0*sqrt(1./(u(2)**2-1))*sin(u(1)) !equation for electron phase 
                                                            !in SGA       
du(2) = -g0*sqrt(1.-1./u(2)**2)*cos(u(1)) !equation for electron total energy in SGA

!write (*,*)  du(1), du(2)
!pause
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(2), du(2), pr(5), lam
common /a/ dt, k, n, g0, om, B, B0, a

if (t.ge.k*(dt/2.)) then

 write (2,1) t,u(1),(u(2)-1)*511. 
  k = k + 1
   else
    end if

1 format (3e10.3)

return
end 
!plot the dependencies of the phase and the time, between the energy and the time   


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               


17. Determination the dependencies between the electron momentum and time, electron momentum and coordinate

!  Using the same condition like in program 16 
!  metod Runge-Kuta.
external fct, out
real aux(8,4)
common /a/  dt, k, n, g0, om, B0, B, a, y, rl 

real :: pr(5)=(/0., 2*3.14159*200., 2*3.14159/50, .0001, .0/)
real :: u(4) = (/0., 0., 0., 0./) 
real :: du(4) = (/0.25, 0.25, 0.25, 0.25/)

integer :: nd = 4
real :: dt = 2*3.14159265 

real :: t, m0, B, B0, a
real :: f0 = 0.00
f = 2.45e09
pi = 3.14159265
ez = 4.8e-10
m0 = 9.1e-28
c = 3.0e10
om = 2*pi*f
E0 = 5.
E = E0
a = 0.00034
B0 = m0*c*om/ez
g0 = E0/B0
B = B0
b = B/B0
rl = c/om
u(1) = 0.
u(2) = 0.
u(3) = 0.
u(4) = 0.
write (*,*) g0, B0, b, E, rl
pause


 
n = 1 
k = 1 
 
open (unit=2, file = 'res_autores.dat', status='unknown') 


call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=' ,ih 



end

subroutine fct (t, u, du)

real u(4), du(4), lam 
common /a/ dt, k, n, g0, om, B, B0, a, y, rl
y=sqrt(1+u(1)**2+u(2)**2)
du(1) = -g0*cos(t)- u(2)*(1+a*t)/y       
du(2) = -g0*sin(t)+u(1)*(1+a*t)/y
du(3) = u(1)/y
du(4) = u(2)/y

!write (*,*)  du(1), du(2), du(3), du(4)
!pause
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(4), du(4), pr(5), lam
common /a/ dt, k, n, g0, om, B, B0, a, y, rl

if (t.ge.k*(dt)) then

 write (2,1) t,(y-1)*511. ,u(3)*rl,u(4)*rl 
  k = k + 1
   else
    end if

1 format (5e12.3)

return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(4),DERY(4),AUX(8,4),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39      
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               

18. Electron acceleration in a uniform electric field using speed light

! Using metod Runge-Kuta

external fct, out
real aux(8,4)
common /a/  dt, k, n, vm, A, om, rl  

real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)
real :: u(3) = (/0., 0.5/) 
real :: du(3) = (/0.5, 0.5/)
real :: u(4) = (/0., 0.5/) 
real :: du(4) = (/0.5, 0.5/)
integer :: nd = 2
real :: dt = 2*3.14159265 
! real :: lam = 0.1
real :: t
real :: 
c = 3.0e10
B0 = 200.0
E0 = 5000
eq = 4.8e-10
me = 9.1e-28
om = (eq*B0)/(me*c)
A = 2.0
vm = A*om
n = 1 
k = 1
b = B0/B0
rl = c/om
u(1) = vm/c
u(2) = vm/c
u(3) = x/rl
u(4) = y/rl


 
open (unit=2, file = '1_rk.dat', status='unknown') 


call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=', ih


end

subroutine fct (t, u, du)

real u(3), du(3), u(4), du(4)
common /a/ dt, k, n, vm, A, om, rl  

du(1) = -u(2)
du(2) = u(1)
du(3) = u(1)
du(4) = u(2)              
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(3), du(3), u(4), du(4), pr(5), lam
common /a/ dt, k, n, om, vm, A

if (t.ge.k*(dt/20.)) then

 write (2,1) t, u(1), u(2), u(3), u(4)  ! 
  k = k + 1
   else
    end if

1 format (f8.3, f8.3, f8.3)

return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               


19. Charged particle motion in magnetic mirror trap under ECR

! Using Boris method get the equations of the charge particle motion. And analyze the ECR 
!(electron cyclotron resonance) phenomenon in parabolic approximation of the magnetic 
!field (mirror trap) for different input values

real :: ez = 4.8e-10, c=3.0e10, me=9.1e-28, f=2.4e09, pi=3.14159265

common/a/ gx(250),gy(250),x,y,z,ux,uy,uz,rl,gm,tp,pi,dt,cl2,al,dl,d


open (unit=1, file = 'input_t.dat') 
open (unit=2, file = 'res.dat') 
!
!
  call cpu_time(start)
!
read(1,*) kt,kd,E,al,R,B00,x,y,z,W0,fi
write(2,1) kt,kd,E,al,R,B00,x,y,z,W0,fi
1 format('kt='I7,2x,'kd=',I3,2x,'E(kV/cm)=',f4.1,2x,'al=',e10.3,2x,    &
'R=',f5.3,2x,'B00(%)=',f8.5/'x(cm) =',f4.1,2x,'y(cm) =',f4.1,2x,    &
'z(cm)=',f4.1,2x,'W0(eV)=',e10.3,2x,'fi=',f5.1)
2 format('g0=',e10.3,2x,'B0=',f6.1,2x,'cl=',e10.3,2x,'rl =',f6.2/    &
'wmax=',e10.3,2x,'v0=',e10.3,2x,'w00=',e10.3,2x,'rc=',e10.3)

pi2=2.0*pi

DL=8.0
D=12.0
fi=fi/180.*pi

v0=sqrt(1.-1./(1.+w0/511000.)**2)  

write(*,*) v0*3.0e8               
U0=v0/sqrt(1.0-v0**2)              
W00=511000.*(sqrt(1.0+u0**2)-1.0) 
om=2.0*pi*f
rl=c/om                            
zm=dl/2.0                          
 cl=(zm/rl)/sqrt(R-1.0)
 cl2=cl*cl                          
 zm=zm/rl                          

E=E*1.0e05/3.0e04                  
g0=ez*E/(me*c*om)                  
B0=me*c*om/ez                      
wmax=511.0*2.0*g0**(2.0/3.0)        
rc=v0*3.0e10/om                    
 
!    B0=875                       
!    B=875
 write(2,2) g0, B0, cl, rl, wmax, v0, w00, rc
     B=B/B0                     
        x=x/rl           
        y=y/rl
        z=z/rl
        ux=u0*sin(fi)               
        uy=u0*cos(fi)
        uz=0.                       
        km=0.  
        kp=0

        m=250
        dt=2.*pi/m
          

        do i=1,m
        gx(i)=-g0*cos(pi2*(i-1)/m)*dt/2.0       
        gy(i)=-g0*sin(pi2*(i-1)/m)*dt/2.0       
        end do

        write(2,*) 'wmax=', 2.*g0**(2./3.)*511.
        wm=0


        do k=1,kt
        l=k-kp
        tm=k*dt
 
        tp=-(1.+b00)*dt/2.         
        call motion(l)
if (k.eq.km+kd) then  
        w=(gm-1.0)*511.0
    write(2,3) tm/(2.*pi), x*rl,y*rl,z*rl, w

        km=km+kd
        else
        end if
if (k.eq.kp+m) then  
        kp=kp+m
        else
        end if
        end do
3 format(7f12.5)


     call cpu_time(finish)
        write(2,*) 'cpu-time =', finish-start
        end
                                             
subroutine motion(l)    
common/a/ gx(250),gy(250),x,y,z,ux,uy,uz,rl,gm,tp,pi,dt,cl2,al,dl,d



      r = sqrt(x*x + y*y)
      bxp=-x*z/cl2
      byp=-y*z/cl2
      bzp=1.0+z*z/cl2-r*r/(2.*cl2)

!      bxp=0.0
!      byp=0.0
!      bzp=1.0

      u=sqrt(ux**2+uy**2)

!Boris's scheme of partical motion
      uxm = ux  + gx(l) 
      uym = uy  + gy(l)    
      uzm = uz

      gmn = sqrt(1 + uxm*uxm + uym*uym + uzm*uzm) 

      tg = tp/gmn
      tx = tg * bxp
      ty = tg * byp
      tz = tg * bzp

      uxr = uxm + uym*tz - uzm*ty 
      uyr = uym + uzm*tx - uxm*tz
      uzr = uzm + uxm*ty - uym*tx
      gmn = sqrt(1 + uxm*uxm + uym*uym + uzm*uzm) 

      txyzr= 1.+ tx*tx + ty*ty + tz*tz

      sx = 2*tx/txyzr
      sy = 2*ty/txyzr
      sz = 2*tz/txyzr

      uxp = uxm + uyr*sz - uzr*sy 
      uyp = uym + uzr*sx - uxr*sz
      uzp = uzm + uxr*sy - uyr*sx

      ux = uxp + gx(l) 
      uy = uyp + gy(l)     
      uz = uzp

      gm = sqrt (1. + ux**2 + uy**2 + uz**2)
      gt=dt/gm
      x = x + ux*gt
      y = y + uy*gt
      z = z + uz*gt

      if (abs(z*rl).ge.dl/2.0.or.r*rl.ge.d/2.0) then 
                                                      
      write(*,*) 'out',x*rl,y*rl,z*rl
      stop
      else
      endif

      return 
      end


20. Solving nonlinear equation

!using method Newton-Raphson 
    program heat1

    implicit none

    ! Variables
    real a,b,c,d,dx,dt,x_i,t_j
    integer n,m,i,j
    
    real, allocatable:: U(:,:)
    
    open(10,file='result_out1.txt')
    a=1
    b=4
    c=1
    d=15
    write(0,*) "vvedite melkost razbijenija po x"
    read(*,*) n
    
   ! do i=1,n
     
     write(0,*) "vvedite melkost razbijenija po t"
        read(*,*) m
        dx=(b-a)/n
         dt=(d-c)/m
        allocate(U(n,m))
        
        
       do j=1,m     
        t_j=j*dt
        U(1,j)=3
        U(n,j)=1
        end do
        
        do i=2,n-1   
         x_i=i*dx
        U(i,1)=1
        end do
        
        
       
         do j=1,m-1       
         do i=2,n-1       
         U(i,j+1)=U(i,j)+0.005*(dt/(dx)**2)*(U(i+1,j)-2*U(i,j)+U(i-1,j))
         print *,i,j+1,U(i,j+1)
         end do
         end do
   
   
        !print *, x_i
    !end do
   
        !do j=1,m
        
    
   
    !print *, t_j
   
    
    do i=1,n
    write(10,*) i*dx,U(i,m/100),U(i,m/2),U(i,m)
    enddo
    

    end program heat1

21. Poisson equations solver

! Using swap method in one dimensional case, determine the electric filed of the 
! electron layer
! Get the plots of the dependencies between the density and the electric filed with   
! coordinate z 
      program poisn1
real d,n0
!   j - number of grid points, al - cr  - boudary coefficients
!   q(1000),fi(1000) - charge and potentials at grid points
!
      common/psn/j,al,bl,cl,ar,br,cr,q(1000),fi(1000)
!
!     Ez(1000) - electric field strength at grid points 
!
      dimension Ez(1000)
!
!  open files: p_in.dat - initial data;  in-dstr.dat - initial charge distribution 
!  fi_Ez.dat - values of potential and electric field strength at grid points
!
      open(5,file='p_in.dat',status='old')
      open(6,file='in-dstr.dat',status='new')
      open(7,file='fi_Ez.dat',status='new')
!
!                           number of grid points, other initial parameters
!

      read(5,1) j,al,bl,cl,ar,br,cr
    1 format(5x,i5/(5x,f20.8))
      write(6,3) j,al,bl,cl,ar,br,cr
    3 format(5x,'j=',i5,//2x,6e10.3//)
      pi=3.14159265
      n0=10e+08
      d=0.4               !  thick of a layer
      jz=400              !  number of filled cells
      del=d/jz            ! step (in meters)
      !z=2 
     
! charge dencity 
     
      rho=n0*abs(sin(2*pi*k/400))            
       qz=rho
!
!   spatial initial distribution
!
      do k=1,j
      q(k)=0
      if(k.ge.50.and.k.le.350) q(k)=n0*abs(sin(2*pi*k/400)) 
      write(*,*) qz
      end do
!      
!    Poisson equation solver
!      
      call pnsolv
!
!    Calculation of electric field strength at grid points
!
      do k=2,j-1
      Ez(k)=(fi(k-1)-fi(k+1))/(2.*del)
      write(7,7) k,q(k),fi(k),Ez(k)     !   (V/m)
      end do
!
!    Test (analytic solution)
!
!      Ean=               
      write(7,8) Ean                    !     (V/m)
    7 format(3x,I3,3x,3e15.5)  
    8 format(3x,e15.5)  

!
      stop
      end
!
!     SUBROUTINE PNSOLV FOR POISSON EQUATION SOLVING
!
      subroutine pnsolv
      common/psn/j,al,bl,cl,ar,br,cr,q(1000),fi(1000)
      dimension z(1000),y(1000)
!
!     (calculation  z and y coefficients)
!  
      j1=j+1
      jm1=j-1              
      cm=ar+br
      z(j)=ar/cm
      y(j)=cr/cm
      do i=1,jm1
          m=j1-i
          cm=z(m)-2.
          z(m-1)=-1./cm
          y(m-1)=(q(m-1)-y(m))/cm
      end do
!
!      (calculation  fi(j)  at grid points) 
!
      f0=(al*y(1)-cl)/(al-bl-al*z(1))
      fi(1)=z(1)*f0+y(1)
      do i=1,jm1
          fi(i+1)=fi(i)*z(i+1)+y(i+1)
      end do
      return
      end

22. Distribution of Langmuir waves in a uniform plasma

!
! 1D plasma simulation (demo version)
! using numerical method created by particle cell method
! calculate the wavelength and the velocity of wave propogation

 program plasma_1D
   real n,m0
      common/psn/ nze,j,al,bl,cl,ar,br,cr,q(5000),f(5000),ze(100001)
      common /E/ g0e,dt,rn,Uz(100001),Ez(5000)      
      common/i/ nzi,qi(5000),zi(100001),Uzi(100001)
  open(4,file='E(z).dat',status='new')
  open(7,file='E(t).dat',status='new')
  open(5,file='P-input.dat',status='old')
  open(6,file='res.dat',status='new')
  open(unit=2, file='input.dat')
                                 call cpu_time(start)
      read(2,*)  kt, dt, n, cfe        !  
      read(5,1) j,al,bl,cl,ar,br,cr
    1 format(5x,i5/(5x,f20.8))
      write(6,3) j,al,bl,cl,ar,br,cr
    3 format(5x,'j=',i5,//2x,6e10.3//)
!
!  Input data
!
      pi=3.14159265
      e=1.6E-19
      m0=9.1E-31
      c=3.0E8
      om=sqrt(n*e*e/(m0*8.85e-12))
      write(*,*) 'omega', om
      Tosc=2.0*pi/om
      rn=c/om
      dt=dt*2.0*pi

      d=0.1e-1            !   thick of a layer
      jz=4000             !   number of filled cells
      del=d/jz            !   step (in meters)

      rho=-1.0e16         ! charge dencity
      qz=rho*1.6e-19*(del**2)/8.85e-12
      write(*,*) 'qz =', qz

!      qz=0.0e-7

      do i=1,4000
      q(i)=0.         ! Charge of electrons at grid points
      qi(i)=0.         ! Charge of ions at grid points
      f(i)=0.           ! Potential at grid points
      end do
!
      call ingen_ue   ! call subroutine for initial space electrons distribution 
      call ingen_ui   ! call subroutine for initial space ions distribution 

    7 format(3x,I4,3x,e12.5)
    8 format(//)
!
      do i=1,4000
      q(i)=(q(i)-qi(i))*qz   ! Charge at grid points (electrons+ions)
      end do

!      open(11,file='shift.dat',status='new')
!      do i=1,4000
!      write(11,*) i,q(i),qi(i)
!      end do
!      close(11)
!      stop



!
      write(6,*) om,Tosc
!
      do i=1,nze
      uz(i)=0.          ! initial velocities of electrons
      end do
      do i=1,nzi
      uzi(i)=0.         ! initial velocities of ions
      end do
      kp0=50
      kp=50
!
!     Cycle in time
!
      do k=1,kt
      t=k*dt
      call pnsolv        ! call Poisson's solver
!
      do i=2,3999
      Ez(i)=(f(i-1)-f(i+1)) !/(2.0*del)   ! Calculation of electric field strength
                                     ! at grid points
      end do

      open(11,file='shift.dat',status='new')
      do i=1,4000
      write(11,*) i,f(i),Ez(i)
      end do
      close(11)
      stop


      write(7,*) t/om/1.0e-9,Ez(1750)   ! write down time (nanosec) and E-field 
      do i=1,4000
      q(i)=0.            ! Charge and potential at grid point = 0
      f(i)=0.
      end do
!
      call move         ! Integrating motion equation for electrons 
!
!      open(11,file='shift.dat',status='new')
!      do i=1,4000
!      write(11,*) i,q(i),E(i)
!      end do
!      close(11)
!      stop

!
      do i=1,4000
      q(i)=(q(i)-qi(i))*qz   ! Charge at grid points (electrons+ions)
      end do
      end do
!
!                ! End of time cycle
!
!                Diagnostics: q(i), qi(i), Ez(i) 
!
      open(3,file='El_dstr.dat',status='unknown')
      open(13,file='Ions_dstr.dat',status='unknown')
      do i=1,4000
      write(3,*) (i-1)*0.0005,q(i)
      end do
      close(3)
      do i=1,4000
      write(13,*) (i-1)*0.0005,qi(i)
      end do
      close(13)
      do i=1,3999
      write(4,*) (i-1)*0.0005,Ez(i)
      end do
             call cpu_time(finish)
write(6,*) 'calculaton time ',finish-start,' sec'  
      end
!
!     SUBROUTINE PNSOLV FOR THE POISSON EQUATION SOLVING (Poisson's solver)
!
      subroutine pnsolv
      common/psn/ nze,j,al,bl,cl,ar,br,cr,q(5000),f(5000),ze(100001)
      dimension c(5000),s(5000)

!
!            DOWNWARDS SCANNING
!  
!      j1=j+1
      j1=4001
!      jm1=j-1 
      jm1=3999             
      cm=ar+br
      c(4000)=ar/cm
      s(4000)=cr/cm
      do i=1,jm1
          m=j1-i
          cm=c(m)-2.
          c(m-1)=-1./cm
          s(m-1)=(q(m-1)-s(m))/cm
      end do
!
!            UPWARDS SCANNING
!
      f0=(al*s(1)-cl)/(al-bl-al*c(1))
      f(1)=c(1)*f0+s(1)
      do i=1,jm1
          f(i+1)=f(i)*c(i+1)+s(i+1)
      end do
      return
      end

23. Proton acceleration by the filed of electron leyer


!Using method Runge-Kutta find the max possible acceleration of the proton (maximum 
!gradient of the magnetic filed). Plot the graphics dependencies between the power and 
!time, coordinate difference and time  


external fct, out
common /a/ g0, dt, k, pp2b, gme, gmi, dbdz, pi, rl, d
real aux(8,4),u(4),du(4),pr(5),n 
integer :: nd = 4
real :: pi = 3.14159  
real :: ez = 4.8e-10, c=3.0e10, me=9.1e-28, f=2.4e09

open (unit=1, file = 'input.txt') 
open (unit=2, file = 'res.txt') 
read(1,*) kt,dt,n,d,B0,W0,dbdz
write(2,1) kt,dt,n,d,B0,W0,dbdz
1 format('kt=',I5,2x,'dt=',f4.1,2x,'n(cm-3)=',e10.3,2x,'d(cm) =',f4.1,2x, &
'B(Gs)=',f7.1/'W0(keV)=',e10.3,2x,'dBdz(cm-1) =',e10.3)

   om=ez*B0/(me*c)
     rl=c/om
      E0=2.*pi*n*ez*d
       g0=ez*E0/(me*om*c)
        d=d/rl
      dbdz=dbdz*rl   
       b=B0/B0
       dt=2*pi/10.0
write(*,*) 'g0, om, b, dbdz', g0,om,b,dbdz
pr(1)=0. 
pr(2)=20.0*pi*kt 
pr(3)=2.0*pi/100.
pr(4)=1.0e-03
pr(5)=0.

v0=sqrt(1.-1./(1.+w0/511.)**2)
U0=v0/sqrt(1.0-v0**2)
W00=511.*(sqrt(1.0+u0**2)-1.0)


gme = sqrt(1.0+u0**2)  

write(*,*) 'v0, u0, w00, gme', v0, u0,w00,gme
write(*,*) 'rl', rl

u(1)=0.
u(2)=0.000/rl
u(3)=0.0
u(4)=-0.04/rl
dt=50.0*pi

k = 1
do i=1,4
du(i)=0.25
end do

call rkgs(pr, u, du, nd, ih, fct, out, aux)
write(*,*) 'ih=', ih
end

subroutine fct (t, u, du)
real u(4), du(4) 
common /a/ g0, dt, k, pp2b, gme, gmi, dbdz, pi, rl, d  
if(u(2)*rl.gt.210.) then
write(*,*) 'dbzu =',1.0-dbdz*u(2),dbdz,u(2)  
stop
else
end if

if(gme**2 - 1.0 - u(1)**2.le.0.) then
write(*,*) 'utr =', gme**2 - 1.0 - u(1)**2,gme,u(1),t
stop
else
end if
 utr=sqrt(gme**2 - 1.0 - u(1)**2)

            du(1) = 0.5*utr**2/gme*dbdz                 ! norm
             du(2) = u(1)/gme
              du(3) = g0*(u(2)-u(4))/(d/2.)/1836.
               du(4) = u(3)
if(1.0-dbdz*u(2).le.0.) then
write(*,*) 'B is zero?=',1.0-dbdz*u(2) 
stop
else
end if
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(4), du(4), pr(4)
common /a/ g0, dt, k, pp2b, gme, gmi, dbdz, pi, rl, d
if (t.ge.k*2*pi/5.0) then
 write (2,1) t/(2.*pi), u(2)*rl, u(4)*rl, (u(2)-u(4))*rl,  &
             9.1e-28*1836.*(u(3)*3.0e10)**2/2./(1.6e-12)/1.0e06  
k = k + 1

if(abs(u(2)-u(4))*rl.gt.d/2.0*rl) then
write(*,*) u(2)*rl, u(4)*rl, d/2.*rl
write(*,*) 'out of layer'
stop
else
end if

 if(u(2)*rl.gt.200.) then
write(*,*) u(2)*rl, u(4)*rl, d/2.*rl
write(2,*) 'end of acceleration length'
stop
else
end if

else
end if

1 format (f6.1,2x,6e11.3)
return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(4),DERY(4),AUX(8,4),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                              
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               

24. Electron motion in uniform electric filed

!using method Euler
real v,x,dt,t,dv,Fm
open(2,file='y2.dat', status='unknown')
q1=1.0e-9
q2=1.0e-9
r1=0.1
r2=0.2
d=0.5
t0=d/v0
n=100
pi=3.141593
E0=0.0
v0=0.0
dt=t0/n
x=x/d
write(*,*) dt
open(1,file='xv.dat', status='unknown')
Fm=q1/me/(vo**2)
print*, Fm
do k=1,n
    t=k*dt
    v=v+Fm*sin(2.0*pi*x)*dt
    x=x+v*dt
    
write(1,*) t*t0/1.0e-9, v*v0, x*d, sin(2.0*pi*x)
if(x.ge.1.0) then
write(1,*) t*t0/1.0e-9, v*vo, x*d
stop
else
end if
end do
end