c-------------------------------------------------------------------- 
c                    Set Initial Conditions for N-body simulations
c                    Uses VERY simple method
c
c  
c-------------------------------------------------------------------- 
      INCLUDE 'nbody1.h'
      REAL*8    Ekin,Epot,Etot,yc(3)

                          ! open files
      Open(1,file='DATA/particles.dat',form='unformatted')
      Open(3,file='DATA/energies.dat')
      Open(8,file='DATA/analysis.dat')
      Open(11,file='DATA/initial.dat')


                          ! set parameters
      Np         = 200
      time       = 0.
      dt         = 3.d-5
      epsilon    = 1.d-3
      RunTime    = 10.
      eps2       = epsilon**2
      v0         = 0.3465
      R0         = 1.
      do i=1,3         ! offset for the center of snapshot
         yc(i)   = 0.
      EndDo 
            write(3,30)             ! write headers for files
            write(8,40)

      Call Initial(v0,R0)
      Call Energies(Ekin,Epot)
            Etot = Ekin+Epot
            write(3,20) time,iStep,Ekin,Epot,Etot,Ekin/Etot

      Call Analyze
      Call WriteSnapshot(yc)
      Call SaveMoment

      write(*,'(80("-"))')
      write (*,*)'-----------   INITIAL CONDITIONS      -------------'
      write(*,*) ' Number of particles=',Np,' Time=',time
      write(*,*) ' Time step          =',dt,' Step=',iStep
      write(*,*) ' Force softening    =',epsilon
      write(*,30)
      write(*,20) time,iStep,Ekin,Epot,Etot,Ekin/Etot
 20   format(f8.3,i8,6g12.5)
 30   format(3x,'Time',T13,'Step',T19,'Ekin',T30,'Epot',
     &          T42,'Etot',T53,'Ekin/Etot')
 40   format(3x,'Time',T13,'Step',T19,'TotalM',
     &          T28,'MCentral',
     &          T37,'Rmax',T44,'Raverage',
     &          15(' Mass   '))
      Stop
      End
c-------------------------------------------------------------------- 
c               save current state of the system
      SUBROUTINE Initial(v0,R0)
c-------------------------------------------------------------------- 
      INCLUDE 'nbody1.h'

      Nseed = 121071
         Do i=1,Np
            Coords(10,i) =1./Np
 3          do k=1,3
               Coords(k,i) =R0*(2.*RANDD(Nseed)-1.)
            EndDo 
            rr = sqrt(Coords(1,i)**2+Coords(2,i)**2+Coords(3,i)**2)
            if(rr.gt.R0)goto 3
            do k=4,6
               Coords(k,i) =v0*GAUSS(Nseed)
            EndDo 
         EndDo 
 50      format(i8,10g11.4)
      Return
      End

C--------------------------------------
C		normal random numbers
      FUNCTION GAUSS(M)
C                    Uses RANDd() + 5 random numbers with correction
C                    Old code used in many simulations
C   N=1e8  sigma= 0.99556     mean= 0.10542E-03
c  sigma Frac(>sig) Frac(<-sig)   True       n(>)    n(<)
c  1.00 0.1589     0.1589     0.1587     15894145 15889594
c  1.50 0.6704E-01 0.6702E-01 0.6681E-01  6703736  6702430
c  2.00 0.2258E-01 0.2256E-01 0.2275E-01  2258038  2256030
c  2.50 0.5841E-02 0.5843E-02 0.6210E-02   584148   584264
c  3.00 0.1031E-02 0.1028E-02 0.1350E-02   103141   102783
c  4.00 0.5900E-06 0.6100E-06 0.3170E-04       59       61
C
C--------------------------------------
      X=0.
      DO  I=1,5
         X=X+RANDd(M)
      EndDo
      X2   =1.5491933*(X-2.5)
      GAUSS=X2*(1.-0.01*(3.-X2*X2))
      RETURN
      END
C--------------------------------------
C		normal random numbers
      FUNCTION GAUSS2(M)
C                        Uses ranndd for homogenous rand numbers
C                          sum of 12 random numbers + corrections
C                          It is  slow: very (10 times than GAUSS3)
C                          Quality is very good: upto 4 sigma
c  sigma Frac(>sig) Frac(<-sig)   True       n(>)    n(<)
c  1.00 0.1586     0.1586     0.1587     15857116 15857711
c  1.50 0.6678E-01 0.6680E-01 0.6681E-01  6678110  6679725
c  2.00 0.2277E-01 0.2275E-01 0.2275E-01  2277430  2274868
c  2.50 0.6197E-02 0.6189E-02 0.6210E-02   619732   618895
c  3.00 0.1327E-02 0.1319E-02 0.1350E-02   132713   131929
c  4.00 0.2664E-04 0.2727E-04 0.3170E-04     2664     2727
C                   
c-------------------------------------------------------------------- 
      X=0.
      DO  I=1,12
         X=X+Randd(M)
      EndDo
      X2 = X-6.0
      GAUSS2   =X2 *(1.-0.0045*(3.-X2*X2))
      RETURN
      END

C------------------------------------------------
C                                              random number generator
      FUNCTION RANDd(M)
C------------------------------------------------
      DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4
     +	/453815927,2147483648.,2147483647,536870912,131072,256,
     +	 16777216,4,16384,8388608,128/
      ML=M/K1*K1
      M1=(M-ML)*L1
      ML=M/K2*K2
      M2=(M-ML)*L2
      ML=M/K3*K3
      M3=(M-ML)*L3
      ML=M/K4*K4
      M4=(M-ML)*L4
      M5=KI-M
      IF(M1.GE.M5)M1=M1-KI-1
      ML=M+M1
      M5=KI-ML
      IF(M2.GE.M5)M2=M2-KI-1
      ML=ML+M2
      M5=KI-ML
      IF(M3.GE.M5)M3=M3-KI-1
      ML=ML+M3
      M5=KI-ML
      IF(M4.GE.M5)M4=M4-KI-1
      ML=ML+M4
      M5=KI-ML
      IF(LC.GE.M5)ML=ML-KI-1
      M=ML+LC
      RANDd=M/AM
      RETURN
      END