Changes between Version 6 and Version 7 of fortran


Ignore:
Timestamp:
Oct 18, 2016, 1:42:44 PM (8 years ago)
Author:
ivasileska
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • fortran

    v6 v7  
    22= FORTRAN 95 examples =
    33
    4 == Integral function calculation using trapeze, rectangle and Simson method ==
     4== 1. Integral function calculation using trapeze, rectangle and Simpson method ==
    55{{{
    66#!fortran
    77real x,dx,s
    8 !open(1,file='y1.dat', status='unknown')
    98a=-0.5
    109b=0.5
    … …  
    1817       S2=S2+ f(i*dx+a)*dx !integral which is calculate with rectangle method
    1918if(mod(i,2)==0) then
    20        S3=S3+dx/3*2*f(i*dx+a) !integral which is calculate by Simson method
     19       S3=S3+dx/3*2*f(i*dx+a) !integral which is calculate by Simpson method
    2120   else
    2221       S3=S3+dx/3*4*f(i*dx+a)
    … …  
    3534}}}
    3635
    37 == ==
     36== 2. Integral function calculation using Simspon method ==
    3837{{{
    3938#!fortran
    4039real x,dx,s
    41 !open(1,file='y1.dat', status='unknown')
    4240a=0
    4341b=1.57
    … …  
    4745do i=0,n
    4846    if(mod(i,2)==0) then
    49        S=S+dx/3*2*f(i*dx+a)
     47       S=S+dx/3*2*f(i*dx+a) !integral which is calculate by Simpson method
    5048   else
    5149       S=S+dx/3*4*f(i*dx+a)
    … …  
    5957
    6058function f(x)
    61 f=sin(x)*sin(2*x)*sin(3*x)
     59f=sin(x)*sin(2*x)*sin(3*x) !integral function
    6260end function f
    6361end
    6462}}}
    6563
    66 == Console3.f90 ==
     64== 3. Calculation dependence between integral function and coordinates ==
    6765{{{
    6866#!fortran
     67! Get a graph of the dependence in Origin, differentiate and integrate the graph and
     68!compare the resultants with the integrals which are calculated in the first program
    6969real x,dx, a, b
    70 open(1,file='y2.dat', status='unknown')
     70open(1,file='y2.dat', status='unknown')!the document y2 is in the project folder
    7171a=-0.5
    7272b=0.5
    … …  
    7474dx=(b-a)/n
    7575do i=1,n
    76     write(1,*) i*dx + a, 1/sqrt(1-(i*dx+a)**2)
     76    write(1,*) i*dx + a, 1/sqrt(1-(i*dx+a)**2)!integral function
    7777end do
    7878end
    7979}}}
    8080
    81 == Console4.f90 ==
     81== 4. Calculation dependence between integral function and coordinates ==
    8282{{{
    8383#!fortran
     84! Get a graph of the dependence in Origin, differentiate and integrate the graph and
     85!compare the resultant with the integral which is calculated in the second program
    8486real x,dx, a, b
    85 open(1,file='y2.dat', status='unknown')
     87open(1,file='y2.dat', status='unknown')!the document y2 is in the project folder
    8688a=0
    8789b=1.14
    … …  
    8991dx=(b-a)/n
    9092do i=1,n
    91     write(1,*) i*dx + a, sin(i*dx+a)*sin(2*i*dx+a)*sin(3*i*dx+a)
     93    write(1,*) i*dx + a, sin(i*dx+a)*sin(2*i*dx+a)*sin(3*i*dx+a)!integral function
    9294end domethod of
    9395end
    9496}}}
    9597
    96 == Console5.f90 ==
     98== 5. Simpson formula for Simpson method of integration ==
    9799{{{
    98100#!fortran
    99 program Simp
     101program Simpson
    100102
    101103integer n
    … …  
    106108   
    107109        x(i)=-1+h(i-1)
    108 
    109110        y1(i)=f1(x(i))
    110111        y2(i)=f2(x(i),h)
    111112        y3(i)=f3(x(i),h)
    112113       
    113          a=h/6*(y1+4*y2+y3)
     114         a=h/6*(y1+4*y2+y3)!Simpson's formula
    114115         s=sum(a)
    115116end do
    … …  
    118119real function f1(x)
    119120    real x
    120     f1=x*(5-4*x)**(-0.5)
     121    f1=x*(5-4*x)**(-0.5) !first function
    121122  end function
    122123real function f2(x,h)
    123124    real x, h
    124     f2=(x+h/2)*(5-4*(x+h/2))**(-0.5)
     125    f2=(x+h/2)*(5-4*(x+h/2))**(-0.5) !second function
    125126  end function
    126127real function f3(x,h)
    127128    real x, h
    128     f3=(x+h)*(5-4*(x+h))**(-0.5)
     129    f3=(x+h)*(5-4*(x+h))**(-0.5) !third function
    129130
    130131 end function
    131 print *, smethod of
     132print *, smethodof !Simpson method
    132133end program
    133134}}}
    134 
    135 == Console6.f90 ==
     135== 6. Characteristics of the electron motion, which is located between charged coaxial rings ==
    136136{{{
    137137#!fortran
    138 program kolco
    139 
    140 integer :: n=1000
    141 real :: r1=1.0, r2=2.0, q=3.0, d=5.0, eq=-4.8e-10, m=9.1e-28 ! ������� ��������� � ���-�������
    142 real ::  dt, t, x, vx
    143 
    144 t=1.0e-9
    145 dt=t/float(n)
    146 vx=0
    147 x=2.5
    148 
    149 
    150 open (1,file='kolco',status='unknown')
    151 
    152 do i=1,n
    153    E=q*x/(x**2+r1**2)**(1.5)- q*(5-x)/((5-x)**2+r2**2)**(1.5)
    154    vx=vx+(eq*E/m)*dt
    155    x=x+vx*dt
    156    write(1,*) vx, x, i*dt
    157 end do
    158 
    159 end program
    160 }}}
    161 
    162 == Console7.f90 ==
    163 {{{
    164 #!fortran
    165 program kolco
     138!Between two coaxial rings, charged with q1=q2=1nC and their radius are equal to r1=1cm
     139!and r2=2cm, are located at distance 5 cm between each other. At halfway between the rings is
     140!located electron with initial energy equal to 0. Find the dependencies between the
     141!velocity(energy)and the time,and electron coordinate with time. Also plot the phase
     142!trajectory V(x). And explain the resultants.   
     143program rings
    166144
    167145integer :: n=10000
    168 real(8) :: r1=1.0, r2=2.0, q=3.0, d=5.0, eq=4.8e-10, m=9.1e-28 ! ������� ��������� � ���-�������
     146real(8) :: r1=1.0, r2=2.0, q=3.0, d=5.0, eq=4.8e-10, m=9.1e-28 !the unis are in CGS system
    169147real(8) ::  dt, t, x, vx
    170148
    … …  
    175153
    176154
    177 open (1,file='kolco',status='unknown')
     155open (1,file='rings',status='unknown')
    178156
    179157do i=1,n
    180    E=q*x/sqrt((x**2+r1**2)**3)- q*(5-x)/sqrt(((5-x)**2+r2**2)**3)
    181    vx=vx+(eq*E/m)*dt
    182    x=x+vx*dt
     158   E=q*x/sqrt((x**2+r1**2)**3)- q*(5-x)/sqrt(((5-x)**2+r2**2)**3) !energy
     159   vx=vx+(eq*E/m)*dt !phase trajectory
     160   x=x+vx*dt !coordinate
    183161   write(1,*) x, i*dt, vx
    184162end do
    … …  
    187165}}}
    188166
    189 == Console8.f90 ==
     167== 7. Electron motion only in magnetic field ==
    190168{{{
    191169#!fortran
     170!Using the equation of electron motion find the phase trajectory. Electron filed E=0
     171!Part 1. Using only dependence of B0
     172integer :: n=1000
     173real(8) :: dt, t, vm, b, rl, om, x, y, vx, vy, B0
     174real(8) :: eq=-4.8e-10, me=9.1e-28, c = 3.0e10
     175!B0=me*c*om/ez
     176B0 = 200.0
     177!E0 = 5000
     178!A = 1e-30
     179x=1.0
     180y=1.0
     181vx=1000.0
     182vy=1000.0
     183a=100
     184!b = B0/B0
     185!om = eq*B0/me*c
     186!vm = A*om
     187!print *, vm
     188!rl = c/om
     189!vx = vm/vm
     190!vy = vm/vm
     191!x = x/rl
     192!y = y/rl
     193!print *, rl
     194t = 1.0e-8
     195dt = t/float(n)
     196
     197open (1,file='1',status='unknown')
     198
     199do i=1,n
     200     vx = vx+(B0*eq)/(c*me)*vy*dt
     201     vy = vy-(B0*eq)/(c*me)*vx*dt
     202     x = x+vx*dt
     203     y = y+vy*dt
     204   write(1,*) x, vx
     205end do
     206
     207pause
     208end program
     209}}}
     210
     211
     212== 8. Electron motion in electric and magnetic filed.  ==
     213{{{
     214#!fortran
     215!Like in program 7.Part 1
     216integer :: n=1000
     217real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b
     218real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0
     219x = 1.0
     220y = 1.0
     221a = 0.5
     222
     223vx = 1.0e09
     224vy = 1.0e09
     225E0=0.0
     226B0 = 200.0
     227
     228b=(2*3.14*6)/(1.0e09)
     229
     230!B0 = 1e-25
     231!E0 = 5000
     232!A = 1e-30
     233!b = B0/B0
     234!om = eq*B0/me*c
     235!vm = A*om
     236!print *, vm
     237!rl = c/om
     238!vx = vm/vm
     239!vy = vm/vm
     240!x = x/rl
     241!y = y/rl
     242!print *, rl
     243t = 1.00e-08
     244dt = t/float(n)
     245
     246
     247open (1,file='1',status='unknown')
     248
     249do i=1,n
     250   
     251   Bz = B0*(1+a*y)
     252   !Bz = B0
     253   x = x + vx*dt
     254   y = y + vy*dt
     255
     256   !vx = vx + (vy*eq*Bz)/(c*me)*dt
     257   !vy = vy - (vx*eq*Bz)/(c*me)*dt
     258
     259   
     260   if(i*dt.gt.b) then
     261E0 = 5.0
     262end if
     263   
     264  vx = vx +(eq*E0)/(me)*dt+(vy*eq*Bz)/(c*me)*dt
     265   vy = vy + (eq*E0)/(me)*dt-(vx*eq*Bz)/(c*me)*dt
     266  ! x = x + vx*dt
     267   !y = y + vy*dt
     268
     269   !print *, vx, vy
     270
     271
     272
     273 
     274
     275          write (1,*) x, y
     276
     277    end do
     278 pause 
     279end program
     280}}}
     281== 9. Electron motion in electric and magnetic filed ==
     282{{{
     283#!fortran
     284!Like in program 8. Part 2. the differnece is that here we are using different functions
     285!for vx and vy
     286integer :: n=1000
     287real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b
     288real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0
     289x = 1.0
     290y = 1.0
     291a = 0.5
     292
     293vx = 1.0e09
     294vy = 1.0e09
     295E0 = 0.0
     296B0 = 200.0
     297
     298b = (2*3.14)/(1.0e09)
     299
     300t = 1.00e-08
     301dt = t/float(n)
     302
     303
     304open (1,file='1',status='unknown')
     305
     306do i=1,n
     307   
     308   Bz = B0*(1+a*y)
     309   !Bz = B0
     310   x = x + vx*dt
     311   y = y + vy*dt
     312
     313   vx = vx + (vy*eq*Bz)/(c*me)*dt
     314   vy = vy - (vx*eq*Bz)/(c*me)*dt
     315
     316   
     317   if(i*dt.gt.b) then
     318E0 =  3.0
     319   end if
     320   
     321  vx = vx + (vy*eq*Bz)/(c*me)*dt
     322  vy = vy + (eq*E0)/(me)*dt  - (vx*eq*Bz)/(c*me)*dt
     323  x = x - vx*dt
     324  y = y + vy*dt
     325
     326   !print *, vx, vy
     327
     328
     329
     330 
     331
     332          write (1,*) x, y
     333
     334    end do
     335 pause 
     336end program
     337}}}
     338== 10. Electron motion only in magnetic field ==
     339{{{
     340#!fortran
     341! Like in program 7. Part 2 here we have dependence of magnetic filed in z coordinate
     342integer :: N = 1500
     343real(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
     344real(8) vx, vy, v, dt, t
     345
     346v = sqrt(2*We/m)
     347vx = sin(pi/4.0)*v
     348vy = cos(pi/4.0)*v
     349t = 5e-9
     350dt=t/float(N)
     351!Norma
     352!tau = t*10e8
     353!me = m*1e28
     354!qe = e*1e10
     355!vl = c*1e-10
     356!ve = v/c
     357!vx = sin(pi/4.0)*v
     358!vy = cos(pi/4.0)*v
     359!dt = tau/float(N)
     360
     361open (1,file='mov',status='replace')
     362do i=1,N
     363   Bz = B0*(1-b**2*(x**2+y**2))
     364   x = x + vx*dt
     365   y = y + vy*dt
     366   vx = vx + vy/m*Bz*e/c*dt
     367   vy = vy - vx/m*Bz*e/c*dt
     368   print *, x, y
     369   pause
     370   write(1,*) y, x
     371end do
     372end program
     373}}}
     374
     375 
     376
     377== 11. Determination the proton energy and power  ==
     378{{{
     379#!fortran
     380integer :: n=1000
     381real(8) :: dt, t, x, y, vx, vy, v, E, W
     382real(8) ::  mp = 1.67e-24, q = 4.8e-10
     383x = 1.0
     384y = 1.0
     385d = 1.0
     386a = 30
     387W0 = 1000
     388
     389v = sqrt(2*W0/mp)
     390
     391
     392t = 1.00e-05
     393dt = t/float(n)
     394
     395
     396open (1,file='1',status='unknown')
     397
     398do i=1,n
     399   
     400  vx = v*sin(a)
     401  vy = v*cos(a)
     402  x = x + vx*dt
     403   y = y + vy*dt
     404 
     405  E = (4*3.14*mp*y)/d
     406  W = (E*q)/ x
     407   
     408
     409   
     410
     411   !print *, vx, vy
     412
     413
     414
     415 
     416          write (1,*) W, t
     417
     418    end do
     419 pause 
     420end program
     421}}}
     422== 12. Runge-Kutta method for harmonic oscillations ==
     423{{{
     424#!fortran
     425!  metod Runge-Kuta. Oscillator. Part 1
     426
     427
    192428external fct, out
    193429real aux(8,2)
    194 common /a/ dt, k, n, A, Vm, om
    195 real :: pr(5) = (/0.,2*3.14159*10., 2*3.14159/50, .0001, .0/)
     430common /a/ b, lam, f0, dt, k, n
     431
     432real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)
     433real :: u(2) = (/0., 0.5/)
    196434real :: du(2) = (/0.5, 0.5/)
     435real :: b=0.003
    197436integer :: nd = 2
    198 real :: dt=2*3.14159265
    199 real :: t
    200 A = 2.0
    201 om = 1.0
    202 Vm = A*om
    203 u(1) = Vm/Vm
    204 u(2) = 0
    205 k = 1
    206 n=1
    207 
    208 open (unit=2, file = 'osc_rk0.dat' ,status='unknown')
    209 write(*,*) 'x=' , u(1), 'v=', u(2)
     437         
     438real :: dt = 2*3.14159265
     439   
     440! real :: lam = 0.1
     441real :: f0 = 0.00
     442 
     443open (unit=2, file = 'osc_rk.dat', status='new')
     444
     445write(*,*) 'x=', u(1), 'v=', u(2)
    210446pause
    211447call rkgs(pr, u, du, nd, ih, fct, out, aux)
     448
    212449write(2,*) 'ih=', ih
     450
    213451stop
    214452end
    215 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    216453
    217454subroutine fct (t, u, du)
    218 real u(2), du(2)
    219 common /a/ dt, k, n, A, Vm, om
    220 du(1) = -u(2)/(1+0.1*sin(om*dt))
     455
     456real u(2), du(2), lam
     457common /a/ b, lam, f0, dt, k, n
     458
     459du(1) = -u(2)-2.*b*u(1)                 
    221460du(2) = u(1)
    222 return
    223 end
    224 
    225 
    226 subroutine out (t, u, du, ih, nd, pr)
    227 real u(2), du(2), pr(5)
    228 common /a/ dt, k, n, A, Vm, om
    229 
    230 if(t.ge.k*(dt/50.)) then
    231 write (2,1) t/om, u(1)*Vm, u(2)*A
    232 k=k+1
    233 else
    234 end if
    235 
    236 1 format (f8.3, f8.3, f8.3)
    237 return
    238 end
    239 
    240 
    241 
    242 }}}
    243 
    244 == Console9.f90 ==
    245 {{{
    246 #!fortran
    247 !  Console9.f90
    248 !
    249 !  FUNCTIONS:
    250 !  Console9 - Entry point of console application.
    251 !
    252 
    253 !****************************************************************************
    254 !
    255 !  PROGRAM: Console9
    256 !
    257 !  PURPOSE:  Entry point for the console application.
    258 !
    259 !****************************************************************************
    260 
    261     program Console9
    262 
    263     implicit none
    264 
    265     ! Variables
    266 
    267     ! Body of Console9
    268     print *, 'Hello World'
    269 
    270     end program Console9
    271 
    272 }}}
    273 
    274 
    275 == Console10.f90 ==
    276 {{{
    277 #!fortran
    278 !  Osc_rk_0.f90
    279 !
    280 !  Lab #3
    281 !  metod Runge-Kuta. Oscillator.
    282 
    283 external fct, out
    284 real aux(8,4)
    285 common /a/  dt, k, n, vm, A, om, rl, x, y, r, u(1), u(2), u(3), u(4)
    286 
    287 real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)
    288 
    289 
    290 integer :: nd = 2
    291 real :: dt = 2*3.14159265
    292 real :: t
    293 real :: c = 3.0e10
    294 !B0=me*c*om/ez
    295 B0 = 200.0
    296 E0 = 5000
    297 eq = 4.8e-10
    298 me = 9.1e-28
    299 om = eq*B0/me*c
    300 A = 2.0
    301 vm = A*om
    302 n = 1
    303 k = 1
    304 b = B0/B0
    305 rl = c/om
    306 u(1) = vm/vm
    307 u(2) = vm/vm
    308 u(3) = x/rl
    309 u(4) = y/rl
    310 
    311 
    312 
    313  
    314 open (unit=2, file = '1_rk.dat', status='unknown')
    315 
    316 write(*,*) 'x=', u(3), 'vx=', u(1), 'y=', u(4), 'vy=', u(2)
    317 call rkgs(pr, u, du, nd, ih, fct, out, aux)
    318 
    319 write(2,*) 'ih=', ih
    320 
    321 
    322 end
    323 
    324 subroutine fct (t, u, du)
    325 
    326 real u(3), du(3), u(4), du(4)
    327 common /a/ dt, k, n, vm, A, om, rl, x, y, r,  u(1), u(2), u(3), u(4)
    328  
    329 du(1) = -u(3)
    330 du(2) = u(4)
    331 du(3) = u(1)
    332 du(4) = u(2)             
    333461                               
    334462return
    … …  
    336464
    337465subroutine  out (t, u, du, ih, nd, pr)
    338 
    339 common /a/ dt, k, n, om, vm, A, rl, x, y, r,  u(1), u(2), u(3), u(4)
     466real u(2), du(2), pr(5), lam
     467common /a/ b, lam, f0, dt, k, n
    340468
    341469if (t.ge.k*(dt/20.)) then
    342470
    343  write (4,1) t, u(1), u(2), u(3), u(4) 
     471 write (2,1) t,u(1),u(2)  !
    344472  k = k + 1
    345473   else
    … …  
    503631      END                                                               
    504632
    505 
    506633}}}
    507634
    508 == Console11.f90 ==
     635== 13. Runge-Kutta method for harmonic oscillations ==
    509636{{{
    510637#!fortran
    511 !  Console11.f90
    512 !
    513 !  FUNCTIONS:
    514 !  Console11 - Entry point of console application.
    515 !
    516 
    517 !****************************************************************************
    518 !
    519 !  PROGRAM: Console11
    520 !
    521 !  PURPOSE:  Entry point for the console application.
    522 !
    523 !****************************************************************************
    524 
    525     program Console11
    526 
    527     implicit none
    528 
    529     ! Variables
    530 
    531     ! Body of Console11
    532     print *, 'Hello World'
    533 
    534     end program Console11
    535 
    536 }}}
    537 
    538 == Console12.f90 ==
    539 {{{
    540 #!fortran
    541 !  Console12.f90
    542 !
    543 !  FUNCTIONS:
    544 !  Console12 - Entry point of console application.
    545 !
    546 
    547 !****************************************************************************
    548 !
    549 !  PROGRAM: Console12
    550 !
    551 !  PURPOSE:  Entry point for the console application.
    552 !
    553 !****************************************************************************
    554 
    555     program Console12
    556 
    557     implicit none
    558 
    559     ! Variables
    560 
    561     ! Body of Console12
    562     print *, 'Hello World'
    563 
    564     end program Console12
    565 
    566 }}}
    567 
    568 == Console13.f90 ==
    569 {{{
    570 #!fortran
    571 integer :: n=1000
    572 real(8) :: dt, t, vm, b, rl, om, x, y, vx, vy, B0
    573 real(8) :: eq=-4.8e-10, me=9.1e-28, c = 3.0e10
    574 !B0=me*c*om/ez
    575 B0 = 200.0
    576 !E0 = 5000
    577 !A = 1e-30
    578 x=1.0
    579 y=1.0
    580 vx=1000.0
    581 vy=1000.0
    582 a=100
    583 !b = B0/B0
    584 !om = eq*B0/me*c
    585 !vm = A*om
    586 !print *, vm
    587 !rl = c/om
    588 !vx = vm/vm
    589 !vy = vm/vm
    590 !x = x/rl
    591 !y = y/rl
    592 !print *, rl
    593 t = 1.0e-8
    594 dt = t/float(n)
    595 
    596 open (1,file='1',status='unknown')
    597 
    598 do i=1,n
    599      vx = vx+(B0*eq)/(c*me)*vy*dt
    600      vy = vy-(B0*eq)/(c*me)*vx*dt
    601      x = x+vx*dt
    602      y = y+vy*dt
    603    write(1,*) x, vx
    604 end do
    605 
    606 pause
    607 end program
    608 }}}
    609 
    610 == Console15.f90 ==
    611 {{{
    612 #!fortran
    613 x = (/1,5/)
    614 print *, x
    615 end
    616 }}}
    617 
    618 == Console16.f90 ==
    619 {{{
    620 #!fortran
    621 integer :: n=1000
    622 real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b
    623 real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0
    624 x = 1.0
    625 y = 1.0
    626 a = 0.5
    627 
    628 vx = 1.0e09
    629 vy = 1.0e09
    630 E0=0.0
    631 B0 = 200.0
    632 
    633 b=(2*3.14*6)/(1.0e09)
    634 
    635 !B0 = 1e-25
    636 !E0 = 5000
    637 !A = 1e-30
    638 !b = B0/B0
    639 !om = eq*B0/me*c
    640 !vm = A*om
    641 !print *, vm
    642 !rl = c/om
    643 !vx = vm/vm
    644 !vy = vm/vm
    645 !x = x/rl
    646 !y = y/rl
    647 !print *, rl
    648 t = 1.00e-08
    649 dt = t/float(n)
    650 
    651 
    652 open (1,file='1',status='unknown')
    653 
    654 do i=1,n
    655    
    656    Bz = B0*(1+a*y)
    657    !Bz = B0
    658    x = x + vx*dt
    659    y = y + vy*dt
    660 
    661    !vx = vx + (vy*eq*Bz)/(c*me)*dt
    662    !vy = vy - (vx*eq*Bz)/(c*me)*dt
    663 
    664    
    665    if(i*dt.gt.b) then
    666 E0 = 5.0
    667 end if
    668    
    669   vx = vx +(eq*E0)/(me)*dt+(vy*eq*Bz)/(c*me)*dt
    670    vy = vy + (eq*E0)/(me)*dt-(vx*eq*Bz)/(c*me)*dt
    671   ! x = x + vx*dt
    672    !y = y + vy*dt
    673 
    674    !print *, vx, vy
    675 
    676 
    677 
    678  
    679 
    680           write (1,*) x, y
    681 
    682     end do
    683  pause 
    684 end program
    685 }}}
    686 
    687 == Console17.f90 ==
    688 {{{
    689 #!fortran
    690 integer :: N = 1500
    691 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
    692 real(8) vx, vy, v, dt, t
    693 
    694 v = sqrt(2*We/m)
    695 vx = sin(pi/4.0)*v
    696 vy = cos(pi/4.0)*v
    697 t = 5e-9
    698 dt=t/float(N)
    699 !Norma
    700 !tau = t*10e8
    701 !me = m*1e28
    702 !qe = e*1e10
    703 !vl = c*1e-10
    704 !ve = v/c
    705 !vx = sin(pi/4.0)*v
    706 !vy = cos(pi/4.0)*v
    707 !dt = tau/float(N)
    708 
    709 open (1,file='mov',status='replace')
    710 do i=1,N
    711    Bz = B0*(1-b**2*(x**2+y**2))
    712    x = x + vx*dt
    713    y = y + vy*dt
    714    vx = vx + vy/m*Bz*e/c*dt
    715    vy = vy - vx/m*Bz*e/c*dt
    716    print *, x, y
    717    pause
    718    write(1,*) y, x
    719 end do
    720 end program
    721 }}}
    722 
    723 == Console18.f90 ==
    724 {{{
    725 #!fortran
    726 !  Osc_rk_0.f90
    727 !
    728 !  Lab #3
    729 !  metod Runge-Kuta. Oscillator.
     638
     639!  metod Runge-Kuta. Oscillator. Part 2 (the differnece between part 1 is that here we are
     640!  using onother function for du(1))
    730641
    731642external fct, out
    732643real aux(8,2)
    733 common /a/ b, lam, f0, dt, k, n
    734 
    735 real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)
     644common /a/ b, lam, f0, dt, k, n, c
     645
     646real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/)
    736647real :: u(2) = (/0., 0.5/)
    737648real :: du(2) = (/0.5, 0.5/)
    738649real :: b=0.003
    739650integer :: nd = 2
    740          
     651real::  c = 2.0
    741652real :: dt = 2*3.14159265
     653
    742654   
    743655! real :: lam = 0.1
    … …  
    761673
    762674real u(2), du(2), lam
    763 common /a/ b, lam, f0, dt, k, n
    764 
    765 du(1) = -u(2)-2.*b*u(1)                 
     675common /a/ b, lam, f0, dt, k, n, c
     676
     677du(1) = c*cos(0.5*t)-u(2)-2.*b*u(1)                 
    766678du(2) = u(1)
    767679                               
    … …  
    771683subroutine  out (t, u, du, ih, nd, pr)
    772684real u(2), du(2), pr(5), lam
    773 common /a/ b, lam, f0, dt, k, n
     685common /a/ b, lam, f0, dt, k, n, c
    774686
    775687if (t.ge.k*(dt/20.)) then
    … …  
    939851}}}
    940852
    941 == Console19.f90 ==
     853== 14. Runge-Kutta method for harmonic oscillations ==
    942854{{{
    943855#!fortran
    944 !  Osc_rk_0.f90
    945 !
    946 !  Lab #3
    947 !  metod Runge-Kuta. Oscillator.
    948 
     856!Part 3. Another example and using another function for du(1)
    949857external fct, out
    950858real aux(8,2)
    951 common /a/ b, lam, f0, dt, k, n, c
    952 
    953 real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/)
    954 real :: u(2) = (/0., 0.5/)
     859common /a/ dt, k, n, A, Vm, om
     860real :: pr(5) = (/0.,2*3.14159*10., 2*3.14159/50, .0001, .0/)
    955861real :: du(2) = (/0.5, 0.5/)
    956 real :: b=0.003
    957862integer :: nd = 2
    958 real::  c = 2.0
    959 real :: dt = 2*3.14159265
    960 
    961    
    962 ! real :: lam = 0.1
    963 real :: f0 = 0.00
    964 
    965 n = 1
     863real :: dt=2*3.14159265
     864real :: t
     865A = 2.0
     866om = 1.0
     867Vm = A*om
     868u(1) = Vm/Vm
     869u(2) = 0
    966870k = 1
    967  
    968 open (unit=2, file = 'osc_rk.dat', status='new')
    969 
    970 write(*,*) 'x=', u(1), 'v=', u(2)
     871n=1
     872
     873open (unit=2, file = 'osc_rk0.dat' ,status='unknown')
     874write(*,*) 'x=' , u(1), 'v=', u(2)
    971875pause
    972876call rkgs(pr, u, du, nd, ih, fct, out, aux)
    973 
    974877write(2,*) 'ih=', ih
    975 
    976878stop
    977879end
    978 
    979880subroutine fct (t, u, du)
    980 
    981 real u(2), du(2), lam
    982 common /a/ b, lam, f0, dt, k, n, c
    983 
    984 du(1) = c*cos(0.5*t)-u(2)-2.*b*u(1)                 
     881real u(2), du(2)
     882common /a/ dt, k, n, A, Vm, om
     883du(1) = -u(2)/(1+0.1*sin(om*dt))
    985884du(2) = u(1)
    986                                
    987 return
    988 end 
    989 
    990 subroutine  out (t, u, du, ih, nd, pr)
    991 real u(2), du(2), pr(5), lam
    992 common /a/ b, lam, f0, dt, k, n, c
    993 
    994 if (t.ge.k*(dt/20.)) then
    995 
    996  write (2,1) t,u(1),u(2)  !
    997   k = k + 1
    998    else
    999     end if
     885return
     886end
     887
     888
     889subroutine out (t, u, du, ih, nd, pr)
     890real u(2), du(2), pr(5)
     891common /a/ dt, k, n, A, Vm, om
     892
     893if(t.ge.k*(dt/50.)) then
     894write (2,1) t/om, u(1)*Vm, u(2)*A
     895k=k+1
     896else
     897end if
    1000898
    10018991 format (f8.3, f8.3, f8.3)
    1002 
    1003900return
    1004 end
    1005 
     901end
    1006902
    1007903      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
    … …  
    11571053
    11581054}}}
    1159 
    1160 
    1161 == Console20.f90 ==
     1055== 15. Electron acceleration in uniform electric field ==
    11621056{{{
    11631057#!fortran
    1164 !  Console20.f90
    1165 !
    1166 !  FUNCTIONS:
    1167 !  Console20 - Entry point of console application.
    1168 !
    1169 
    1170 !****************************************************************************
    1171 !
    1172 !  PROGRAM: Console20
    1173 !
    1174 !  PURPOSE:  Entry point for the console application.
    1175 !
    1176 !****************************************************************************
    1177 
    1178     program Console20
    1179 
    1180     implicit none
    1181 
    1182     ! Variables
    1183 
    1184     ! Body of Console20
    1185     print *, 'Hello World'
    1186 
    1187     end program Console20
    1188 
    1189 }}}
    1190 
    1191 == Console21.f90 ==
    1192 {{{
    1193 #!fortran
    1194 integer :: n=1000
    1195 real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b
    1196 real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0
    1197 x = 1.0
    1198 y = 1.0
    1199 a = 0.5
    1200 
    1201 vx = 1.0e09
    1202 vy = 1.0e09
    1203 E0 = 0.0
    1204 B0 = 200.0
    1205 
    1206 b = (2*3.14)/(1.0e09)
    1207 
    1208 t = 1.00e-08
    1209 dt = t/float(n)
    1210 
    1211 
    1212 open (1,file='1',status='unknown')
    1213 
    1214 do i=1,n
    1215    
    1216    Bz = B0*(1+a*y)
    1217    !Bz = B0
    1218    x = x + vx*dt
    1219    y = y + vy*dt
    1220 
    1221    vx = vx + (vy*eq*Bz)/(c*me)*dt
    1222    vy = vy - (vx*eq*Bz)/(c*me)*dt
    1223 
    1224    
    1225    if(i*dt.gt.b) then
    1226 E0 =  3.0
    1227    end if
    1228    
    1229   vx = vx + (vy*eq*Bz)/(c*me)*dt
    1230   vy = vy + (eq*E0)/(me)*dt  - (vx*eq*Bz)/(c*me)*dt
    1231   x = x - vx*dt
    1232   y = y + vy*dt
    1233 
    1234    !print *, vx, vy
    1235 
    1236 
    1237 
    1238  
    1239 
    1240           write (1,*) x, y
    1241 
    1242     end do
    1243  pause 
    1244 end program
    1245 }}}
    1246 
    1247 == Console22.f90 ==
    1248 {{{
    1249 #!fortran
    1250 integer :: n=1000
    1251 real(8) :: dt, t, x, y, vx, vy, v, E, W
    1252 real(8) ::  mp = 1.67e-24, q = 4.8e-10
    1253 x = 1.0
    1254 y = 1.0
    1255 d = 1.0
    1256 a = 30
    1257 W0 = 1000
    1258 
    1259 v = sqrt(2*W0/mp)
    1260 
    1261 
    1262 t = 1.00e-05
    1263 dt = t/float(n)
    1264 
    1265 
    1266 open (1,file='1',status='unknown')
    1267 
    1268 do i=1,n
    1269    
    1270   vx = v*sin(a)
    1271   vy = v*cos(a)
    1272   x = x + vx*dt
    1273    y = y + vy*dt
    1274  
    1275   E = (4*3.14*mp*y)/d
    1276   W = (E*q)/ x
    1277    
    1278 
    1279    
    1280 
    1281    !print *, vx, vy
    1282 
    1283 
    1284 
    1285  
    1286           write (1,*) W, t
    1287 
    1288     end do
    1289  pause 
    1290 end program
    1291 }}}
    1292 
    1293 == Console24.f90 ==
    1294 {{{
    1295 #!fortran
    1296 
    1297 !  Lab #1
    1298 !  metod Runge-Kuta. Rezonans i avtorezonans.
     1058!  Using metod Runge-Kuta. .
    12991059
    13001060external fct, out
    1301 real aux(8,2)
    1302 common /a/  dt, k, n, g0, om, B0, B, a
    1303 
    1304 real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/)
    1305 real :: u(2) = (/0., 0.5/)
    1306 real :: du(2) = (/0.5, 0.5/)
     1061real aux(8,4)
     1062common /a/  dt, k, n, vm, A, om, rl, x, y, r, u(1), u(2), u(3), u(4)
     1063
     1064real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)
     1065
    13071066
    13081067integer :: nd = 2
    13091068real :: dt = 2*3.14159265
    1310 
    1311 real :: t, m0, B, B0, a
    1312 real :: f0 = 0.00
    1313 f = 2.45e09
    1314 pi = 3.14159265
    1315 ez = 4.8e-10
    1316 m0 = 9.1e-28
    1317 c = 3.0e10
    1318 om = 2*pi*f
    1319 E = 2.
    1320 w = 1000.
    1321 a = 0.09
    1322 g0 = ez*E/(m0*c*om)
    1323 B0 = m0*c*om/ez
    1324 B = B0*(1+a*t)
    1325 b = B/B0
    1326 u(1) = pi
    1327 u(2) = w/511000.+1
    1328 write (*,*) g0
    1329 pause
    1330 
    1331 
     1069real :: t
     1070real :: c = 3.0e10
     1071!B0=me*c*om/ez
     1072B0 = 200.0
     1073E0 = 5000
     1074eq = 4.8e-10
     1075me = 9.1e-28
     1076om = eq*B0/me*c
     1077A = 2.0
     1078vm = A*om
     1079n = 1
     1080k = 1
     1081b = B0/B0
     1082rl = c/om
     1083u(1) = vm/vm
     1084u(2) = vm/vm
     1085u(3) = x/rl
     1086u(4) = y/rl
    13321087 
    1333 n = 1
    1334 k = 1
    1335  
    1336 open (unit=2, file = 'rez_avtorez.dat', status='unknown')
    1337 
    1338 
     1088open (unit=2, file = '1_rk.dat', status='unknown')
     1089
     1090write(*,*) 'x=', u(3), 'vx=', u(1), 'y=', u(4), 'vy=', u(2)
    13391091call rkgs(pr, u, du, nd, ih, fct, out, aux)
    13401092
    1341 write(2,*) 'ih=' ,ih
    1342 
    1343 
     1093write(2,*) 'ih=', ih
    13441094
    13451095end
    … …  
    13471097subroutine fct (t, u, du)
    13481098
    1349 real u(2), du(2), lam
    1350 common /a/ dt, k, n, g0, om, B, B0, a
    1351 
    1352 du(1) = (a*t+1-u(2))/u(2)+g0*sqrt(1./(u(2)**2-1))*sin(u(1))         
    1353 du(2) = -g0*sqrt(1.-1./u(2)**2)*cos(u(1))
    1354 
    1355 !write (*,*)  du(1), du(2)
    1356 !pause
     1099real u(3), du(3), u(4), du(4)
     1100common /a/ dt, k, n, vm, A, om, rl, x, y, r,  u(1), u(2), u(3), u(4)
     1101 
     1102du(1) = -u(3)
     1103du(2) = u(4)
     1104du(3) = u(1)
     1105du(4) = u(2)             
    13571106                               
    13581107return
    … …  
    13601109
    13611110subroutine  out (t, u, du, ih, nd, pr)
    1362 real u(2), du(2), pr(5), lam
    1363 common /a/ dt, k, n, g0, om, B, B0, a
    1364 
    1365 if (t.ge.k*(dt/2.)) then
    1366 
    1367  write (2,1) t,u(1),(u(2)-1)*511.
     1111
     1112common /a/ dt, k, n, om, vm, A, rl, x, y, r,  u(1), u(2), u(3), u(4)
     1113
     1114if (t.ge.k*(dt/20.)) then
     1115
     1116 write (4,1) t, u(1), u(2), u(3), u(4) 
    13681117  k = k + 1
    13691118   else
    13701119    end if
    13711120
    1372 1 format (3e10.3)
     11211 format (f8.3, f8.3, f8.3)
    13731122
    13741123return
    … …  
    15301279}}}
    15311280
    1532 == Console26.f90 ==
     1281
     1282
     1283
     1284
     1285== 16. Determination the conditions of electron capture in SGA mode and the bunch time ==
    15331286{{{
    15341287#!fortran
    1535 
    1536 !  Lab #1
    1537 !  metod Runge-Kuta. Rezonans i avtorezonans.
     1288!  Using metod Runge-Kuta. Resonant and auto resonant.
     1289! SGA (synchrotron gyro-magnetic auto resonant) is a self sustaining ECR plasma in
     1290!magnetic field which is rising in time.
    15381291
    15391292external fct, out
    1540 real aux(8,4)
    1541 common /a/  dt, k, n, g0, om, B0, B, a, y, rl
    1542 
    1543 real :: pr(5)=(/0., 2*3.14159*200., 2*3.14159/50, .0001, .0/)
    1544 real :: u(4) = (/0., 0., 0., 0./)
    1545 real :: du(4) = (/0.25, 0.25, 0.25, 0.25/)
    1546 
    1547 integer :: nd = 4
     1293real aux(8,2)
     1294common /a/  dt, k, n, g0, om, B0, B, a
     1295
     1296real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/)
     1297real :: u(2) = (/0., 0.5/)
     1298real :: du(2) = (/0.5, 0.5/)
     1299
     1300integer :: nd = 2
    15481301real :: dt = 2*3.14159265
    15491302
    … …  
    15561309c = 3.0e10
    15571310om = 2*pi*f
    1558 E0 = 5.
    1559 E = E0
    1560 a = 0.00034
     1311E = 2.
     1312w = 1000.
     1313a = 0.09
     1314g0 = ez*E/(m0*c*om)
    15611315B0 = m0*c*om/ez
    1562 g0 = E0/B0
    1563 B = B0
     1316B = B0*(1+a*t)
    15641317b = B/B0
    1565 rl = c/om
    1566 u(1) = 0.
    1567 u(2) = 0.
    1568 u(3) = 0.
    1569 u(4) = 0.
    1570 write (*,*) g0, B0, b, E, rl
     1318u(1) = pi
     1319u(2) = w/511000.+1
     1320write (*,*) g0
    15711321pause
    15721322
    … …  
    15761326k = 1
    15771327 
    1578 open (unit=2, file = 'rez_avtorez.dat', status='unknown')
     1328open (unit=2, file = 'res_autorer.dat', status='unknown')
    15791329
    15801330
    … …  
    15891339subroutine fct (t, u, du)
    15901340
    1591 real u(4), du(4), lam
    1592 common /a/ dt, k, n, g0, om, B, B0, a, y, rl
    1593 y=sqrt(1+u(1)**2+u(2)**2)
    1594 du(1) = -g0*cos(t)- u(2)*(1+a*t)/y       
    1595 du(2) = -g0*sin(t)+u(1)*(1+a*t)/y
    1596 du(3) = u(1)/y
    1597 du(4) = u(2)/y
    1598 
    1599 !write (*,*)  du(1), du(2), du(3), du(4)
     1341real u(2), du(2), lam
     1342common /a/ dt, k, n, g0, om, B, B0, a
     1343
     1344du(1) = (a*t+1-u(2))/u(2)+g0*sqrt(1./(u(2)**2-1))*sin(u(1)) !equation for electron phase
     1345                                                            !in SGA       
     1346du(2) = -g0*sqrt(1.-1./u(2)**2)*cos(u(1)) !equation for electron total energy in SGA
     1347
     1348!write (*,*)  du(1), du(2)
    16001349!pause
    16011350                               
    … …  
    16041353
    16051354subroutine  out (t, u, du, ih, nd, pr)
    1606 real u(4), du(4), pr(5), lam
    1607 common /a/ dt, k, n, g0, om, B, B0, a, y, rl
    1608 
    1609 if (t.ge.k*(dt)) then
    1610 
    1611  write (2,1) t,(y-1)*511. ,u(3)*rl,u(4)*rl
     1355real u(2), du(2), pr(5), lam
     1356common /a/ dt, k, n, g0, om, B, B0, a
     1357
     1358if (t.ge.k*(dt/2.)) then
     1359
     1360 write (2,1) t,u(1),(u(2)-1)*511.
    16121361  k = k + 1
    16131362   else
    16141363    end if
    16151364
    1616 1 format (5e12.3)
     13651 format (3e10.3)
    16171366
    16181367return
    16191368end
     1369!plot the dependencies of the phase and the time, between the energy and the time   
    16201370
    16211371
    … …  
    16231373!                                                                       
    16241374!                                                                       
    1625       DIMENSION Y(4),DERY(4),AUX(8,4),A(4),B(4),C(4),PRMT(5)           
     1375      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)           
    16261376      DO 1 I=1,NDIM                                                     
    16271377    1 AUX(8,I)=.06666667*DERY(I)                                       
    … …  
    17631513   36 IHLF=11                                                           
    17641514      CALL FCT(X,Y,DERY)                                               
    1765       GOTO 39     
     1515      GOTO 39                                                           
     1516   37 IHLF=12                                                           
     1517      GOTO 39                                                           
    17661518   38 IHLF=13                                                           
    17671519   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                               
    17681520   40 RETURN                                                           
    17691521      END                                                               
     1522
     1523
    17701524}}}
    1771  
    1772 
    1773 == Console29.f90 ==
     1525
     1526== 17. Determination the dependencies between the electron momentum and time, electron momentum and coordinate==
    17741527{{{
    17751528#!fortran
    1776 ! ! program El in a mirro trap,  parabolic approximation of the magnetic field
    1777 ! for a mirror trap. ECR
    1778 
    1779 real :: ez = 4.8e-10, c=3.0e10, me=9.1e-28, f=2.4e09, pi=3.14159265
    1780 
    1781 common/a/ gx(250),gy(250),x,y,z,ux,uy,uz,rl,gm,tp,pi,dt,cl2,al,dl,d
    1782 
    1783 
    1784 open (unit=1, file = 'input_t.dat')
    1785 open (unit=2, file = 'res.dat')
    1786 !
    1787 !
    1788   call cpu_time(start)
    1789 !
    1790 read(1,*) kt,kd,E,al,R,B00,x,y,z,W0,fi
    1791 write(2,1) kt,kd,E,al,R,B00,x,y,z,W0,fi
    1792 1 format('kt='I7,2x,'kd=',I3,2x,'E(kV/cm)=',f4.1,2x,'al=',e10.3,2x,    &
    1793 'R=',f5.3,2x,'B00(%)=',f8.5/'x(cm) =',f4.1,2x,'y(cm) =',f4.1,2x,    &
    1794 'z(cm)=',f4.1,2x,'W0(eV)=',e10.3,2x,'fi=',f5.1)
    1795 2 format('g0=',e10.3,2x,'B0=',f6.1,2x,'cl=',e10.3,2x,'rl =',f6.2/    &
    1796 'wmax=',e10.3,2x,'v0=',e10.3,2x,'w00=',e10.3,2x,'rc=',e10.3)
    1797 
    1798 pi2=2.0*pi
    1799 
    1800 DL=8.0
    1801 D=12.0
    1802 fi=fi/180.*pi
    1803 
    1804 v0=sqrt(1.-1./(1.+w0/511000.)**2)  ! ��������� �������� ��������� � �������� �
    1805 
    1806 write(*,*) v0*3.0e8                ! ��������� �������� ��������� � m/c
    1807 U0=v0/sqrt(1.0-v0**2)              ! ��������� ������� ��������� � �������� mc
    1808 W00=511000.*(sqrt(1.0+u0**2)-1.0)  ! ���� ��� ��������� �������
    1809 om=2.0*pi*f
    1810 rl=c/om                            ! �������������� ������ ������������� ��������
    1811 zm=dl/2.0                          ! ������� ����������
    1812 
    1813                                    ! ������ ��������� cl, ������������� �������� ���������� ����
    1814                                    ! �������, ����� �������� �������� ���������� ���������
    1815  cl=(zm/rl)/sqrt(R-1.0)
    1816  cl2=cl*cl                         
    1817  zm=zm/rl                           ! ����������
    1818 
    1819 E=E*1.0e05/3.0e04                  ! ��������� ��� ���� � ����
    1820 g0=ez*E/(me*c*om)                  ! ������������ ��������� ��� ����
    1821 
    1822 B0=me*c*om/ez                      ! ����������� �������� �������� ���������� ���� - ECR ��� ��������� � ������ �����
    1823 wmax=511.0*2.0*g0**(2.0/3.0)        ! ������������ �������� ������� ��� ��� � ���������� ������������� � ��������� ����
    1824 rc=v0*3.0e10/om                    !
     1529
     1530!  Using the same condition like in program 16
     1531!  metod Runge-Kuta.
     1532external fct, out
     1533real aux(8,4)
     1534common /a/  dt, k, n, g0, om, B0, B, a, y, rl
     1535
     1536real :: pr(5)=(/0., 2*3.14159*200., 2*3.14159/50, .0001, .0/)
     1537real :: u(4) = (/0., 0., 0., 0./)
     1538real :: du(4) = (/0.25, 0.25, 0.25, 0.25/)
     1539
     1540integer :: nd = 4
     1541real :: dt = 2*3.14159265
     1542
     1543real :: t, m0, B, B0, a
     1544real :: f0 = 0.00
     1545f = 2.45e09
     1546pi = 3.14159265
     1547ez = 4.8e-10
     1548m0 = 9.1e-28
     1549c = 3.0e10
     1550om = 2*pi*f
     1551E0 = 5.
     1552E = E0
     1553a = 0.00034
     1554B0 = m0*c*om/ez
     1555g0 = E0/B0
     1556B = B0
     1557b = B/B0
     1558rl = c/om
     1559u(1) = 0.
     1560u(2) = 0.
     1561u(3) = 0.
     1562u(4) = 0.
     1563write (*,*) g0, B0, b, E, rl
     1564pause
     1565
     1566
    18251567 
    1826 !    B0=875                        ! ��������� ���� � ������� � �������������� ������ �������
    1827 !    B=875
    1828  write(2,2) g0, B0, cl, rl, wmax, v0, w00, rc
    1829      B=B/B0                     ! ���������������� �������� ���������� ����
    1830         x=x/rl           
    1831         y=y/rl
    1832         z=z/rl
    1833         ux=u0*sin(fi)               ! ����������� ������� �������� � ��������� ������ �������
    1834         uy=u0*cos(fi)
    1835         uz=0.                        ! � ����������� Z ��������� ������� ����� ����.
    1836         km=0. 
    1837         kp=0
    1838 
    1839 !    ������ ���� ��������������
    1840         m=250
    1841         dt=2.*pi/m
    1842          
    1843 !    ������ ���������������� ����
    1844         do i=1,m
    1845         gx(i)=-g0*cos(pi2*(i-1)/m)*dt/2.0       
    1846         gy(i)=-g0*sin(pi2*(i-1)/m)*dt/2.0       
    1847         end do
    1848 
    1849         write(2,*) 'wmax=', 2.*g0**(2./3.)*511.
    1850         wm=0
    1851 
    1852 !  ��������� ����
    1853         do k=1,kt
    1854         l=k-kp
    1855         tm=k*dt
     1568n = 1
     1569k = 1
    18561570 
    1857         tp=-(1.+b00)*dt/2.         
    1858         call motion(l)
    1859 if (k.eq.km+kd) then   ! ����� ����������� ����� kd ��������� �����
    1860         w=(gm-1.0)*511.0
    1861     write(2,3) tm/(2.*pi), x*rl,y*rl,z*rl, w
    1862 
    1863         km=km+kd
    1864         else
    1865         end if
    1866 if (k.eq.kp+m) then  ! ������� ��������� ����� ��� ��� ����
    1867         kp=kp+m
    1868         else
    1869         end if
    1870         end do
    1871 3 format(7f12.5)
    1872 
    1873 
    1874      call cpu_time(finish)
    1875         write(2,*) 'cpu-time =', finish-start
    1876         end
    1877                                              
    1878 subroutine motion(l)   
    1879 common/a/ gx(250),gy(250),x,y,z,ux,uy,uz,rl,gm,tp,pi,dt,cl2,al,dl,d
    1880 
    1881 ! �������� ��������� ���� ������� ���������� ���� (���� 2)
    1882 
    1883       r = sqrt(x*x + y*y)
    1884       bxp=-x*z/cl2
    1885       byp=-y*z/cl2
    1886       bzp=1.0+z*z/cl2-r*r/(2.*cl2)
    1887 
    1888 !      bxp=0.0
    1889 !      byp=0.0
    1890 !      bzp=1.0
    1891 ! ����� ������
    1892       u=sqrt(ux**2+uy**2)
    1893 
    1894 ! �������� ��� ������������� ���� (���� 3)
    1895       uxm = ux  + gx(l)
    1896       uym = uy  + gy(l)   
    1897       uzm = uz
    1898 
    1899       gmn = sqrt(1 + uxm*uxm + uym*uym + uzm*uzm)
    1900 
    1901       tg = tp/gmn
    1902       tx = tg * bxp
    1903       ty = tg * byp
    1904       tz = tg * bzp
    1905 
    1906       uxr = uxm + uym*tz - uzm*ty
    1907       uyr = uym + uzm*tx - uxm*tz
    1908       uzr = uzm + uxm*ty - uym*tx
    1909       gmn = sqrt(1 + uxm*uxm + uym*uym + uzm*uzm)
    1910 
    1911       txyzr= 1.+ tx*tx + ty*ty + tz*tz
    1912 
    1913       sx = 2*tx/txyzr
    1914       sy = 2*ty/txyzr
    1915       sz = 2*tz/txyzr
    1916 
    1917       uxp = uxm + uyr*sz - uzr*sy
    1918       uyp = uym + uzr*sx - uxr*sz
    1919       uzp = uzm + uxr*sy - uyr*sx
    1920 
    1921       ux = uxp + gx(l)
    1922       uy = uyp + gy(l)     
    1923       uz = uzp
    1924 
    1925       gm = sqrt (1. + ux**2 + uy**2 + uz**2)
    1926       gt=dt/gm
    1927       x = x + ux*gt
    1928       y = y + uy*gt
    1929       z = z + uz*gt
    1930 
    1931       if (abs(z*rl).ge.dl/2.0.or.r*rl.ge.d/2.0) then ! ������� ���������
    1932                                                       ! ��������� �� ������ ������
    1933       write(*,*) 'out',x*rl,y*rl,z*rl
    1934       stop
    1935       else
    1936       endif
    1937 
    1938       return
    1939       end
    1940 
    1941 
    1942 }}}
    1943 
    1944 
    1945 == Console30.f90 ==
    1946 {{{
    1947 #!fortran
    1948 !  Osc_rk_0.f90
    1949 !
    1950 !  Lab #3
    1951 !  metod Runge-Kuta. Oscillator.
    1952 
    1953 external fct, out
    1954 real aux(8,2)
    1955 common /a/  dt, k, n, vm, A, om
    1956 
    1957 real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)
    1958 real :: u(2) = (/0., 0.5/)
    1959 real :: du(2) = (/0.5, 0.5/)
    1960 integer :: nd = 2
    1961 real :: dt = 2*3.14159265
    1962 ! real :: lam = 0.1
    1963 real :: t
    1964 real :: om=1
    1965 A=2.0
    1966 vm=A*om
    1967 n = 1
    1968 k = 1
    1969 u(1)=vm/vm
    1970 u(2)=0
    1971  
    1972 open (unit=2, file = 'osc_rk.dat', status='unknown')
     1571open (unit=2, file = 'res_autores.dat', status='unknown')
    19731572
    19741573
    19751574call rkgs(pr, u, du, nd, ih, fct, out, aux)
    19761575
    1977 write(2,*) 'ih=', ih
     1576write(2,*) 'ih=' ,ih
     1577
    19781578
    19791579
    … …  
    19821582subroutine fct (t, u, du)
    19831583
    1984 real u(2), du(2), lam
    1985 common /a/ dt, k, n, g, vm, A, om
    1986 
    1987 du(1) = -u(2)/(1+1.2*sin(t*2))               
    1988 du(2) = u(1)
     1584real u(4), du(4), lam
     1585common /a/ dt, k, n, g0, om, B, B0, a, y, rl
     1586y=sqrt(1+u(1)**2+u(2)**2)
     1587du(1) = -g0*cos(t)- u(2)*(1+a*t)/y       
     1588du(2) = -g0*sin(t)+u(1)*(1+a*t)/y
     1589du(3) = u(1)/y
     1590du(4) = u(2)/y
     1591
     1592!write (*,*)  du(1), du(2), du(3), du(4)
     1593!pause
    19891594                               
    19901595return
    … …  
    19921597
    19931598subroutine  out (t, u, du, ih, nd, pr)
    1994 real u(2), du(2), pr(5), lam
    1995 common /a/ dt, k, n, om, vm, A
    1996 
    1997 if (t.ge.k*(dt/20.)) then
    1998 
    1999  write (2,1) t,u(1),u(2)  !
     1599real u(4), du(4), pr(5), lam
     1600common /a/ dt, k, n, g0, om, B, B0, a, y, rl
     1601
     1602if (t.ge.k*(dt)) then
     1603
     1604 write (2,1) t,(y-1)*511. ,u(3)*rl,u(4)*rl
    20001605  k = k + 1
    20011606   else
    20021607    end if
    20031608
    2004 1 format (f8.3, f8.3, f8.3)
     16091 format (5e12.3)
    20051610
    20061611return
    … …  
    20111616!                                                                       
    20121617!                                                                       
    2013       DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)           
     1618      DIMENSION Y(4),DERY(4),AUX(8,4),A(4),B(4),C(4),PRMT(5)           
    20141619      DO 1 I=1,NDIM                                                     
    20151620    1 AUX(8,I)=.06666667*DERY(I)                                       
    … …  
    21511756   36 IHLF=11                                                           
    21521757      CALL FCT(X,Y,DERY)                                               
    2153       GOTO 39                                                           
    2154    37 IHLF=12                                                           
    2155       GOTO 39                                                           
     1758      GOTO 39     
    21561759   38 IHLF=13                                                           
    21571760   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                               
    21581761   40 RETURN                                                           
    21591762      END                                                               
    2160 
    2161 
    21621763}}}
    21631764
    2164 == Console31.f90 ==
     1765 
     1766
     1767== 18. Electron acceleration in a uniform electric field using speed light ==
    21651768{{{
    21661769#!fortran
    2167     program heat1
    2168 
    2169     implicit none
    2170 
    2171     ! Variables
    2172     real a,b,c,d,dx,dt,x_i,t_j
    2173     integer n,m,i,j
    2174    
    2175     real, allocatable:: U(:,:)
    2176    
    2177     open(10,file='result_out1.txt')
    2178     a=1
    2179     b=4
    2180     c=1
    2181     d=15
    2182     write(0,*) "vvedite melkost razbijenija po x"
    2183     read(*,*) n
    2184    
    2185    ! do i=1,n
    2186      
    2187      write(0,*) "vvedite melkost razbijenija po t"
    2188         read(*,*) m
    2189         dx=(b-a)/n
    2190          dt=(d-c)/m
    2191         allocate(U(n,m))
    2192        
    2193        
    2194        do j=1,m      ! ��������� �������
    2195         t_j=j*dt
    2196         U(1,j)=3
    2197         U(n,j)=1
    2198         end do
    2199        
    2200         do i=2,n-1   ! ��������� ������� (�������  ����� �� ������, �.�. ���� ������ ���������)
    2201          x_i=i*dx
    2202         U(i,1)=1
    2203         end do
    2204        
    2205        
    2206         ! �������� ����
    2207          do j=1,m-1        ! ���� �� �������
    2208          do i=2,n-1        ! ���� �� ������������
    2209          U(i,j+1)=U(i,j)+0.005*(dt/(dx)**2)*(U(i+1,j)-2*U(i,j)+U(i-1,j))
    2210          print *,i,j+1,U(i,j+1)
    2211          end do
    2212          end do
    2213    
    2214    
    2215         !print *, x_i
    2216     !end do
    2217    
    2218         !do j=1,m
    2219        
    2220    
    2221    
    2222     !print *, t_j
    2223    
    2224    
    2225     do i=1,n
    2226     write(10,*) i*dx,U(i,m/100),U(i,m/2),U(i,m)
    2227     enddo
    2228    
    2229 
    2230     end program heat1
    2231 
    2232 }}}
    2233 
    2234 == Console32.f90 ==
    2235 {{{
    2236 #!fortran
    2237 !  Osc_rk_0.f90
    2238 !
    2239 !  Lab #3
    2240 !  metod Runge-Kuta. Oscillator.
     1770! Using metod Runge-Kuta
    22411771
    22421772external fct, out
    … …  
    24671997}}}
    24681998
    2469 == Console33.f90 - Poisson solver ==
     1999
     2000== 19. Charged particle motion in magnetic mirror trap under ECR ==
    24702001{{{
    24712002#!fortran
    2472 
     2003! Using Boris method get the equations of the charge particle motion. And analyze the ECR
     2004!(electron cyclotron resonance) phenomenon in parabolic approximation of the magnetic
     2005!field (mirror trap) for different input values
     2006
     2007real :: ez = 4.8e-10, c=3.0e10, me=9.1e-28, f=2.4e09, pi=3.14159265
     2008
     2009common/a/ gx(250),gy(250),x,y,z,ux,uy,uz,rl,gm,tp,pi,dt,cl2,al,dl,d
     2010
     2011
     2012open (unit=1, file = 'input_t.dat')
     2013open (unit=2, file = 'res.dat')
     2014!
     2015!
     2016  call cpu_time(start)
     2017!
     2018read(1,*) kt,kd,E,al,R,B00,x,y,z,W0,fi
     2019write(2,1) kt,kd,E,al,R,B00,x,y,z,W0,fi
     20201 format('kt='I7,2x,'kd=',I3,2x,'E(kV/cm)=',f4.1,2x,'al=',e10.3,2x,    &
     2021'R=',f5.3,2x,'B00(%)=',f8.5/'x(cm) =',f4.1,2x,'y(cm) =',f4.1,2x,    &
     2022'z(cm)=',f4.1,2x,'W0(eV)=',e10.3,2x,'fi=',f5.1)
     20232 format('g0=',e10.3,2x,'B0=',f6.1,2x,'cl=',e10.3,2x,'rl =',f6.2/    &
     2024'wmax=',e10.3,2x,'v0=',e10.3,2x,'w00=',e10.3,2x,'rc=',e10.3)
     2025
     2026pi2=2.0*pi
     2027
     2028DL=8.0
     2029D=12.0
     2030fi=fi/180.*pi
     2031
     2032v0=sqrt(1.-1./(1.+w0/511000.)**2) 
     2033
     2034write(*,*) v0*3.0e8               
     2035U0=v0/sqrt(1.0-v0**2)             
     2036W00=511000.*(sqrt(1.0+u0**2)-1.0)
     2037om=2.0*pi*f
     2038rl=c/om                           
     2039zm=dl/2.0                         
     2040 cl=(zm/rl)/sqrt(R-1.0)
     2041 cl2=cl*cl                         
     2042 zm=zm/rl                         
     2043
     2044E=E*1.0e05/3.0e04                 
     2045g0=ez*E/(me*c*om)                 
     2046B0=me*c*om/ez                     
     2047wmax=511.0*2.0*g0**(2.0/3.0)       
     2048rc=v0*3.0e10/om                   
     2049 
     2050!    B0=875                       
     2051!    B=875
     2052 write(2,2) g0, B0, cl, rl, wmax, v0, w00, rc
     2053     B=B/B0                     
     2054        x=x/rl           
     2055        y=y/rl
     2056        z=z/rl
     2057        ux=u0*sin(fi)               
     2058        uy=u0*cos(fi)
     2059        uz=0.                       
     2060        km=0. 
     2061        kp=0
     2062
     2063        m=250
     2064        dt=2.*pi/m
     2065         
     2066
     2067        do i=1,m
     2068        gx(i)=-g0*cos(pi2*(i-1)/m)*dt/2.0       
     2069        gy(i)=-g0*sin(pi2*(i-1)/m)*dt/2.0       
     2070        end do
     2071
     2072        write(2,*) 'wmax=', 2.*g0**(2./3.)*511.
     2073        wm=0
     2074
     2075
     2076        do k=1,kt
     2077        l=k-kp
     2078        tm=k*dt
     2079 
     2080        tp=-(1.+b00)*dt/2.         
     2081        call motion(l)
     2082if (k.eq.km+kd) then 
     2083        w=(gm-1.0)*511.0
     2084    write(2,3) tm/(2.*pi), x*rl,y*rl,z*rl, w
     2085
     2086        km=km+kd
     2087        else
     2088        end if
     2089if (k.eq.kp+m) then 
     2090        kp=kp+m
     2091        else
     2092        end if
     2093        end do
     20943 format(7f12.5)
     2095
     2096
     2097     call cpu_time(finish)
     2098        write(2,*) 'cpu-time =', finish-start
     2099        end
     2100                                             
     2101subroutine motion(l)   
     2102common/a/ gx(250),gy(250),x,y,z,ux,uy,uz,rl,gm,tp,pi,dt,cl2,al,dl,d
     2103
     2104
     2105
     2106      r = sqrt(x*x + y*y)
     2107      bxp=-x*z/cl2
     2108      byp=-y*z/cl2
     2109      bzp=1.0+z*z/cl2-r*r/(2.*cl2)
     2110
     2111!      bxp=0.0
     2112!      byp=0.0
     2113!      bzp=1.0
     2114
     2115      u=sqrt(ux**2+uy**2)
     2116
     2117!Boris's scheme of partical motion
     2118      uxm = ux  + gx(l)
     2119      uym = uy  + gy(l)   
     2120      uzm = uz
     2121
     2122      gmn = sqrt(1 + uxm*uxm + uym*uym + uzm*uzm)
     2123
     2124      tg = tp/gmn
     2125      tx = tg * bxp
     2126      ty = tg * byp
     2127      tz = tg * bzp
     2128
     2129      uxr = uxm + uym*tz - uzm*ty
     2130      uyr = uym + uzm*tx - uxm*tz
     2131      uzr = uzm + uxm*ty - uym*tx
     2132      gmn = sqrt(1 + uxm*uxm + uym*uym + uzm*uzm)
     2133
     2134      txyzr= 1.+ tx*tx + ty*ty + tz*tz
     2135
     2136      sx = 2*tx/txyzr
     2137      sy = 2*ty/txyzr
     2138      sz = 2*tz/txyzr
     2139
     2140      uxp = uxm + uyr*sz - uzr*sy
     2141      uyp = uym + uzr*sx - uxr*sz
     2142      uzp = uzm + uxr*sy - uyr*sx
     2143
     2144      ux = uxp + gx(l)
     2145      uy = uyp + gy(l)     
     2146      uz = uzp
     2147
     2148      gm = sqrt (1. + ux**2 + uy**2 + uz**2)
     2149      gt=dt/gm
     2150      x = x + ux*gt
     2151      y = y + uy*gt
     2152      z = z + uz*gt
     2153
     2154      if (abs(z*rl).ge.dl/2.0.or.r*rl.ge.d/2.0) then
     2155                                                     
     2156      write(*,*) 'out',x*rl,y*rl,z*rl
     2157      stop
     2158      else
     2159      endif
     2160
     2161      return
     2162      end
     2163
     2164
     2165}}}
     2166
     2167
     2168== 20. Solving nonlinear equation ==
     2169{{{
     2170#!fortran
     2171!using method Newton-Raphson
     2172    program heat1
     2173
     2174    implicit none
     2175
     2176    ! Variables
     2177    real a,b,c,d,dx,dt,x_i,t_j
     2178    integer n,m,i,j
     2179   
     2180    real, allocatable:: U(:,:)
     2181   
     2182    open(10,file='result_out1.txt')
     2183    a=1
     2184    b=4
     2185    c=1
     2186    d=15
     2187    write(0,*) "vvedite melkost razbijenija po x"
     2188    read(*,*) n
     2189   
     2190   ! do i=1,n
     2191     
     2192     write(0,*) "vvedite melkost razbijenija po t"
     2193        read(*,*) m
     2194        dx=(b-a)/n
     2195         dt=(d-c)/m
     2196        allocate(U(n,m))
     2197       
     2198       
     2199       do j=1,m     
     2200        t_j=j*dt
     2201        U(1,j)=3
     2202        U(n,j)=1
     2203        end do
     2204       
     2205        do i=2,n-1   
     2206         x_i=i*dx
     2207        U(i,1)=1
     2208        end do
     2209       
     2210       
     2211       
     2212         do j=1,m-1       
     2213         do i=2,n-1       
     2214         U(i,j+1)=U(i,j)+0.005*(dt/(dx)**2)*(U(i+1,j)-2*U(i,j)+U(i-1,j))
     2215         print *,i,j+1,U(i,j+1)
     2216         end do
     2217         end do
     2218   
     2219   
     2220        !print *, x_i
     2221    !end do
     2222   
     2223        !do j=1,m
     2224       
     2225   
     2226   
     2227    !print *, t_j
     2228   
     2229   
     2230    do i=1,n
     2231    write(10,*) i*dx,U(i,m/100),U(i,m/2),U(i,m)
     2232    enddo
     2233   
     2234
     2235    end program heat1
     2236
     2237}}}
     2238
     2239== 21. Poisson equations solver ==
     2240{{{
     2241#!fortran
     2242! Using swap method in one dimensional case, determine the electric filed of the
     2243! electron layer
     2244! Get the plots of the dependencies between the density and the electric filed with   
     2245! coordinate z
    24732246      program poisn1
    24742247real d,n0
    … …  
    24982271      pi=3.14159265
    24992272      n0=10e+08
    2500       d=0.4            ! ⮫騭� ᫮� � ������      thick of a layer
    2501       jz=400             ! �������⢮ ���������� �祥�  number of filled cells
    2502       del=d/jz            ! 蠣 ��⪨ � ������   step (in meters)
     2273      d=0.4               !  thick of a layer
     2274      jz=400              !  number of filled cells
     2275      del=d/jz            ! step (in meters)
    25032276      !z=2
    25042277     
    2505      
    2506           ! charge dencity
     2278! charge dencity
    25072279     
    25082280      rho=n0*abs(sin(2*pi*k/400))           
    25092281       qz=rho
    25102282!
    2511 !  ��砫쭮� ����।������  spatial initial distribution
     2283!  spatial initial distribution
    25122284!
    25132285      do k=1,j
    … …  
    25172289      end do
    25182290!     
    2519 !      �맮� ����ணࠬ�� ��襭�� �ࠢ����� ����ᮭ�   Poisson equation solver
     2291!    Poisson equation solver
    25202292!     
    25212293      call pnsolv
    … …  
    25252297      do k=2,j-1
    25262298      Ez(k)=(fi(k-1)-fi(k+1))/(2.*del)
    2527       write(7,7) k,q(k),fi(k),Ez(k)     ! � �/�     (V/m)
     2299      write(7,7) k,q(k),fi(k),Ez(k)     !   (V/m)
    25282300      end do
    25292301!
    … …  
    25312303!
    25322304!      Ean=               
    2533       write(7,8) Ean                    ! � �/�     (V/m)
     2305      write(7,8) Ean                    !     (V/m)
    25342306    7 format(3x,I3,3x,3e15.5) 
    25352307    8 format(3x,e15.5) 
    … …  
    25452317      dimension z(1000),y(1000)
    25462318!
    2547 !            �ண���� ��ࠢ� ������ (calculation  z and y coefficients)
     2319!     (calculation  z and y coefficients)
    25482320! 
    25492321      j1=j+1
    … …  
    25592331      end do
    25602332!
    2561 !            �ண���� ᫥�� ���ࠢ� (calculation  fi(j)  at grid points)
     2333!      (calculation  fi(j)  at grid points)
    25622334!
    25632335      f0=(al*y(1)-cl)/(al-bl-al*z(1))
    … …  
    25702342}}}
    25712343
    2572 == Console35.f90 - 1D plasma simulation ==
     2344== 22. Distribution of Langmuir waves in a uniform plasma ==
    25732345{{{
    25742346#!fortran
    25752347!
    25762348! 1D plasma simulation (demo version)
    2577 !
     2349! using numerical method created by particle cell method
     2350! calculate the wavelength and the velocity of wave propogation
     2351
    25782352 program plasma_1D
    25792353   real n,m0
    … …  
    26052379      dt=dt*2.0*pi
    26062380
    2607       d=0.1e-1            ! ⮫騭� ᫮� � ������      thick of a layer
    2608       jz=4000              ! �������⢮ ���������� �祥�  number of filled cells
    2609       del=d/jz            ! 蠣 ��⪨ � ������   step (in meters)
     2381      d=0.1e-1            !   thick of a layer
     2382      jz=4000             !   number of filled cells
     2383      del=d/jz            !   step (in meters)
    26102384
    26112385      rho=-1.0e16         ! charge dencity
    … …  
    26522426      kp=50
    26532427!
    2654 !           Cycle in time
     2428!     Cycle in time
    26552429!
    26562430      do k=1,kt
    … …  
    27482522
    27492523
    2750 == Console36.f90 - Proton acceleration ==
     2524== 23. Proton acceleration by the filed of electron leyer ==
    27512525{{{
    27522526#!fortran
    27532527
    27542528
    2755 !method Runge-Kutta
     2529!Using method Runge-Kutta find the max possible acceleration of the proton (maximum
     2530!gradient of the magnetic filed). Plot the graphics dependencies between the power and
     2531!time, coordinate difference and time 
    27562532
    27572533
    … …  
    29452721   17 AUX(4,I)=Y(I)                                                     
    29462722      ITEST=1                                                           
    2947       ISTEP=ISTEP+ISTEP-2                          Console2.f90                    
     2723      ISTEP=ISTEP+ISTEP-2                                             
    29482724   18 IHLF=IHLF+1
    29492725      X=X-H                                                             
    … …  
    30212797}}}
    30222798
    3023 == Console37.f90 - Euler ==
     2799== 24. Electron motion in uniform electric filed ==
    30242800{{{
    30252801#!fortran
     2802!using method Euler
    30262803real v,x,dt,t,dv,Fm
    30272804open(2,file='y2.dat', status='unknown')