wiki:fortran

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

--

FORTRAN 95 examples

Integral function calculation using trapeze, rectangle and Simson method

real x,dx,s
!open(1,file='y1.dat', status='unknown')
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 Simson 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

real x,dx,s
!open(1,file='y1.dat', status='unknown')
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)
   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)
end function f
end

Console3.f90

real x,dx, a, b
open(1,file='y2.dat', status='unknown')
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)
end do
end

Console4.f90

real x,dx, a, b
open(1,file='y2.dat', status='unknown')
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)
end domethod of 
end

Console5.f90

program Simp

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)
         s=sum(a)
end do

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

 end function
print *, smethod of 
end program

Console6.f90

program kolco

integer :: n=1000
real :: r1=1.0, r2=2.0, q=3.0, d=5.0, eq=-4.8e-10, m=9.1e-28 ! ������� ��������� � ���-�������
real ::  dt, t, x, vx

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


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

do i=1,n 
   E=q*x/(x**2+r1**2)**(1.5)- q*(5-x)/((5-x)**2+r2**2)**(1.5)
   vx=vx+(eq*E/m)*dt
   x=x+vx*dt
   write(1,*) vx, x, i*dt
end do

end program

Console7.f90

program kolco

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 ! ������� ��������� � ���-�������
real(8) ::  dt, t, x, vx

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


open (1,file='kolco',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)
   vx=vx+(eq*E/m)*dt
   x=x+vx*dt
   write(1,*) x, i*dt, vx
end do

end program

Console8.f90

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



Console9.f90

!  Console9.f90 
!
!  FUNCTIONS:
!  Console9 - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Console9
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

    program Console9

    implicit none

    ! Variables

    ! Body of Console9
    print *, 'Hello World'

    end program Console9

Console10.f90

!  Osc_rk_0.f90 
!
!  Lab #3
!  metod Runge-Kuta. Oscillator.

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                                                               


Console11.f90

!  Console11.f90 
!
!  FUNCTIONS:
!  Console11 - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Console11
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

    program Console11

    implicit none

    ! Variables

    ! Body of Console11
    print *, 'Hello World'

    end program Console11

Console12.f90

!  Console12.f90 
!
!  FUNCTIONS:
!  Console12 - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Console12
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

    program Console12

    implicit none

    ! Variables

    ! Body of Console12
    print *, 'Hello World'

    end program Console12

Console13.f90

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

Console15.f90

x = (/1,5/)
print *, x
end

Console16.f90

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

Console17.f90

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

Console18.f90

!  Osc_rk_0.f90 
!
!  Lab #3
!  metod Runge-Kuta. Oscillator.

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

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

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                                                               

Console19.f90

!  Osc_rk_0.f90 
!
!  Lab #3
!  metod Runge-Kuta. Oscillator.

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                                                               

Console20.f90

!  Console20.f90 
!
!  FUNCTIONS:
!  Console20 - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Console20
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

    program Console20

    implicit none

    ! Variables

    ! Body of Console20
    print *, 'Hello World'

    end program Console20

Console21.f90

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

Console22.f90

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

Console24.f90

!  Lab #1
!  metod Runge-Kuta. Rezonans i avtorezonans.

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 = 'rez_avtorez.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))          
du(2) = -g0*sqrt(1.-1./u(2)**2)*cos(u(1))

!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 


      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                                                               


Console26.f90

!  Lab #1
!  metod Runge-Kuta. Rezonans i avtorezonans.

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 = 'rez_avtorez.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                                                               

Console29.f90

! ! program El in a mirro trap,  parabolic approximation of the magnetic field
! for a mirror trap. ECR

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                ! ��������� �������� ��������� � m/c
U0=v0/sqrt(1.0-v0**2)              ! ��������� ������� ��������� � �������� mc
W00=511000.*(sqrt(1.0+u0**2)-1.0)  ! ���� ��� ��������� �������
om=2.0*pi*f
rl=c/om                            ! �������������� ������ ������������� ��������
zm=dl/2.0                          ! ������� ����������

                                   ! ������ ��������� cl, ������������� �������� ���������� ����
                                   ! �������, ����� �������� �������� ���������� ���������
 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                      ! ����������� �������� �������� ���������� ���� - ECR ��� ��������� � ������ �����
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.                        ! � ����������� Z ��������� ������� ����� ����.
        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   ! ����� ����������� ����� kd ��������� �����
        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

! �������� ��������� ���� ������� ���������� ���� (���� 2)

      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)

! �������� ��� ������������� ���� (���� 3)
      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


Console30.f90

!  Osc_rk_0.f90 
!
!  Lab #3
!  metod Runge-Kuta. Oscillator.

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

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/)
integer :: nd = 2
real :: dt = 2*3.14159265 
! real :: lam = 0.1
real :: t
real :: om=1
A=2.0
vm=A*om
n = 1 
k = 1
u(1)=vm/vm
u(2)=0
 
open (unit=2, file = 'osc_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(2), du(2), lam 
common /a/ dt, k, n, g, vm, A, om

du(1) = -u(2)/(1+1.2*sin(t*2))               
du(2) = u(1)
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(2), du(2), 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)  ! 
  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                                                               


Console31.f90

    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

Console32.f90

!  Osc_rk_0.f90 
!
!  Lab #3
!  metod Runge-Kuta. Oscillator.

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                                                               


Console33.f90 - Poisson solver

      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

Console35.f90 - 1D plasma simulation

!
! 1D plasma simulation (demo version)
!
 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

Console36.f90 - Proton acceleration


!method Runge-Kutta


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                          Console2.f90                     
   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                                                               

Console37.f90 - 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