Context Navigation

Changes between Version 2 and Version 3 of fortran

Timestamp:: Oct 17, 2016, 8:50:12 PM (9 years ago)
Author:: Leon Kos
Comment:: Full set

Legend:

: Unmodified
: Added
: Removed
: Modified

fortran

-                      v2
+                      v3
 }}}
+== Console29/Console29/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
+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)
+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
+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
+}}}
+== Console3/Console3/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
+}}}
+== Console30/Console30/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
+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
+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
+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.
+AUX(6,I)=0.
+      IREC=0
+      H=H+H
+      IHLF=-1
+      ISTEP=0
+      IEND=0
+!
+!
+!     START OF A RUNGE-KUTTA STEP
+IF((X+H-XEND)*H)7,6,5
+H=XEND-X
+IEND=1
+!
+!     RECORDING OF INITIAL VALUES OF THIS STEP
+CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)
+      IF(PRMT(5))40,8,40
+ITEST=0
+ISTEP=ISTEP+1
+!     START OF INNERMOST RUNGE-KUTTA LOOP
+      J=1
+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
+AUX(6,I)=AUX(6,I)+R2-CJ*R1
+      IF(J-4)12,15,15
+J=J+1
+      IF(J-3)13,14,13
+X=X+.5*H
+CALL FCT(X,Y,DERY)
+      GOTO 10
+!     END OF INNERMOST RUNGE-KUTTA LOOP
+!     TEST OF ACCURACY
+IF(ITEST)16,16,20
+!
+!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY
+DO 17 I=1,NDIM
+AUX(4,I)=Y(I)
+      ITEST=1
+      ISTEP=ISTEP+ISTEP-2
+IHLF=IHLF+1
+      X=X-H
+      H=.5*H
+      DO 19 I=1,NDIM
+      Y(I)=AUX(1,I)
+      DERY(I)=AUX(2,I)
+AUX(6,I)=AUX(3,I)
+      GOTO 9
+!
+!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE
+IMOD=ISTEP/2
+      IF(ISTEP-IMOD-IMOD)21,23,21
+CALL FCT(X,Y,DERY)
+      DO 22 I=1,NDIM
+      AUX(5,I)=Y(I)
+AUX(7,I)=DERY(I)
+      GOTO 9
+!
+!     COMPUTATION OF TEST VALUE DELT
+DELT=0.
+      DO 24 I=1,NDIM
+DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))
+      IF(DELT-PRMT(4))28,28,25
+!
+!     ERROR IS TOO GREAT
+IF(IHLF-10)26,36,36
+DO 27 I=1,NDIM
+AUX(4,I)=AUX(5,I)
+      ISTEP=ISTEP+ISTEP-4
+      X=X-H
+      IEND=0
+      GOTO 18
+!
+!     RESULT VALUES ARE GOOD
+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)
+DERY(I)=AUX(7,I)
+      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)
+      IF(PRMT(5))40,30,40
+DO 31 I=1,NDIM
+      Y(I)=AUX(1,I)
+DERY(I)=AUX(2,I)
+      IREC=IHLF
+      IF(IEND)32,32,39
+!
+!     INCREMENT GETS DOUBLED
+IHLF=IHLF-1
+      ISTEP=ISTEP/2
+      H=H+H
+      IF(IHLF)4,33,33
+IMOD=ISTEP/2
+      IF(ISTEP-IMOD-IMOD)4,34,4
+IF(DELT-.02*PRMT(4))35,35,4
+IHLF=IHLF-1
+      ISTEP=ISTEP/2
+      H=H+H
+      GOTO 4
+!
+!     RETURNS TO CALLING PROGRAM
+IHLF=11
+      CALL FCT(X,Y,DERY)
+      GOTO 39
+IHLF=12
+      GOTO 39
+IHLF=13
+CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)
+RETURN
+      END
+}}}
+== Console31/Console31/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/Console32/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
+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
+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
+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.
+AUX(6,I)=0.
+      IREC=0
+      H=H+H
+      IHLF=-1
+      ISTEP=0
+      IEND=0
+!
+!
+!     START OF A RUNGE-KUTTA STEP
+IF((X+H-XEND)*H)7,6,5
+H=XEND-X
+IEND=1
+!
+!     RECORDING OF INITIAL VALUES OF THIS STEP
+CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)
+      IF(PRMT(5))40,8,40
+ITEST=0
+ISTEP=ISTEP+1
+!     START OF INNERMOST RUNGE-KUTTA LOOP
+      J=1
+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
+AUX(6,I)=AUX(6,I)+R2-CJ*R1
+      IF(J-4)12,15,15
+J=J+1
+      IF(J-3)13,14,13
+X=X+.5*H
+CALL FCT(X,Y,DERY)
+      GOTO 10
+!     END OF INNERMOST RUNGE-KUTTA LOOP
+!     TEST OF ACCURACY
+IF(ITEST)16,16,20
+!
+!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY
+DO 17 I=1,NDIM
+AUX(4,I)=Y(I)
+      ITEST=1
+      ISTEP=ISTEP+ISTEP-2
+IHLF=IHLF+1
+      X=X-H
+      H=.5*H
+      DO 19 I=1,NDIM
+      Y(I)=AUX(1,I)
+      DERY(I)=AUX(2,I)
+AUX(6,I)=AUX(3,I)
+      GOTO 9
+!
+!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE
+IMOD=ISTEP/2
+      IF(ISTEP-IMOD-IMOD)21,23,21
+CALL FCT(X,Y,DERY)
+      DO 22 I=1,NDIM
+      AUX(5,I)=Y(I)
+AUX(7,I)=DERY(I)
+      GOTO 9
+!
+!     COMPUTATION OF TEST VALUE DELT
+DELT=0.
+      DO 24 I=1,NDIM
+DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))
+      IF(DELT-PRMT(4))28,28,25
+!
+!     ERROR IS TOO GREAT
+IF(IHLF-10)26,36,36
+DO 27 I=1,NDIM
+AUX(4,I)=AUX(5,I)
+      ISTEP=ISTEP+ISTEP-4
+      X=X-H
+      IEND=0
+      GOTO 18
+!
+!     RESULT VALUES ARE GOOD
+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)
+DERY(I)=AUX(7,I)
+      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)
+      IF(PRMT(5))40,30,40
+DO 31 I=1,NDIM
+      Y(I)=AUX(1,I)
+DERY(I)=AUX(2,I)
+      IREC=IHLF
+      IF(IEND)32,32,39
+!
+!     INCREMENT GETS DOUBLED
+IHLF=IHLF-1
+      ISTEP=ISTEP/2
+      H=H+H
+      IF(IHLF)4,33,33
+IMOD=ISTEP/2
+      IF(ISTEP-IMOD-IMOD)4,34,4
+IF(DELT-.02*PRMT(4))35,35,4
+IHLF=IHLF-1
+      ISTEP=ISTEP/2
+      H=H+H
+      GOTO 4
+!
+!     RETURNS TO CALLING PROGRAM
+IHLF=11
+      CALL FCT(X,Y,DERY)
+      GOTO 39
+IHLF=12
+      GOTO 39
+IHLF=13
+CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)
+RETURN
+      END
+}}}
+== Console33/Console33/Console33.f90 ==
+{{{
+#!fortran
+! ��襭�� �ࠢ����� ����ᮭ� ��⮤�� �ண���� (test)
+! �ᯮ���㥬 ���⥬� ��
+!
+      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
+format(5x,i5/(5x,f20.8))
+      write(6,3) j,al,bl,cl,ar,br,cr
+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)
+format(3x,I3,3x,3e15.5)
+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/Console35/Console35.f90 ==
+{{{
+#!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
+format(5x,i5/(5x,f20.8))
+      write(6,3) j,al,bl,cl,ar,br,cr
+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
+format(3x,I4,3x,e12.5)
+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
+}}}
+== Console4/Console4/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/Console5/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/Console6/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/Console7/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/Console8/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
+format (f8.3, f8.3, f8.3)
+return
+end
+}}}
+== Console9/Console9/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
+}}}
+== Eyler/Eyler/Eyler.f90 ==
+{{{
+#!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
+}}}
+== uskarenie-protona/Console36/Console36.f90 ==
+{{{
+#!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
+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,  &
+.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
+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
+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
+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.
+AUX(6,I)=0.
+      IREC=0
+      H=H+H
+      IHLF=-1
+      ISTEP=0
+      IEND=0
+!
+!
+!     START OF A RUNGE-KUTTA STEP
+IF((X+H-XEND)*H)7,6,5
+H=XEND-X
+IEND=1
+!
+!     RECORDING OF INITIAL VALUES OF THIS STEP
+CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)
+      IF(PRMT(5))40,8,40
+ITEST=0
+ISTEP=ISTEP+1
+!     START OF INNERMOST RUNGE-KUTTA LOOP
+      J=1
+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
+AUX(6,I)=AUX(6,I)+R2-CJ*R1
+      IF(J-4)12,15,15
+J=J+1
+      IF(J-3)13,14,13
+X=X+.5*H
+CALL FCT(X,Y,DERY)
+      GOTO 10
+!     END OF INNERMOST RUNGE-KUTTA LOOP
+!     TEST OF ACCURACY
+IF(ITEST)16,16,20
+!
+!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY
+DO 17 I=1,NDIM
+AUX(4,I)=Y(I)
+      ITEST=1
+      ISTEP=ISTEP+ISTEP-2
+IHLF=IHLF+1
+      X=X-H
+      H=.5*H
+      DO 19 I=1,NDIM
+      Y(I)=AUX(1,I)
+      DERY(I)=AUX(2,I)
+AUX(6,I)=AUX(3,I)
+      GOTO 9
+!
+!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE
+IMOD=ISTEP/2
+      IF(ISTEP-IMOD-IMOD)21,23,21
+CALL FCT(X,Y,DERY)
+      DO 22 I=1,NDIM
+      AUX(5,I)=Y(I)
+AUX(7,I)=DERY(I)
+      GOTO 9
+!
+!     COMPUTATION OF TEST VALUE DELT
+DELT=0.
+      DO 24 I=1,NDIM
+DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))
+      IF(DELT-PRMT(4))28,28,25
+!
+!     ERROR IS TOO GREAT
+IF(IHLF-10)26,36,36
+DO 27 I=1,NDIM
+AUX(4,I)=AUX(5,I)
+      ISTEP=ISTEP+ISTEP-4
+      X=X-H
+      IEND=0
+      GOTO 18
+!
+!     RESULT VALUES ARE GOOD
+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)
+DERY(I)=AUX(7,I)
+      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)
+      IF(PRMT(5))40,30,40
+DO 31 I=1,NDIM
+      Y(I)=AUX(1,I)
+DERY(I)=AUX(2,I)
+      IREC=IHLF
+      IF(IEND)32,32,39
+!
+!     INCREMENT GETS DOUBLED
+IHLF=IHLF-1
+      ISTEP=ISTEP/2
+      H=H+H
+      IF(IHLF)4,33,33
+IMOD=ISTEP/2
+      IF(ISTEP-IMOD-IMOD)4,34,4
+IF(DELT-.02*PRMT(4))35,35,4
+IHLF=IHLF-1
+      ISTEP=ISTEP/2
+      H=H+H
+      GOTO 4
+!
+!     RETURNS TO CALLING PROGRAM
+IHLF=11
+      CALL FCT(X,Y,DERY)
+      GOTO 39
+IHLF=12
+      GOTO 39
+IHLF=13
+CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)
+RETURN
+      END
+}}}