[[PageOutline]] = FORTRAN 95 examples = == Console1.f90 == {{{ #!fortran 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 S2=S2+ f(i*dx+a)*dx if(mod(i,2)==0) then S3=S3+dx/3*2*f(i*dx+a) 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) end function f end }}} == Console2.f90 == {{{ #!fortran 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 == {{{ #!fortran 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 == {{{ #!fortran 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 do end }}} == Console5.f90 == {{{ #!fortran 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 *, s end program }}} == Console6.f90 == {{{ #!fortran 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 == {{{ #!fortran 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 == {{{ #!fortran 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 == {{{ #!fortran ! 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 == {{{ #!fortran ! 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 == {{{ #!fortran ! 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 == {{{ #!fortran ! 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 == {{{ #!fortran 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 == {{{ #!fortran x = (/1,5/) print *, x end }}} == Console16.f90 == {{{ #!fortran 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 == {{{ #!fortran 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 == {{{ #!fortran ! 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 == {{{ #!fortran ! 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 == {{{ #!fortran ! 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 == {{{ #!fortran 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 == {{{ #!fortran 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 == {{{ #!fortran ! 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 == {{{ #!fortran ! 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 == {{{ #!fortran ! ! 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 == {{{ #!fortran ! 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 == {{{ #!fortran 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 == {{{ #!fortran ! 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 == {{{ #!fortran 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 == {{{ #!fortran ! ! 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 == {{{ #!fortran !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 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 == {{{ #!fortran 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 }}}