| Version 2 (modified by , 8 years ago) (diff) |
|---|
FORTRAN 95 examples
Console1/Console1/Console1.f90
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
Console10/Console10/Console10.f90
! Osc_rk_0.f90 ! ! Lab #3 ! metod Runge-Kuta. Oscillator. external fct, out real aux(8,4) common /a/ dt, k, n, vm, A, om, rl, x, y, r, u(1), u(2), u(3), u(4) real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/) integer :: nd = 2 real :: dt = 2*3.14159265 real :: t real :: c = 3.0e10 !B0=me*c*om/ez B0 = 200.0 E0 = 5000 eq = 4.8e-10 me = 9.1e-28 om = eq*B0/me*c A = 2.0 vm = A*om n = 1 k = 1 b = B0/B0 rl = c/om u(1) = vm/vm u(2) = vm/vm u(3) = x/rl u(4) = y/rl open (unit=2, file = '1_rk.dat', status='unknown') write(*,*) 'x=', u(3), 'vx=', u(1), 'y=', u(4), 'vy=', u(2) call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=', ih end subroutine fct (t, u, du) real u(3), du(3), u(4), du(4) common /a/ dt, k, n, vm, A, om, rl, x, y, r, u(1), u(2), u(3), u(4) du(1) = -u(3) du(2) = u(4) du(3) = u(1) du(4) = u(2) return end subroutine out (t, u, du, ih, nd, pr) common /a/ dt, k, n, om, vm, A, rl, x, y, r, u(1), u(2), u(3), u(4) if (t.ge.k*(dt/20.)) then write (4,1) t, u(1), u(2), u(3), u(4) k = k + 1 else end if 1 format (f8.3, f8.3, f8.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
Console11/Console11/Console11.f90
! Console11.f90 ! ! FUNCTIONS: ! Console11 - Entry point of console application. ! !**************************************************************************** ! ! PROGRAM: Console11 ! ! PURPOSE: Entry point for the console application. ! !**************************************************************************** program Console11 implicit none ! Variables ! Body of Console11 print *, 'Hello World' end program Console11
Console12/Console12/Console12.f90
! Console12.f90 ! ! FUNCTIONS: ! Console12 - Entry point of console application. ! !**************************************************************************** ! ! PROGRAM: Console12 ! ! PURPOSE: Entry point for the console application. ! !**************************************************************************** program Console12 implicit none ! Variables ! Body of Console12 print *, 'Hello World' end program Console12
Console13/Console13/Console13.f90
integer :: n=1000 real(8) :: dt, t, vm, b, rl, om, x, y, vx, vy, B0 real(8) :: eq=-4.8e-10, me=9.1e-28, c = 3.0e10 !B0=me*c*om/ez B0 = 200.0 !E0 = 5000 !A = 1e-30 x=1.0 y=1.0 vx=1000.0 vy=1000.0 a=100 !b = B0/B0 !om = eq*B0/me*c !vm = A*om !print *, vm !rl = c/om !vx = vm/vm !vy = vm/vm !x = x/rl !y = y/rl !print *, rl t = 1.0e-8 dt = t/float(n) open (1,file='1',status='unknown') do i=1,n vx = vx+(B0*eq)/(c*me)*vy*dt vy = vy-(B0*eq)/(c*me)*vx*dt x = x+vx*dt y = y+vy*dt write(1,*) x, vx end do pause end program
Console15/Console15/Console15.f90
x = (/1,5/) print *, x end
Console16/Console16/Console16.f90
integer :: n=1000 real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0 x = 1.0 y = 1.0 a = 0.5 vx = 1.0e09 vy = 1.0e09 E0=0.0 B0 = 200.0 b=(2*3.14*6)/(1.0e09) !B0 = 1e-25 !E0 = 5000 !A = 1e-30 !b = B0/B0 !om = eq*B0/me*c !vm = A*om !print *, vm !rl = c/om !vx = vm/vm !vy = vm/vm !x = x/rl !y = y/rl !print *, rl t = 1.00e-08 dt = t/float(n) open (1,file='1',status='unknown') do i=1,n Bz = B0*(1+a*y) !Bz = B0 x = x + vx*dt y = y + vy*dt !vx = vx + (vy*eq*Bz)/(c*me)*dt !vy = vy - (vx*eq*Bz)/(c*me)*dt if(i*dt.gt.b) then E0 = 5.0 end if vx = vx +(eq*E0)/(me)*dt+(vy*eq*Bz)/(c*me)*dt vy = vy + (eq*E0)/(me)*dt-(vx*eq*Bz)/(c*me)*dt ! x = x + vx*dt !y = y + vy*dt !print *, vx, vy write (1,*) x, y end do pause end program
Console17/Console17/Console17.f90
integer :: N = 1500 real(8) :: x = 2.5, y = 0.5, B0 = 500.0, We = 4.0054e-08, b = 2.0, pi = 3.14159, m = 9.1e-28, e = -4.8e-10, c = 2.99e10 real(8) vx, vy, v, dt, t v = sqrt(2*We/m) vx = sin(pi/4.0)*v vy = cos(pi/4.0)*v t = 5e-9 dt=t/float(N) !Norma !tau = t*10e8 !me = m*1e28 !qe = e*1e10 !vl = c*1e-10 !ve = v/c !vx = sin(pi/4.0)*v !vy = cos(pi/4.0)*v !dt = tau/float(N) open (1,file='mov',status='replace') do i=1,N Bz = B0*(1-b**2*(x**2+y**2)) x = x + vx*dt y = y + vy*dt vx = vx + vy/m*Bz*e/c*dt vy = vy - vx/m*Bz*e/c*dt print *, x, y pause write(1,*) y, x end do end program
Console18/Console18/Console18.f90
! Osc_rk_0.f90 ! ! Lab #3 ! metod Runge-Kuta. Oscillator. external fct, out real aux(8,2) common /a/ b, lam, f0, dt, k, n real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/) real :: u(2) = (/0., 0.5/) real :: du(2) = (/0.5, 0.5/) real :: b=0.003 integer :: nd = 2 real :: dt = 2*3.14159265 ! real :: lam = 0.1 real :: f0 = 0.00 n = 1 k = 1 open (unit=2, file = 'osc_rk.dat', status='new') write(*,*) 'x=', u(1), 'v=', u(2) pause call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=', ih stop end subroutine fct (t, u, du) real u(2), du(2), lam common /a/ b, lam, f0, dt, k, n du(1) = -u(2)-2.*b*u(1) du(2) = u(1) return end subroutine out (t, u, du, ih, nd, pr) real u(2), du(2), pr(5), lam common /a/ b, lam, f0, dt, k, n if (t.ge.k*(dt/20.)) then write (2,1) t,u(1),u(2) ! k = k + 1 else end if 1 format (f8.3, f8.3, f8.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
Console19/Console19/Console19.f90
! Osc_rk_0.f90 ! ! Lab #3 ! metod Runge-Kuta. Oscillator. external fct, out real aux(8,2) common /a/ b, lam, f0, dt, k, n, c real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/) real :: u(2) = (/0., 0.5/) real :: du(2) = (/0.5, 0.5/) real :: b=0.003 integer :: nd = 2 real:: c = 2.0 real :: dt = 2*3.14159265 ! real :: lam = 0.1 real :: f0 = 0.00 n = 1 k = 1 open (unit=2, file = 'osc_rk.dat', status='new') write(*,*) 'x=', u(1), 'v=', u(2) pause call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=', ih stop end subroutine fct (t, u, du) real u(2), du(2), lam common /a/ b, lam, f0, dt, k, n, c du(1) = c*cos(0.5*t)-u(2)-2.*b*u(1) du(2) = u(1) return end subroutine out (t, u, du, ih, nd, pr) real u(2), du(2), pr(5), lam common /a/ b, lam, f0, dt, k, n, c if (t.ge.k*(dt/20.)) then write (2,1) t,u(1),u(2) ! k = k + 1 else end if 1 format (f8.3, f8.3, f8.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
Console2/Console2/Console2.f90
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
Console20/Console20/Console20.f90
! Console20.f90 ! ! FUNCTIONS: ! Console20 - Entry point of console application. ! !**************************************************************************** ! ! PROGRAM: Console20 ! ! PURPOSE: Entry point for the console application. ! !**************************************************************************** program Console20 implicit none ! Variables ! Body of Console20 print *, 'Hello World' end program Console20
Console21/Console21/Console21.f90
integer :: n=1000 real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0 x = 1.0 y = 1.0 a = 0.5 vx = 1.0e09 vy = 1.0e09 E0 = 0.0 B0 = 200.0 b = (2*3.14)/(1.0e09) t = 1.00e-08 dt = t/float(n) open (1,file='1',status='unknown') do i=1,n Bz = B0*(1+a*y) !Bz = B0 x = x + vx*dt y = y + vy*dt vx = vx + (vy*eq*Bz)/(c*me)*dt vy = vy - (vx*eq*Bz)/(c*me)*dt if(i*dt.gt.b) then E0 = 3.0 end if vx = vx + (vy*eq*Bz)/(c*me)*dt vy = vy + (eq*E0)/(me)*dt - (vx*eq*Bz)/(c*me)*dt x = x - vx*dt y = y + vy*dt !print *, vx, vy write (1,*) x, y end do pause end program
Console22/Console22/Console22.f90
integer :: n=1000 real(8) :: dt, t, x, y, vx, vy, v, E, W real(8) :: mp = 1.67e-24, q = 4.8e-10 x = 1.0 y = 1.0 d = 1.0 a = 30 W0 = 1000 v = sqrt(2*W0/mp) t = 1.00e-05 dt = t/float(n) open (1,file='1',status='unknown') do i=1,n vx = v*sin(a) vy = v*cos(a) x = x + vx*dt y = y + vy*dt E = (4*3.14*mp*y)/d W = (E*q)/ x !print *, vx, vy write (1,*) W, t end do pause end program
Console24/Console24/Console24.f90
! Lab #1 ! metod Runge-Kuta. Rezonans i avtorezonans. external fct, out real aux(8,2) common /a/ dt, k, n, g0, om, B0, B, a real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/) real :: u(2) = (/0., 0.5/) real :: du(2) = (/0.5, 0.5/) integer :: nd = 2 real :: dt = 2*3.14159265 real :: t, m0, B, B0, a real :: f0 = 0.00 f = 2.45e09 pi = 3.14159265 ez = 4.8e-10 m0 = 9.1e-28 c = 3.0e10 om = 2*pi*f E = 2. w = 1000. a = 0.09 g0 = ez*E/(m0*c*om) B0 = m0*c*om/ez B = B0*(1+a*t) b = B/B0 u(1) = pi u(2) = w/511000.+1 write (*,*) g0 pause n = 1 k = 1 open (unit=2, file = 'rez_avtorez.dat', status='unknown') call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=' ,ih end subroutine fct (t, u, du) real u(2), du(2), lam common /a/ dt, k, n, g0, om, B, B0, a du(1) = (a*t+1-u(2))/u(2)+g0*sqrt(1./(u(2)**2-1))*sin(u(1)) du(2) = -g0*sqrt(1.-1./u(2)**2)*cos(u(1)) !write (*,*) du(1), du(2) !pause return end subroutine out (t, u, du, ih, nd, pr) real u(2), du(2), pr(5), lam common /a/ dt, k, n, g0, om, B, B0, a if (t.ge.k*(dt/2.)) then write (2,1) t,u(1),(u(2)-1)*511. k = k + 1 else end if 1 format (3e10.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 37 IHLF=12 GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
Console26/Console26/Console26.f90
! Lab #1 ! metod Runge-Kuta. Rezonans i avtorezonans. external fct, out real aux(8,4) common /a/ dt, k, n, g0, om, B0, B, a, y, rl real :: pr(5)=(/0., 2*3.14159*200., 2*3.14159/50, .0001, .0/) real :: u(4) = (/0., 0., 0., 0./) real :: du(4) = (/0.25, 0.25, 0.25, 0.25/) integer :: nd = 4 real :: dt = 2*3.14159265 real :: t, m0, B, B0, a real :: f0 = 0.00 f = 2.45e09 pi = 3.14159265 ez = 4.8e-10 m0 = 9.1e-28 c = 3.0e10 om = 2*pi*f E0 = 5. E = E0 a = 0.00034 B0 = m0*c*om/ez g0 = E0/B0 B = B0 b = B/B0 rl = c/om u(1) = 0. u(2) = 0. u(3) = 0. u(4) = 0. write (*,*) g0, B0, b, E, rl pause n = 1 k = 1 open (unit=2, file = 'rez_avtorez.dat', status='unknown') call rkgs(pr, u, du, nd, ih, fct, out, aux) write(2,*) 'ih=' ,ih end subroutine fct (t, u, du) real u(4), du(4), lam common /a/ dt, k, n, g0, om, B, B0, a, y, rl y=sqrt(1+u(1)**2+u(2)**2) du(1) = -g0*cos(t)- u(2)*(1+a*t)/y du(2) = -g0*sin(t)+u(1)*(1+a*t)/y du(3) = u(1)/y du(4) = u(2)/y !write (*,*) du(1), du(2), du(3), du(4) !pause return end subroutine out (t, u, du, ih, nd, pr) real u(4), du(4), pr(5), lam common /a/ dt, k, n, g0, om, B, B0, a, y, rl if (t.ge.k*(dt)) then write (2,1) t,(y-1)*511. ,u(3)*rl,u(4)*rl k = k + 1 else end if 1 format (5e12.3) return end SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ! ! DIMENSION Y(4),DERY(4),AUX(8,4),A(4),B(4),C(4),PRMT(5) DO 1 I=1,NDIM 1 AUX(8,I)=.06666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) ! ! ERROR TEST IF(H*(XEND-X))38,37,2 ! ! PREPARATIONS FOR RUNGE-KUTTA METHOD 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 ! ! PREPARATIONS OF FIRST RUNGE-KUTTA STEP DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0. IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 ! ! ! START OF A RUNGE-KUTTA STEP 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 ! ! RECORDING OF INITIAL VALUES OF THIS STEP 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 ! START OF INNERMOST RUNGE-KUTTA LOOP J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GOTO 10 ! END OF INNERMOST RUNGE-KUTTA LOOP ! TEST OF ACCURACY 15 IF(ITEST)16,16,20 ! ! IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GOTO 9 ! ! IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GOTO 9 ! ! COMPUTATION OF TEST VALUE DELT 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 ! ! ERROR IS TOO GREAT 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GOTO 18 ! ! RESULT VALUES ARE GOOD 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 ! ! INCREMENT GETS DOUBLED 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GOTO 4 ! ! RETURNS TO CALLING PROGRAM 36 IHLF=11 CALL FCT(X,Y,DERY) GOTO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END
