Changes between Version 2 and Version 3 of fortran


Ignore:
Timestamp:
Oct 17, 2016, 8:50:12 PM (4 years ago)
Author:
Leon Kos
Comment:

Full set

Legend:

Unmodified
Added
Removed
Modified
  • fortran

    v2 v3  
    15591559}}}
    15601560 
     1561
     1562== Console29/Console29/Console29.f90 ==
     1563{{{
     1564#!fortran
     1565! ! program El in a mirro trap,  parabolic approximation of the magnetic field
     1566! for a mirror trap. ECR
     1567
     1568real :: ez = 4.8e-10, c=3.0e10, me=9.1e-28, f=2.4e09, pi=3.14159265
     1569
     1570common/a/ gx(250),gy(250),x,y,z,ux,uy,uz,rl,gm,tp,pi,dt,cl2,al,dl,d
     1571
     1572
     1573open (unit=1, file = 'input_t.dat')
     1574open (unit=2, file = 'res.dat')
     1575!
     1576!
     1577  call cpu_time(start)
     1578!
     1579read(1,*) kt,kd,E,al,R,B00,x,y,z,W0,fi
     1580write(2,1) kt,kd,E,al,R,B00,x,y,z,W0,fi
     15811 format('kt='I7,2x,'kd=',I3,2x,'E(kV/cm)=',f4.1,2x,'al=',e10.3,2x,    &
     1582'R=',f5.3,2x,'B00(%)=',f8.5/'x(cm) =',f4.1,2x,'y(cm) =',f4.1,2x,    &
     1583'z(cm)=',f4.1,2x,'W0(eV)=',e10.3,2x,'fi=',f5.1)
     15842 format('g0=',e10.3,2x,'B0=',f6.1,2x,'cl=',e10.3,2x,'rl =',f6.2/    &
     1585'wmax=',e10.3,2x,'v0=',e10.3,2x,'w00=',e10.3,2x,'rc=',e10.3)
     1586
     1587pi2=2.0*pi
     1588
     1589DL=8.0
     1590D=12.0
     1591fi=fi/180.*pi
     1592
     1593v0=sqrt(1.-1./(1.+w0/511000.)**2)  ! ��������� �������� ��������� � �������� �
     1594
     1595write(*,*) v0*3.0e8                ! ��������� �������� ��������� � m/c
     1596U0=v0/sqrt(1.0-v0**2)              ! ��������� ������� ��������� � �������� mc
     1597W00=511000.*(sqrt(1.0+u0**2)-1.0)  ! ���� ��� ��������� �������
     1598om=2.0*pi*f
     1599rl=c/om                            ! �������������� ������ ������������� ��������
     1600zm=dl/2.0                          ! ������� ����������
     1601
     1602                                   ! ������ ��������� cl, ������������� �������� ���������� ����
     1603                                   ! �������, ����� �������� �������� ���������� ���������
     1604 cl=(zm/rl)/sqrt(R-1.0)
     1605 cl2=cl*cl                         
     1606 zm=zm/rl                           ! ����������
     1607
     1608E=E*1.0e05/3.0e04                  ! ��������� ��� ���� � ����
     1609g0=ez*E/(me*c*om)                  ! ������������ ��������� ��� ����
     1610
     1611B0=me*c*om/ez                      ! ����������� �������� �������� ���������� ���� - ECR ��� ��������� � ������ �����
     1612wmax=511.0*2.0*g0**(2.0/3.0)        ! ������������ �������� ������� ��� ��� � ���������� ������������� � ��������� ����
     1613rc=v0*3.0e10/om                    !
     1614 
     1615!    B0=875                        ! ��������� ���� � ������� � �������������� ������ �������
     1616!    B=875
     1617 write(2,2) g0, B0, cl, rl, wmax, v0, w00, rc
     1618     B=B/B0                     ! ���������������� �������� ���������� ����
     1619        x=x/rl           
     1620        y=y/rl
     1621        z=z/rl
     1622        ux=u0*sin(fi)               ! ����������� ������� �������� � ��������� ������ �������
     1623        uy=u0*cos(fi)
     1624        uz=0.                        ! � ����������� Z ��������� ������� ����� ����.
     1625        km=0. 
     1626        kp=0
     1627
     1628!    ������ ���� ��������������
     1629        m=250
     1630        dt=2.*pi/m
     1631         
     1632!    ������ ���������������� ����
     1633        do i=1,m
     1634        gx(i)=-g0*cos(pi2*(i-1)/m)*dt/2.0       
     1635        gy(i)=-g0*sin(pi2*(i-1)/m)*dt/2.0       
     1636        end do
     1637
     1638        write(2,*) 'wmax=', 2.*g0**(2./3.)*511.
     1639        wm=0
     1640
     1641!  ��������� ����
     1642        do k=1,kt
     1643        l=k-kp
     1644        tm=k*dt
     1645 
     1646        tp=-(1.+b00)*dt/2.         
     1647        call motion(l)
     1648if (k.eq.km+kd) then   ! ����� ����������� ����� kd ��������� �����
     1649        w=(gm-1.0)*511.0
     1650    write(2,3) tm/(2.*pi), x*rl,y*rl,z*rl, w
     1651
     1652        km=km+kd
     1653        else
     1654        end if
     1655if (k.eq.kp+m) then  ! ������� ��������� ����� ��� ��� ����
     1656        kp=kp+m
     1657        else
     1658        end if
     1659        end do
     16603 format(7f12.5)
     1661
     1662
     1663     call cpu_time(finish)
     1664        write(2,*) 'cpu-time =', finish-start
     1665        end
     1666                                             
     1667subroutine motion(l)   
     1668common/a/ gx(250),gy(250),x,y,z,ux,uy,uz,rl,gm,tp,pi,dt,cl2,al,dl,d
     1669
     1670! �������� ��������� ���� ������� ���������� ���� (���� 2)
     1671
     1672      r = sqrt(x*x + y*y)
     1673      bxp=-x*z/cl2
     1674      byp=-y*z/cl2
     1675      bzp=1.0+z*z/cl2-r*r/(2.*cl2)
     1676
     1677!      bxp=0.0
     1678!      byp=0.0
     1679!      bzp=1.0
     1680! ����� ������
     1681      u=sqrt(ux**2+uy**2)
     1682
     1683! �������� ��� ������������� ���� (���� 3)
     1684      uxm = ux  + gx(l)
     1685      uym = uy  + gy(l)   
     1686      uzm = uz
     1687
     1688      gmn = sqrt(1 + uxm*uxm + uym*uym + uzm*uzm)
     1689
     1690      tg = tp/gmn
     1691      tx = tg * bxp
     1692      ty = tg * byp
     1693      tz = tg * bzp
     1694
     1695      uxr = uxm + uym*tz - uzm*ty
     1696      uyr = uym + uzm*tx - uxm*tz
     1697      uzr = uzm + uxm*ty - uym*tx
     1698      gmn = sqrt(1 + uxm*uxm + uym*uym + uzm*uzm)
     1699
     1700      txyzr= 1.+ tx*tx + ty*ty + tz*tz
     1701
     1702      sx = 2*tx/txyzr
     1703      sy = 2*ty/txyzr
     1704      sz = 2*tz/txyzr
     1705
     1706      uxp = uxm + uyr*sz - uzr*sy
     1707      uyp = uym + uzr*sx - uxr*sz
     1708      uzp = uzm + uxr*sy - uyr*sx
     1709
     1710      ux = uxp + gx(l)
     1711      uy = uyp + gy(l)     
     1712      uz = uzp
     1713
     1714      gm = sqrt (1. + ux**2 + uy**2 + uz**2)
     1715      gt=dt/gm
     1716      x = x + ux*gt
     1717      y = y + uy*gt
     1718      z = z + uz*gt
     1719
     1720      if (abs(z*rl).ge.dl/2.0.or.r*rl.ge.d/2.0) then ! ������� ���������
     1721                                                      ! ��������� �� ������ ������
     1722      write(*,*) 'out',x*rl,y*rl,z*rl
     1723      stop
     1724      else
     1725      endif
     1726
     1727      return
     1728      end
     1729
     1730
     1731}}}
     1732
     1733== Console3/Console3/Console3.f90 ==
     1734{{{
     1735#!fortran
     1736real x,dx, a, b
     1737open(1,file='y2.dat', status='unknown')
     1738a=-0.5
     1739b=0.5
     1740n=100
     1741dx=(b-a)/n
     1742do i=1,n
     1743    write(1,*) i*dx + a, 1/sqrt(1-(i*dx+a)**2)
     1744end do
     1745end
     1746}}}
     1747
     1748== Console30/Console30/Console30.f90 ==
     1749{{{
     1750#!fortran
     1751!  Osc_rk_0.f90
     1752!
     1753!  Lab #3
     1754!  metod Runge-Kuta. Oscillator.
     1755
     1756external fct, out
     1757real aux(8,2)
     1758common /a/  dt, k, n, vm, A, om
     1759
     1760real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)
     1761real :: u(2) = (/0., 0.5/)
     1762real :: du(2) = (/0.5, 0.5/)
     1763integer :: nd = 2
     1764real :: dt = 2*3.14159265
     1765! real :: lam = 0.1
     1766real :: t
     1767real :: om=1
     1768A=2.0
     1769vm=A*om
     1770n = 1
     1771k = 1
     1772u(1)=vm/vm
     1773u(2)=0
     1774 
     1775open (unit=2, file = 'osc_rk.dat', status='unknown')
     1776
     1777
     1778call rkgs(pr, u, du, nd, ih, fct, out, aux)
     1779
     1780write(2,*) 'ih=', ih
     1781
     1782
     1783end
     1784
     1785subroutine fct (t, u, du)
     1786
     1787real u(2), du(2), lam
     1788common /a/ dt, k, n, g, vm, A, om
     1789
     1790du(1) = -u(2)/(1+1.2*sin(t*2))               
     1791du(2) = u(1)
     1792                               
     1793return
     1794end 
     1795
     1796subroutine  out (t, u, du, ih, nd, pr)
     1797real u(2), du(2), pr(5), lam
     1798common /a/ dt, k, n, om, vm, A
     1799
     1800if (t.ge.k*(dt/20.)) then
     1801
     1802 write (2,1) t,u(1),u(2)  !
     1803  k = k + 1
     1804   else
     1805    end if
     1806
     18071 format (f8.3, f8.3, f8.3)
     1808
     1809return
     1810end
     1811
     1812
     1813      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
     1814!                                                                       
     1815!                                                                       
     1816      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)           
     1817      DO 1 I=1,NDIM                                                     
     1818    1 AUX(8,I)=.06666667*DERY(I)                                       
     1819      X=PRMT(1)                                                         
     1820      XEND=PRMT(2)                                                     
     1821      H=PRMT(3)                                                         
     1822      PRMT(5)=0.                                                       
     1823      CALL FCT(X,Y,DERY)                                               
     1824!                                                                       
     1825!     ERROR TEST                                                       
     1826      IF(H*(XEND-X))38,37,2                                             
     1827!                                                                       
     1828!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
     1829    2 A(1)=.5                                                           
     1830      A(2)=.2928932                                                     
     1831      A(3)=1.707107                                                     
     1832      A(4)=.1666667                                                     
     1833      B(1)=2.                                                           
     1834      B(2)=1.                                                           
     1835      B(3)=1.                                                           
     1836      B(4)=2.                                                           
     1837      C(1)=.5                                                           
     1838      C(2)=.2928932                                                     
     1839      C(3)=1.707107                                                     
     1840      C(4)=.5                                                           
     1841!                                                                       
     1842!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                           
     1843      DO 3 I=1,NDIM                                                     
     1844      AUX(1,I)=Y(I)                                                     
     1845      AUX(2,I)=DERY(I)                                                 
     1846      AUX(3,I)=0.                                                       
     1847    3 AUX(6,I)=0.                                                       
     1848      IREC=0                                                           
     1849      H=H+H                                                             
     1850      IHLF=-1                                                           
     1851      ISTEP=0                                                           
     1852      IEND=0                                                           
     1853!                                                                     
     1854!                                                                   
     1855!     START OF A RUNGE-KUTTA STEP                                       
     1856    4 IF((X+H-XEND)*H)7,6,5                                             
     1857    5 H=XEND-X                                                         
     1858    6 IEND=1                                                           
     1859!                                                                     
     1860!     RECORDING OF INITIAL VALUES OF THIS STEP                         
     1861    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                               
     1862      IF(PRMT(5))40,8,40                                               
     1863    8 ITEST=0                                                           
     1864    9 ISTEP=ISTEP+1                                                     
     1865!     START OF INNERMOST RUNGE-KUTTA LOOP                               
     1866      J=1                                                               
     1867   10 AJ=A(J)                                                           
     1868      BJ=B(J)                                                           
     1869      CJ=C(J)                                                           
     1870      DO 11 I=1,NDIM                                                   
     1871      R1=H*DERY(I)                                                     
     1872      R2=AJ*(R1-BJ*AUX(6,I))                                           
     1873      Y(I)=Y(I)+R2                                                     
     1874      R2=R2+R2+R2                                                       
     1875   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                       
     1876      IF(J-4)12,15,15                                                   
     1877   12 J=J+1                                                             
     1878      IF(J-3)13,14,13                                                   
     1879   13 X=X+.5*H                                                         
     1880   14 CALL FCT(X,Y,DERY)                                               
     1881      GOTO 10                                                           
     1882!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
     1883!     TEST OF ACCURACY                                                 
     1884   15 IF(ITEST)16,16,20                                                 
     1885!                                                                       
     1886!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
     1887   16 DO 17 I=1,NDIM                                                   
     1888   17 AUX(4,I)=Y(I)                                                     
     1889      ITEST=1                                                           
     1890      ISTEP=ISTEP+ISTEP-2                                               
     1891   18 IHLF=IHLF+1
     1892      X=X-H                                                             
     1893      H=.5*H                                                           
     1894      DO 19 I=1,NDIM                                                   
     1895      Y(I)=AUX(1,I)                                                     
     1896      DERY(I)=AUX(2,I)                                                 
     1897   19 AUX(6,I)=AUX(3,I)                                                 
     1898      GOTO 9                                                           
     1899!                                                                       
     1900!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
     1901   20 IMOD=ISTEP/2                                                     
     1902      IF(ISTEP-IMOD-IMOD)21,23,21                                       
     1903   21 CALL FCT(X,Y,DERY)                                               
     1904      DO 22 I=1,NDIM                                                   
     1905      AUX(5,I)=Y(I)                                                     
     1906   22 AUX(7,I)=DERY(I)                                                 
     1907      GOTO 9                                                           
     1908!                                                                       
     1909!     COMPUTATION OF TEST VALUE DELT                                   
     1910   23 DELT=0.                                                           
     1911      DO 24 I=1,NDIM                                                   
     1912   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
     1913      IF(DELT-PRMT(4))28,28,25                                         
     1914!                                                                       
     1915!     ERROR IS TOO GREAT                                               
     1916   25 IF(IHLF-10)26,36,36                                               
     1917   26 DO 27 I=1,NDIM                                                   
     1918   27 AUX(4,I)=AUX(5,I)                                                 
     1919      ISTEP=ISTEP+ISTEP-4                                               
     1920      X=X-H                                                             
     1921      IEND=0                                                           
     1922      GOTO 18                                                           
     1923!                                                                       
     1924!     RESULT VALUES ARE GOOD                                           
     1925   28 CALL FCT(X,Y,DERY)                                               
     1926      DO 29 I=1,NDIM                                                   
     1927      AUX(1,I)=Y(I)                                                     
     1928      AUX(2,I)=DERY(I)                                                 
     1929      AUX(3,I)=AUX(6,I)                                                 
     1930      Y(I)=AUX(5,I)                                                     
     1931   29 DERY(I)=AUX(7,I)                                                 
     1932      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                             
     1933      IF(PRMT(5))40,30,40                                               
     1934   30 DO 31 I=1,NDIM                                                   
     1935      Y(I)=AUX(1,I)                                                     
     1936   31 DERY(I)=AUX(2,I)                                                 
     1937      IREC=IHLF                                                         
     1938      IF(IEND)32,32,39                                                 
     1939!                                                                       
     1940!     INCREMENT GETS DOUBLED                                           
     1941   32 IHLF=IHLF-1                                                       
     1942      ISTEP=ISTEP/2                                                     
     1943      H=H+H                                                             
     1944      IF(IHLF)4,33,33                                                   
     1945   33 IMOD=ISTEP/2                                                     
     1946      IF(ISTEP-IMOD-IMOD)4,34,4                                         
     1947   34 IF(DELT-.02*PRMT(4))35,35,4                                       
     1948   35 IHLF=IHLF-1                                                       
     1949      ISTEP=ISTEP/2                                                     
     1950      H=H+H                                                             
     1951      GOTO 4                                                           
     1952!                                                                       
     1953!     RETURNS TO CALLING PROGRAM                                       
     1954   36 IHLF=11                                                           
     1955      CALL FCT(X,Y,DERY)                                               
     1956      GOTO 39                                                           
     1957   37 IHLF=12                                                           
     1958      GOTO 39                                                           
     1959   38 IHLF=13                                                           
     1960   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                               
     1961   40 RETURN                                                           
     1962      END                                                               
     1963
     1964
     1965}}}
     1966
     1967== Console31/Console31/Console31.f90 ==
     1968{{{
     1969#!fortran
     1970    program heat1
     1971
     1972    implicit none
     1973
     1974    ! Variables
     1975    real a,b,c,d,dx,dt,x_i,t_j
     1976    integer n,m,i,j
     1977   
     1978    real, allocatable:: U(:,:)
     1979   
     1980    open(10,file='result_out1.txt')
     1981    a=1
     1982    b=4
     1983    c=1
     1984    d=15
     1985    write(0,*) "vvedite melkost razbijenija po x"
     1986    read(*,*) n
     1987   
     1988   ! do i=1,n
     1989     
     1990     write(0,*) "vvedite melkost razbijenija po t"
     1991        read(*,*) m
     1992        dx=(b-a)/n
     1993         dt=(d-c)/m
     1994        allocate(U(n,m))
     1995       
     1996       
     1997       do j=1,m      ! ��������� �������
     1998        t_j=j*dt
     1999        U(1,j)=3
     2000        U(n,j)=1
     2001        end do
     2002       
     2003        do i=2,n-1   ! ��������� ������� (�������  ����� �� ������, �.�. ���� ������ ���������)
     2004         x_i=i*dx
     2005        U(i,1)=1
     2006        end do
     2007       
     2008       
     2009        ! �������� ����
     2010         do j=1,m-1        ! ���� �� �������
     2011         do i=2,n-1        ! ���� �� ������������
     2012         U(i,j+1)=U(i,j)+0.005*(dt/(dx)**2)*(U(i+1,j)-2*U(i,j)+U(i-1,j))
     2013         print *,i,j+1,U(i,j+1)
     2014         end do
     2015         end do
     2016   
     2017   
     2018        !print *, x_i
     2019    !end do
     2020   
     2021        !do j=1,m
     2022       
     2023   
     2024   
     2025    !print *, t_j
     2026   
     2027   
     2028    do i=1,n
     2029    write(10,*) i*dx,U(i,m/100),U(i,m/2),U(i,m)
     2030    enddo
     2031   
     2032
     2033    end program heat1
     2034
     2035}}}
     2036
     2037== Console32/Console32/Console32.f90 ==
     2038{{{
     2039#!fortran
     2040!  Osc_rk_0.f90
     2041!
     2042!  Lab #3
     2043!  metod Runge-Kuta. Oscillator.
     2044
     2045external fct, out
     2046real aux(8,4)
     2047common /a/  dt, k, n, vm, A, om, rl 
     2048
     2049real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)
     2050real :: u(3) = (/0., 0.5/)
     2051real :: du(3) = (/0.5, 0.5/)
     2052real :: u(4) = (/0., 0.5/)
     2053real :: du(4) = (/0.5, 0.5/)
     2054integer :: nd = 2
     2055real :: dt = 2*3.14159265
     2056! real :: lam = 0.1
     2057real :: t
     2058real ::
     2059c = 3.0e10
     2060B0 = 200.0
     2061E0 = 5000
     2062eq = 4.8e-10
     2063me = 9.1e-28
     2064om = (eq*B0)/(me*c)
     2065A = 2.0
     2066vm = A*om
     2067n = 1
     2068k = 1
     2069b = B0/B0
     2070rl = c/om
     2071u(1) = vm/c
     2072u(2) = vm/c
     2073u(3) = x/rl
     2074u(4) = y/rl
     2075
     2076
     2077 
     2078open (unit=2, file = '1_rk.dat', status='unknown')
     2079
     2080
     2081call rkgs(pr, u, du, nd, ih, fct, out, aux)
     2082
     2083write(2,*) 'ih=', ih
     2084
     2085
     2086end
     2087
     2088subroutine fct (t, u, du)
     2089
     2090real u(3), du(3), u(4), du(4)
     2091common /a/ dt, k, n, vm, A, om, rl 
     2092
     2093du(1) = -u(2)
     2094du(2) = u(1)
     2095du(3) = u(1)
     2096du(4) = u(2)             
     2097                               
     2098return
     2099end 
     2100
     2101subroutine  out (t, u, du, ih, nd, pr)
     2102real u(3), du(3), u(4), du(4), pr(5), lam
     2103common /a/ dt, k, n, om, vm, A
     2104
     2105if (t.ge.k*(dt/20.)) then
     2106
     2107 write (2,1) t, u(1), u(2), u(3), u(4)  !
     2108  k = k + 1
     2109   else
     2110    end if
     2111
     21121 format (f8.3, f8.3, f8.3)
     2113
     2114return
     2115end
     2116
     2117
     2118      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
     2119!                                                                       
     2120!                                                                       
     2121      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)           
     2122      DO 1 I=1,NDIM                                                     
     2123    1 AUX(8,I)=.06666667*DERY(I)                                       
     2124      X=PRMT(1)                                                         
     2125      XEND=PRMT(2)                                                     
     2126      H=PRMT(3)                                                         
     2127      PRMT(5)=0.                                                       
     2128      CALL FCT(X,Y,DERY)                                               
     2129!                                                                       
     2130!     ERROR TEST                                                       
     2131      IF(H*(XEND-X))38,37,2                                             
     2132!                                                                       
     2133!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
     2134    2 A(1)=.5                                                           
     2135      A(2)=.2928932                                                     
     2136      A(3)=1.707107                                                     
     2137      A(4)=.1666667                                                     
     2138      B(1)=2.                                                           
     2139      B(2)=1.                                                           
     2140      B(3)=1.                                                           
     2141      B(4)=2.                                                           
     2142      C(1)=.5                                                           
     2143      C(2)=.2928932                                                     
     2144      C(3)=1.707107                                                     
     2145      C(4)=.5                                                           
     2146!                                                                       
     2147!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                           
     2148      DO 3 I=1,NDIM                                                     
     2149      AUX(1,I)=Y(I)                                                     
     2150      AUX(2,I)=DERY(I)                                                 
     2151      AUX(3,I)=0.                                                       
     2152    3 AUX(6,I)=0.                                                       
     2153      IREC=0                                                           
     2154      H=H+H                                                             
     2155      IHLF=-1                                                           
     2156      ISTEP=0                                                           
     2157      IEND=0                                                           
     2158!                                                                     
     2159!                                                                   
     2160!     START OF A RUNGE-KUTTA STEP                                       
     2161    4 IF((X+H-XEND)*H)7,6,5                                             
     2162    5 H=XEND-X                                                         
     2163    6 IEND=1                                                           
     2164!                                                                     
     2165!     RECORDING OF INITIAL VALUES OF THIS STEP                         
     2166    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                               
     2167      IF(PRMT(5))40,8,40                                               
     2168    8 ITEST=0                                                           
     2169    9 ISTEP=ISTEP+1                                                     
     2170!     START OF INNERMOST RUNGE-KUTTA LOOP                               
     2171      J=1                                                               
     2172   10 AJ=A(J)                                                           
     2173      BJ=B(J)                                                           
     2174      CJ=C(J)                                                           
     2175      DO 11 I=1,NDIM                                                   
     2176      R1=H*DERY(I)                                                     
     2177      R2=AJ*(R1-BJ*AUX(6,I))                                           
     2178      Y(I)=Y(I)+R2                                                     
     2179      R2=R2+R2+R2                                                       
     2180   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                       
     2181      IF(J-4)12,15,15                                                   
     2182   12 J=J+1                                                             
     2183      IF(J-3)13,14,13                                                   
     2184   13 X=X+.5*H                                                         
     2185   14 CALL FCT(X,Y,DERY)                                               
     2186      GOTO 10                                                           
     2187!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
     2188!     TEST OF ACCURACY                                                 
     2189   15 IF(ITEST)16,16,20                                                 
     2190!                                                                       
     2191!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
     2192   16 DO 17 I=1,NDIM                                                   
     2193   17 AUX(4,I)=Y(I)                                                     
     2194      ITEST=1                                                           
     2195      ISTEP=ISTEP+ISTEP-2                                               
     2196   18 IHLF=IHLF+1
     2197      X=X-H                                                             
     2198      H=.5*H                                                           
     2199      DO 19 I=1,NDIM                                                   
     2200      Y(I)=AUX(1,I)                                                     
     2201      DERY(I)=AUX(2,I)                                                 
     2202   19 AUX(6,I)=AUX(3,I)                                                 
     2203      GOTO 9                                                           
     2204!                                                                       
     2205!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
     2206   20 IMOD=ISTEP/2                                                     
     2207      IF(ISTEP-IMOD-IMOD)21,23,21                                       
     2208   21 CALL FCT(X,Y,DERY)                                               
     2209      DO 22 I=1,NDIM                                                   
     2210      AUX(5,I)=Y(I)                                                     
     2211   22 AUX(7,I)=DERY(I)                                                 
     2212      GOTO 9                                                           
     2213!                                                                       
     2214!     COMPUTATION OF TEST VALUE DELT                                   
     2215   23 DELT=0.                                                           
     2216      DO 24 I=1,NDIM                                                   
     2217   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
     2218      IF(DELT-PRMT(4))28,28,25                                         
     2219!                                                                       
     2220!     ERROR IS TOO GREAT                                               
     2221   25 IF(IHLF-10)26,36,36                                               
     2222   26 DO 27 I=1,NDIM                                                   
     2223   27 AUX(4,I)=AUX(5,I)                                                 
     2224      ISTEP=ISTEP+ISTEP-4                                               
     2225      X=X-H                                                             
     2226      IEND=0                                                           
     2227      GOTO 18                                                           
     2228!                                                                       
     2229!     RESULT VALUES ARE GOOD                                           
     2230   28 CALL FCT(X,Y,DERY)                                               
     2231      DO 29 I=1,NDIM                                                   
     2232      AUX(1,I)=Y(I)                                                     
     2233      AUX(2,I)=DERY(I)                                                 
     2234      AUX(3,I)=AUX(6,I)                                                 
     2235      Y(I)=AUX(5,I)                                                     
     2236   29 DERY(I)=AUX(7,I)                                                 
     2237      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                             
     2238      IF(PRMT(5))40,30,40                                               
     2239   30 DO 31 I=1,NDIM                                                   
     2240      Y(I)=AUX(1,I)                                                     
     2241   31 DERY(I)=AUX(2,I)                                                 
     2242      IREC=IHLF                                                         
     2243      IF(IEND)32,32,39                                                 
     2244!                                                                       
     2245!     INCREMENT GETS DOUBLED                                           
     2246   32 IHLF=IHLF-1                                                       
     2247      ISTEP=ISTEP/2                                                     
     2248      H=H+H                                                             
     2249      IF(IHLF)4,33,33                                                   
     2250   33 IMOD=ISTEP/2                                                     
     2251      IF(ISTEP-IMOD-IMOD)4,34,4                                         
     2252   34 IF(DELT-.02*PRMT(4))35,35,4                                       
     2253   35 IHLF=IHLF-1                                                       
     2254      ISTEP=ISTEP/2                                                     
     2255      H=H+H                                                             
     2256      GOTO 4                                                           
     2257!                                                                       
     2258!     RETURNS TO CALLING PROGRAM                                       
     2259   36 IHLF=11                                                           
     2260      CALL FCT(X,Y,DERY)                                               
     2261      GOTO 39                                                           
     2262   37 IHLF=12                                                           
     2263      GOTO 39                                                           
     2264   38 IHLF=13                                                           
     2265   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                               
     2266   40 RETURN                                                           
     2267      END                                                               
     2268
     2269
     2270}}}
     2271
     2272== Console33/Console33/Console33.f90 ==
     2273{{{
     2274#!fortran
     2275! ��襭�� �ࠢ����� ����ᮭ� ��⮤�� �ண���� (test)
     2276! �ᯮ���㥬 ���⥬� ��
     2277!
     2278      program poisn1
     2279real d,n0
     2280!   j - number of grid points, al - cr  - boudary coefficients
     2281!   q(1000),fi(1000) - charge and potentials at grid points
     2282!
     2283      common/psn/j,al,bl,cl,ar,br,cr,q(1000),fi(1000)
     2284!
     2285!     Ez(1000) - electric field strength at grid points
     2286!
     2287      dimension Ez(1000)
     2288!
     2289!  open files: p_in.dat - initial data;  in-dstr.dat - initial charge distribution
     2290!  fi_Ez.dat - values of potential and electric field strength at grid points
     2291!
     2292      open(5,file='p_in.dat',status='old')
     2293      open(6,file='in-dstr.dat',status='new')
     2294      open(7,file='fi_Ez.dat',status='new')
     2295!
     2296!   �⥭�� �� 䠩�� ������: ��᫮ 㧫�� ��⪨, ��ࠬ���� ��� ��砫���� �᫮���
     2297!                           number of grid points, other initial parameters
     2298!
     2299
     2300      read(5,1) j,al,bl,cl,ar,br,cr
     2301    1 format(5x,i5/(5x,f20.8))
     2302      write(6,3) j,al,bl,cl,ar,br,cr
     2303    3 format(5x,'j=',i5,//2x,6e10.3//)
     2304      pi=3.14159265
     2305      n0=10e+08
     2306      d=0.4            ! ⮫騭� ᫮� � ������      thick of a layer
     2307      jz=400             ! �������⢮ ���������� �祥�  number of filled cells
     2308      del=d/jz            ! 蠣 ��⪨ � ������   step (in meters)
     2309      !z=2
     2310     
     2311     
     2312          ! charge dencity
     2313     
     2314      rho=n0*abs(sin(2*pi*k/400))           
     2315       qz=rho
     2316!
     2317!  ��砫쭮� ����।������   spatial initial distribution
     2318!
     2319      do k=1,j
     2320      q(k)=0
     2321      if(k.ge.50.and.k.le.350) q(k)=n0*abs(sin(2*pi*k/400))
     2322      write(*,*) qz
     2323      end do
     2324!     
     2325!      �맮� ����ணࠬ�� ��襭�� �ࠢ����� ����ᮭ�   Poisson equation solver
     2326!     
     2327      call pnsolv
     2328!
     2329!    Calculation of electric field strength at grid points
     2330!
     2331      do k=2,j-1
     2332      Ez(k)=(fi(k-1)-fi(k+1))/(2.*del)
     2333      write(7,7) k,q(k),fi(k),Ez(k)     ! � �/�     (V/m)
     2334      end do
     2335!
     2336!    Test (analytic solution)
     2337!
     2338!      Ean=               
     2339      write(7,8) Ean                    ! � �/�     (V/m)
     2340    7 format(3x,I3,3x,3e15.5) 
     2341    8 format(3x,e15.5) 
     2342
     2343!
     2344      stop
     2345      end
     2346!
     2347!     SUBROUTINE PNSOLV FOR POISSON EQUATION SOLVING
     2348!
     2349      subroutine pnsolv
     2350      common/psn/j,al,bl,cl,ar,br,cr,q(1000),fi(1000)
     2351      dimension z(1000),y(1000)
     2352!
     2353!            �ண���� ��ࠢ� ������ (calculation  z and y coefficients)
     2354
     2355      j1=j+1
     2356      jm1=j-1             
     2357      cm=ar+br
     2358      z(j)=ar/cm
     2359      y(j)=cr/cm
     2360      do i=1,jm1
     2361          m=j1-i
     2362          cm=z(m)-2.
     2363          z(m-1)=-1./cm
     2364          y(m-1)=(q(m-1)-y(m))/cm
     2365      end do
     2366!
     2367!            �ண���� ᫥�� ���ࠢ� (calculation  fi(j)  at grid points)
     2368!
     2369      f0=(al*y(1)-cl)/(al-bl-al*z(1))
     2370      fi(1)=z(1)*f0+y(1)
     2371      do i=1,jm1
     2372          fi(i+1)=fi(i)*z(i+1)+y(i+1)
     2373      end do
     2374      return
     2375      end
     2376}}}
     2377
     2378== Console35/Console35/Console35.f90 ==
     2379{{{
     2380#!fortran
     2381!
     2382! 1D plasma simulation (demo version)
     2383!
     2384 program plasma_1D
     2385   real n,m0
     2386      common/psn/ nze,j,al,bl,cl,ar,br,cr,q(5000),f(5000),ze(100001)
     2387      common /E/ g0e,dt,rn,Uz(100001),Ez(5000)     
     2388      common/i/ nzi,qi(5000),zi(100001),Uzi(100001)
     2389  open(4,file='E(z).dat',status='new')
     2390  open(7,file='E(t).dat',status='new')
     2391  open(5,file='P-input.dat',status='old')
     2392  open(6,file='res.dat',status='new')
     2393  open(unit=2, file='input.dat')
     2394                                 call cpu_time(start)
     2395      read(2,*)  kt, dt, n, cfe        ! 
     2396      read(5,1) j,al,bl,cl,ar,br,cr
     2397    1 format(5x,i5/(5x,f20.8))
     2398      write(6,3) j,al,bl,cl,ar,br,cr
     2399    3 format(5x,'j=',i5,//2x,6e10.3//)
     2400!
     2401!  Input data
     2402!
     2403      pi=3.14159265
     2404      e=1.6E-19
     2405      m0=9.1E-31
     2406      c=3.0E8
     2407      om=sqrt(n*e*e/(m0*8.85e-12))
     2408      write(*,*) 'omega', om
     2409      Tosc=2.0*pi/om
     2410      rn=c/om
     2411      dt=dt*2.0*pi
     2412
     2413      d=0.1e-1            ! ⮫騭� ᫮� � ������      thick of a layer
     2414      jz=4000              ! �������⢮ ���������� �祥�  number of filled cells
     2415      del=d/jz            ! 蠣 ��⪨ � ������   step (in meters)
     2416
     2417      rho=-1.0e16         ! charge dencity
     2418      qz=rho*1.6e-19*(del**2)/8.85e-12
     2419      write(*,*) 'qz =', qz
     2420
     2421!      qz=0.0e-7
     2422
     2423      do i=1,4000
     2424      q(i)=0.         ! Charge of electrons at grid points
     2425      qi(i)=0.         ! Charge of ions at grid points
     2426      f(i)=0.           ! Potential at grid points
     2427      end do
     2428!
     2429      call ingen_ue   ! call subroutine for initial space electrons distribution
     2430      call ingen_ui   ! call subroutine for initial space ions distribution
     2431
     2432    7 format(3x,I4,3x,e12.5)
     2433    8 format(//)
     2434!
     2435      do i=1,4000
     2436      q(i)=(q(i)-qi(i))*qz   ! Charge at grid points (electrons+ions)
     2437      end do
     2438
     2439!      open(11,file='shift.dat',status='new')
     2440!      do i=1,4000
     2441!      write(11,*) i,q(i),qi(i)
     2442!      end do
     2443!      close(11)
     2444!      stop
     2445
     2446
     2447
     2448!
     2449      write(6,*) om,Tosc
     2450!
     2451      do i=1,nze
     2452      uz(i)=0.          ! initial velocities of electrons
     2453      end do
     2454      do i=1,nzi
     2455      uzi(i)=0.         ! initial velocities of ions
     2456      end do
     2457      kp0=50
     2458      kp=50
     2459!
     2460!           Cycle in time
     2461!
     2462      do k=1,kt
     2463      t=k*dt
     2464      call pnsolv        ! call Poisson's solver
     2465!
     2466      do i=2,3999
     2467      Ez(i)=(f(i-1)-f(i+1)) !/(2.0*del)   ! Calculation of electric field strength
     2468                                     ! at grid points
     2469      end do
     2470
     2471      open(11,file='shift.dat',status='new')
     2472      do i=1,4000
     2473      write(11,*) i,f(i),Ez(i)
     2474      end do
     2475      close(11)
     2476      stop
     2477
     2478
     2479      write(7,*) t/om/1.0e-9,Ez(1750)   ! write down time (nanosec) and E-field
     2480      do i=1,4000
     2481      q(i)=0.            ! Charge and potential at grid point = 0
     2482      f(i)=0.
     2483      end do
     2484!
     2485      call move         ! Integrating motion equation for electrons
     2486!
     2487!      open(11,file='shift.dat',status='new')
     2488!      do i=1,4000
     2489!      write(11,*) i,q(i),E(i)
     2490!      end do
     2491!      close(11)
     2492!      stop
     2493
     2494!
     2495      do i=1,4000
     2496      q(i)=(q(i)-qi(i))*qz   ! Charge at grid points (electrons+ions)
     2497      end do
     2498      end do
     2499!
     2500!                ! End of time cycle
     2501!
     2502!                Diagnostics: q(i), qi(i), Ez(i)
     2503!
     2504      open(3,file='El_dstr.dat',status='unknown')
     2505      open(13,file='Ions_dstr.dat',status='unknown')
     2506      do i=1,4000
     2507      write(3,*) (i-1)*0.0005,q(i)
     2508      end do
     2509      close(3)
     2510      do i=1,4000
     2511      write(13,*) (i-1)*0.0005,qi(i)
     2512      end do
     2513      close(13)
     2514      do i=1,3999
     2515      write(4,*) (i-1)*0.0005,Ez(i)
     2516      end do
     2517             call cpu_time(finish)
     2518write(6,*) 'calculaton time ',finish-start,' sec' 
     2519      end
     2520!
     2521!     SUBROUTINE PNSOLV FOR THE POISSON EQUATION SOLVING (Poisson's solver)
     2522!
     2523      subroutine pnsolv
     2524      common/psn/ nze,j,al,bl,cl,ar,br,cr,q(5000),f(5000),ze(100001)
     2525      dimension c(5000),s(5000)
     2526
     2527!
     2528!            DOWNWARDS SCANNING
     2529
     2530!      j1=j+1
     2531      j1=4001
     2532!      jm1=j-1
     2533      jm1=3999             
     2534      cm=ar+br
     2535      c(4000)=ar/cm
     2536      s(4000)=cr/cm
     2537      do i=1,jm1
     2538          m=j1-i
     2539          cm=c(m)-2.
     2540          c(m-1)=-1./cm
     2541          s(m-1)=(q(m-1)-s(m))/cm
     2542      end do
     2543!
     2544!            UPWARDS SCANNING
     2545!
     2546      f0=(al*s(1)-cl)/(al-bl-al*c(1))
     2547      f(1)=c(1)*f0+s(1)
     2548      do i=1,jm1
     2549          f(i+1)=f(i)*c(i+1)+s(i+1)
     2550      end do
     2551      return
     2552      end
     2553}}}
     2554
     2555== Console4/Console4/Console4.f90 ==
     2556{{{
     2557#!fortran
     2558real x,dx, a, b
     2559open(1,file='y2.dat', status='unknown')
     2560a=0
     2561b=1.14
     2562n=100
     2563dx=(b-a)/n
     2564do i=1,n
     2565    write(1,*) i*dx + a, sin(i*dx+a)*sin(2*i*dx+a)*sin(3*i*dx+a)
     2566end do
     2567end
     2568}}}
     2569
     2570== Console5/Console5/Console5.f90 ==
     2571{{{
     2572#!fortran
     2573program Simp
     2574
     2575integer n
     2576real y1, y2, y3, x, h, a, s
     2577n=100
     2578h=2/n
     2579do i=1,n
     2580   
     2581        x(i)=-1+h(i-1)
     2582
     2583        y1(i)=f1(x(i))
     2584        y2(i)=f2(x(i),h)
     2585        y3(i)=f3(x(i),h)
     2586       
     2587         a=h/6*(y1+4*y2+y3)
     2588         s=sum(a)
     2589end do
     2590
     2591contains
     2592real function f1(x)
     2593    real x
     2594    f1=x*(5-4*x)**(-0.5)
     2595  end function
     2596real function f2(x,h)
     2597    real x, h
     2598    f2=(x+h/2)*(5-4*(x+h/2))**(-0.5)
     2599  end function
     2600real function f3(x,h)
     2601    real x, h
     2602    f3=(x+h)*(5-4*(x+h))**(-0.5)
     2603
     2604 end function
     2605print *, s
     2606end program
     2607}}}
     2608
     2609== Console6/Console6/Console6.f90 ==
     2610{{{
     2611#!fortran
     2612program kolco
     2613
     2614integer :: n=1000
     2615real :: r1=1.0, r2=2.0, q=3.0, d=5.0, eq=-4.8e-10, m=9.1e-28 ! ������� ��������� � ���-�������
     2616real ::  dt, t, x, vx
     2617
     2618t=1.0e-9
     2619dt=t/float(n)
     2620vx=0
     2621x=2.5
     2622
     2623
     2624open (1,file='kolco',status='unknown')
     2625
     2626do i=1,n
     2627   E=q*x/(x**2+r1**2)**(1.5)- q*(5-x)/((5-x)**2+r2**2)**(1.5)
     2628   vx=vx+(eq*E/m)*dt
     2629   x=x+vx*dt
     2630   write(1,*) vx, x, i*dt
     2631end do
     2632
     2633end program
     2634}}}
     2635
     2636== Console7/Console7/Console7.f90 ==
     2637{{{
     2638#!fortran
     2639program kolco
     2640
     2641integer :: n=10000
     2642real(8) :: r1=1.0, r2=2.0, q=3.0, d=5.0, eq=4.8e-10, m=9.1e-28 ! ������� ��������� � ���-�������
     2643real(8) ::  dt, t, x, vx
     2644
     2645t=1.0e-6
     2646dt=t/float(n)
     2647vx=0
     2648x=2.5
     2649
     2650
     2651open (1,file='kolco',status='unknown')
     2652
     2653do i=1,n
     2654   E=q*x/sqrt((x**2+r1**2)**3)- q*(5-x)/sqrt(((5-x)**2+r2**2)**3)
     2655   vx=vx+(eq*E/m)*dt
     2656   x=x+vx*dt
     2657   write(1,*) x, i*dt, vx
     2658end do
     2659
     2660end program
     2661}}}
     2662
     2663== Console8/Console8/Console8.f90 ==
     2664{{{
     2665#!fortran
     2666external fct, out
     2667real aux(8,2)
     2668common /a/ dt, k, n, A, Vm, om
     2669real :: pr(5) = (/0.,2*3.14159*10., 2*3.14159/50, .0001, .0/)
     2670real :: du(2) = (/0.5, 0.5/)
     2671integer :: nd = 2
     2672real :: dt=2*3.14159265
     2673real :: t
     2674A = 2.0
     2675om = 1.0
     2676Vm = A*om
     2677u(1) = Vm/Vm
     2678u(2) = 0
     2679k = 1
     2680n=1
     2681
     2682open (unit=2, file = 'osc_rk0.dat' ,status='unknown')
     2683write(*,*) 'x=' , u(1), 'v=', u(2)
     2684pause
     2685call rkgs(pr, u, du, nd, ih, fct, out, aux)
     2686write(2,*) 'ih=', ih
     2687stop
     2688end
     2689!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2690
     2691subroutine fct (t, u, du)
     2692real u(2), du(2)
     2693common /a/ dt, k, n, A, Vm, om
     2694du(1) = -u(2)/(1+0.1*sin(om*dt))
     2695du(2) = u(1)
     2696return
     2697end
     2698
     2699
     2700subroutine out (t, u, du, ih, nd, pr)
     2701real u(2), du(2), pr(5)
     2702common /a/ dt, k, n, A, Vm, om
     2703
     2704if(t.ge.k*(dt/50.)) then
     2705write (2,1) t/om, u(1)*Vm, u(2)*A
     2706k=k+1
     2707else
     2708end if
     2709
     27101 format (f8.3, f8.3, f8.3)
     2711return
     2712end
     2713
     2714
     2715
     2716}}}
     2717
     2718== Console9/Console9/Console9.f90 ==
     2719{{{
     2720#!fortran
     2721!  Console9.f90
     2722!
     2723!  FUNCTIONS:
     2724!  Console9 - Entry point of console application.
     2725!
     2726
     2727!****************************************************************************
     2728!
     2729!  PROGRAM: Console9
     2730!
     2731!  PURPOSE:  Entry point for the console application.
     2732!
     2733!****************************************************************************
     2734
     2735    program Console9
     2736
     2737    implicit none
     2738
     2739    ! Variables
     2740
     2741    ! Body of Console9
     2742    print *, 'Hello World'
     2743
     2744    end program Console9
     2745
     2746}}}
     2747
     2748== Eyler/Eyler/Eyler.f90 ==
     2749{{{
     2750#!fortran
     2751real v,x,dt,t,dv,Fm
     2752open(2,file='y2.dat', status='unknown')
     2753q1=1.0e-9
     2754q2=1.0e-9
     2755r1=0.1
     2756r2=0.2
     2757d=0.5
     2758t0=d/v0
     2759n=100
     2760pi=3.141593
     2761E0=0.0
     2762v0=0.0
     2763dt=t0/n
     2764x=x/d
     2765write(*,*) dt
     2766open(1,file='xv.dat', status='unknown')
     2767Fm=q1/me/(vo**2)
     2768print*, Fm
     2769do k=1,n
     2770    t=k*dt
     2771    v=v+Fm*sin(2.0*pi*x)*dt
     2772    x=x+v*dt
     2773   
     2774write(1,*) t*t0/1.0e-9, v*v0, x*d, sin(2.0*pi*x)
     2775if(x.ge.1.0) then
     2776write(1,*) t*t0/1.0e-9, v*vo, x*d
     2777stop
     2778else
     2779end if
     2780end do
     2781end
     2782}}}
     2783
     2784== uskarenie-protona/Console36/Console36.f90 ==
     2785{{{
     2786#!fortran
     2787
     2788! ��� �����
     2789!������������ ������ ���������� ������� ����� ������������ �����
     2790!method Runge-Kutta
     2791!����� ����������� ��������� ���� ��������� ������� (������������ �������� !���������� ����)
     2792
     2793
     2794external fct, out
     2795common /a/ g0, dt, k, pp2b, gme, gmi, dbdz, pi, rl, d
     2796real aux(8,4),u(4),du(4),pr(5),n
     2797integer :: nd = 4
     2798real :: pi = 3.14159 
     2799real :: ez = 4.8e-10, c=3.0e10, me=9.1e-28, f=2.4e09
     2800
     2801open (unit=1, file = 'input.txt')
     2802open (unit=2, file = 'res.txt')
     2803read(1,*) kt,dt,n,d,B0,W0,dbdz
     2804write(2,1) kt,dt,n,d,B0,W0,dbdz
     28051 format('kt=',I5,2x,'dt=',f4.1,2x,'n(cm-3)=',e10.3,2x,'d(cm) =',f4.1,2x, &
     2806'B(Gs)=',f7.1/'W0(keV)=',e10.3,2x,'dBdz(cm-1) =',e10.3)
     2807
     2808   om=ez*B0/(me*c)
     2809     rl=c/om
     2810      E0=2.*pi*n*ez*d
     2811       g0=ez*E0/(me*om*c)
     2812        d=d/rl
     2813      dbdz=dbdz*rl   
     2814       b=B0/B0
     2815       dt=2*pi/10.0
     2816write(*,*) 'g0, om, b, dbdz', g0,om,b,dbdz
     2817pr(1)=0.
     2818pr(2)=20.0*pi*kt
     2819pr(3)=2.0*pi/100.
     2820pr(4)=1.0e-03
     2821pr(5)=0.
     2822
     2823v0=sqrt(1.-1./(1.+w0/511.)**2)
     2824U0=v0/sqrt(1.0-v0**2)
     2825W00=511.*(sqrt(1.0+u0**2)-1.0)
     2826
     2827
     2828gme = sqrt(1.0+u0**2) 
     2829
     2830write(*,*) 'v0, u0, w00, gme', v0, u0,w00,gme
     2831write(*,*) 'rl', rl
     2832
     2833u(1)=0.
     2834u(2)=0.000/rl
     2835u(3)=0.0
     2836u(4)=-0.04/rl
     2837dt=50.0*pi
     2838
     2839k = 1
     2840do i=1,4
     2841du(i)=0.25
     2842end do
     2843
     2844call rkgs(pr, u, du, nd, ih, fct, out, aux)
     2845write(*,*) 'ih=', ih
     2846end
     2847
     2848subroutine fct (t, u, du)
     2849real u(4), du(4)
     2850common /a/ g0, dt, k, pp2b, gme, gmi, dbdz, pi, rl, d 
     2851if(u(2)*rl.gt.210.) then
     2852write(*,*) 'dbzu =',1.0-dbdz*u(2),dbdz,u(2) 
     2853stop
     2854else
     2855end if
     2856
     2857if(gme**2 - 1.0 - u(1)**2.le.0.) then
     2858write(*,*) 'utr =', gme**2 - 1.0 - u(1)**2,gme,u(1),t
     2859stop
     2860else
     2861end if
     2862 utr=sqrt(gme**2 - 1.0 - u(1)**2)
     2863
     2864            du(1) = 0.5*utr**2/gme*dbdz                 ! norm
     2865             du(2) = u(1)/gme
     2866              du(3) = g0*(u(2)-u(4))/(d/2.)/1836.
     2867               du(4) = u(3)
     2868if(1.0-dbdz*u(2).le.0.) then
     2869write(*,*) 'B is zero?=',1.0-dbdz*u(2)
     2870stop
     2871else
     2872end if
     2873return
     2874end 
     2875
     2876subroutine  out (t, u, du, ih, nd, pr)
     2877real u(4), du(4), pr(4)
     2878common /a/ g0, dt, k, pp2b, gme, gmi, dbdz, pi, rl, d
     2879if (t.ge.k*2*pi/5.0) then
     2880 write (2,1) t/(2.*pi), u(2)*rl, u(4)*rl, (u(2)-u(4))*rl,  &
     2881             9.1e-28*1836.*(u(3)*3.0e10)**2/2./(1.6e-12)/1.0e06 
     2882k = k + 1
     2883
     2884if(abs(u(2)-u(4))*rl.gt.d/2.0*rl) then
     2885write(*,*) u(2)*rl, u(4)*rl, d/2.*rl
     2886write(*,*) 'out of layer'
     2887stop
     2888else
     2889end if
     2890
     2891 if(u(2)*rl.gt.200.) then
     2892write(*,*) u(2)*rl, u(4)*rl, d/2.*rl
     2893write(2,*) 'end of acceleration length'
     2894stop
     2895else
     2896end if
     2897
     2898else
     2899end if
     2900
     29011 format (f6.1,2x,6e11.3)
     2902return
     2903end
     2904
     2905
     2906      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
     2907!                                                                       
     2908!                                                                       
     2909      DIMENSION Y(4),DERY(4),AUX(8,4),A(4),B(4),C(4),PRMT(5)           
     2910      DO 1 I=1,NDIM                                                     
     2911    1 AUX(8,I)=.06666667*DERY(I)                                       
     2912      X=PRMT(1)                                                         
     2913      XEND=PRMT(2)                                                     
     2914      H=PRMT(3)                                                         
     2915      PRMT(5)=0.                                                       
     2916      CALL FCT(X,Y,DERY)                                               
     2917!                                                                       
     2918!     ERROR TEST                                                       
     2919      IF(H*(XEND-X))38,37,2                                             
     2920!                                                                       
     2921!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
     2922    2 A(1)=.5                                                           
     2923      A(2)=.2928932                                                     
     2924      A(3)=1.707107                                                     
     2925      A(4)=.1666667                                                     
     2926      B(1)=2.                                                           
     2927      B(2)=1.                                                           
     2928      B(3)=1.                                                           
     2929      B(4)=2.                                                           
     2930      C(1)=.5                                                           
     2931      C(2)=.2928932                                                     
     2932      C(3)=1.707107                                                     
     2933      C(4)=.5                                                           
     2934!                                                                       
     2935!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                           
     2936      DO 3 I=1,NDIM                                                     
     2937      AUX(1,I)=Y(I)                                                     
     2938      AUX(2,I)=DERY(I)                                                 
     2939      AUX(3,I)=0.                                                       
     2940    3 AUX(6,I)=0.                                                       
     2941      IREC=0                                                           
     2942      H=H+H                                                             
     2943      IHLF=-1                                                           
     2944      ISTEP=0                                                           
     2945      IEND=0                                                           
     2946!                                                                     
     2947!                                                                   
     2948!     START OF A RUNGE-KUTTA STEP                                       
     2949    4 IF((X+H-XEND)*H)7,6,5                                             
     2950    5 H=XEND-X                                                         
     2951    6 IEND=1                                                           
     2952!                                                                     
     2953!     RECORDING OF INITIAL VALUES OF THIS STEP                         
     2954    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                               
     2955      IF(PRMT(5))40,8,40                                               
     2956    8 ITEST=0                                                           
     2957    9 ISTEP=ISTEP+1                                                     
     2958!     START OF INNERMOST RUNGE-KUTTA LOOP                               
     2959      J=1                                                               
     2960   10 AJ=A(J)                                                           
     2961      BJ=B(J)                                                           
     2962      CJ=C(J)                                                           
     2963      DO 11 I=1,NDIM                                                   
     2964      R1=H*DERY(I)                                                     
     2965      R2=AJ*(R1-BJ*AUX(6,I))                                           
     2966      Y(I)=Y(I)+R2                                                     
     2967      R2=R2+R2+R2                                                       
     2968   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                       
     2969      IF(J-4)12,15,15                                                   
     2970   12 J=J+1                                                             
     2971      IF(J-3)13,14,13                                                   
     2972   13 X=X+.5*H                                                         
     2973   14 CALL FCT(X,Y,DERY)                                               
     2974      GOTO 10                                                           
     2975!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
     2976!     TEST OF ACCURACY                                                 
     2977   15 IF(ITEST)16,16,20                                                 
     2978!                                                                       
     2979!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
     2980   16 DO 17 I=1,NDIM                                                   
     2981   17 AUX(4,I)=Y(I)                                                     
     2982      ITEST=1                                                           
     2983      ISTEP=ISTEP+ISTEP-2                                               
     2984   18 IHLF=IHLF+1
     2985      X=X-H                                                             
     2986      H=.5*H                                                           
     2987      DO 19 I=1,NDIM                                                   
     2988      Y(I)=AUX(1,I)                                                     
     2989      DERY(I)=AUX(2,I)                                                 
     2990   19 AUX(6,I)=AUX(3,I)                                                 
     2991      GOTO 9                                                           
     2992!                                                                       
     2993!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
     2994   20 IMOD=ISTEP/2                                                     
     2995      IF(ISTEP-IMOD-IMOD)21,23,21                                       
     2996   21 CALL FCT(X,Y,DERY)                                               
     2997      DO 22 I=1,NDIM                                                   
     2998      AUX(5,I)=Y(I)                                                     
     2999   22 AUX(7,I)=DERY(I)                                                 
     3000      GOTO 9                                                           
     3001!                                                                       
     3002!     COMPUTATION OF TEST VALUE DELT                                   
     3003   23 DELT=0.                                                           
     3004      DO 24 I=1,NDIM                                                   
     3005   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
     3006      IF(DELT-PRMT(4))28,28,25                                         
     3007!                                                                       
     3008!     ERROR IS TOO GREAT                                               
     3009   25 IF(IHLF-10)26,36,36                                               
     3010   26 DO 27 I=1,NDIM                                                   
     3011   27 AUX(4,I)=AUX(5,I)                                                 
     3012      ISTEP=ISTEP+ISTEP-4                                               
     3013      X=X-H                                                             
     3014      IEND=0                                                           
     3015      GOTO 18                                                           
     3016!                                                                       
     3017!     RESULT VALUES ARE GOOD                                           
     3018   28 CALL FCT(X,Y,DERY)                                               
     3019      DO 29 I=1,NDIM                                                   
     3020      AUX(1,I)=Y(I)                                                     
     3021      AUX(2,I)=DERY(I)                                                 
     3022      AUX(3,I)=AUX(6,I)                                                 
     3023      Y(I)=AUX(5,I)                                                     
     3024   29 DERY(I)=AUX(7,I)                                                 
     3025      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                             
     3026      IF(PRMT(5))40,30,40                                               
     3027   30 DO 31 I=1,NDIM                                                   
     3028      Y(I)=AUX(1,I)                                                     
     3029   31 DERY(I)=AUX(2,I)                                                 
     3030      IREC=IHLF                                                         
     3031      IF(IEND)32,32,39                                                 
     3032!                                                                       
     3033!     INCREMENT GETS DOUBLED                                           
     3034   32 IHLF=IHLF-1                                                       
     3035      ISTEP=ISTEP/2                                                     
     3036      H=H+H                                                             
     3037      IF(IHLF)4,33,33                                                   
     3038   33 IMOD=ISTEP/2                                                     
     3039      IF(ISTEP-IMOD-IMOD)4,34,4                                         
     3040   34 IF(DELT-.02*PRMT(4))35,35,4                                       
     3041   35 IHLF=IHLF-1                                                       
     3042      ISTEP=ISTEP/2                                                     
     3043      H=H+H                                                             
     3044      GOTO 4                                                           
     3045!                                                                       
     3046!     RETURNS TO CALLING PROGRAM                                       
     3047   36 IHLF=11                                                           
     3048      CALL FCT(X,Y,DERY)                                               
     3049      GOTO 39                                                           
     3050   37 IHLF=12                                                           
     3051      GOTO 39                                                           
     3052   38 IHLF=13                                                           
     3053   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                               
     3054   40 RETURN                                                           
     3055      END                                                               
     3056
     3057}}}