!-------------------------------------------------------------------- ! Set Initial Conditions for N-body simulations ! Uses VERY simple method ! ! !-------------------------------------------------------------------- Program INITIALIZE use DataStructures REAL*8 Ekin,Epot,Etot,center(3),v0,R0,eccentr,am1,am2 Character*80 :: FileName='DATA/particlesA.dat' ! open files Open(3,file='DATA/energies.dat') Open(8,file='DATA/analysis.dat') ! set parameters Np = 500 write(3,30) ! write headers for files write(8,40) write(*,*) 'Np=',Np ALLOCATE (xc(Np),yc(Np), zc(Np)) ALLOCATE (vx(Np),vy(Np), vz(Np)) ALLOCATE (gx(Np),gy(Np), gz(Np)) ALLOCATE (LevPart(Np),aMass(Np)) LevPart =0. ! initialize arays: vector operations gx =0.; gy =0.; gz=0. v0 = 0.2755 R0 = 1. Call Initial(v0,R0) ! eccentr =0.9 ; am1 =1.; am2 =1.; R0 =1. ! Call Initial2(eccentr,am1,am2,R0) Call Energies(Ekin,Epot) Etot = Ekin+Epot write(3,20) time,iStep,Ekin,Epot,Etot,Ekin/Etot,dtstep,eps,iMultiple,aStep,Alph indx =0 Call Analyze(indx,center) Call WriteSnapshot(center) Call SaveMoment('DATA/particlesA.dat') write(*,'(80("-"))') write (*,*)'----------- INITIAL CONDITIONS -------------' write(*,*) ' Number of particles=',Np,' Time=',time write(*,*) ' Time step =',dtstep,' Step=',iStep write(*,*) ' Force softening =',eps write(*,30) write(*,20) time,iStep,Ekin,Epot,Etot,Ekin/Etot 20 format(f8.3,i8,2x,3g14.7,f7.4,1P,2g10.2,i5,3x,2g10.2) 30 format(3x,'Time',T13,'Step',T20,'Ekin',T33,'Epot', & T47,'Etot',T59,'Ekin/Etot',5x,'dt',7x,'eps',4x,' Multiple',4x,'aStep',5x,'Alph') 40 format(3x,'Time',T13,'Step',T19,'TotalM', & T28,'MCentral', & T37,'Rmax',T44,'Raverage', & 15(' Mass ')) Stop End PROGRAM INITIALIZE !-------------------------------------------------------------------- ! Generate initial conditions for ! homogeneous sphere of radius R0 ! RMS velocities are equal to V0 SUBROUTINE Initial(v0,R0) !-------------------------------------------------------------------- use DataStructures REAL*8 v0,R0 Nseed = 121071 Do i=1,Np aMass(i) =1./Np 3 xc(i) =R0*(2.*RANDD(Nseed)-1.) yc(i) =R0*(2.*RANDD(Nseed)-1.) zc(i) =R0*(2.*RANDD(Nseed)-1.) rr = sqrt(xc(i)**2+yc(i)**2+zc(i)**2) if(rr.gt.R0)goto 3 vx(i) =v0*GAUSS(Nseed) vy(i) =v0*GAUSS(Nseed) vz(i) =v0*GAUSS(Nseed) EndDo Return End SUBROUTINE Initial !-------------------------------------------------------------------- ! Generate initial conditions for two particels in keplerian motion ! masses are m1 m2 ! R0 = x1 is the distance of the first mass from the center ! eccentricity eccentr SUBROUTINE Initial2(eccentr,am1,am2,R0) !-------------------------------------------------------------------- use DataStructures REAL*8 eccentr,am1,am2,R0,Mtot If(Np /= 2) stop ' Wrong number of particles. It must be Np=2' Mtot = am1+am2 ! total mass of the system R = Mtot/am2*R0 ! this is the distance between particles V = sqrt(Mtot/R*(1.-eccentr)) ! '-' is for the appocenter position write(*,*) ' eccentr=', eccentr write(*,*) ' masses=', Mtot write(*,*) ' period=', 2.*3.141592653*sqrt(R**3/Mtot/(1.+eccentr)**3) xc(1) = R0 ; yc(1) = 0. ; zc(1) = 0. xc(2) = -am1/am2*R0 ; yc(2) = 0. ; zc(2) = 0. aMass(1) = am1; aMass(2) = am2 vx(1) = 0. ; vy(1) = am2/Mtot*V ; vz(1) = 0. vx(2) = 0. ; vy(2) = -am1/Mtot*V ; vz(2) = 0. write(*,10) aMass(1),xc(1),yc(1),zc(1),vx(1),vy(1),vz(1) write(*,10) aMass(2),xc(2),yc(2),zc(2),vx(2),vy(2),vz(2) 10 format(g12.3,T20,3g12.3,5x,3g12.3) Return End SUBROUTINE Initial2 !-------------------------------------- ! normal random numbers FUNCTION GAUSS(M) ! Uses RANDd() + 5 random numbers with correction ! Old code used in many simulations ! N=1e8 sigma= 0.99556 mean= 0.10542E-03 ! sigma Frac(>sig) Frac(<-sig) True n(>) n(<) ! 1.00 0.1589 0.1589 0.1587 15894145 15889594 ! 1.50 0.6704E-01 0.6702E-01 0.6681E-01 6703736 6702430 ! 2.00 0.2258E-01 0.2256E-01 0.2275E-01 2258038 2256030 ! 2.50 0.5841E-02 0.5843E-02 0.6210E-02 584148 584264 ! 3.00 0.1031E-02 0.1028E-02 0.1350E-02 103141 102783 ! 4.00 0.5900E-06 0.6100E-06 0.3170E-04 59 61 ! !-------------------------------------- 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 FUNCTION GAUSS !-------------------------------------- ! normal random numbers FUNCTION GAUSS2(M) ! Uses ranndd for homogenous rand numbers ! sum of 12 random numbers + corrections ! It is slow: very (10 times than GAUSS3) ! Quality is very good: upto 4 sigma ! sigma Frac(>sig) Frac(<-sig) True n(>) n(<) ! 1.00 0.1586 0.1586 0.1587 15857116 15857711 ! 1.50 0.6678E-01 0.6680E-01 0.6681E-01 6678110 6679725 ! 2.00 0.2277E-01 0.2275E-01 0.2275E-01 2277430 2274868 ! 2.50 0.6197E-02 0.6189E-02 0.6210E-02 619732 618895 ! 3.00 0.1327E-02 0.1319E-02 0.1350E-02 132713 131929 ! 4.00 0.2664E-04 0.2727E-04 0.3170E-04 2664 2727 ! !-------------------------------------------------------------------- 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 FUNCTION GAUSS2 !------------------------------------------------ ! random number generator FUNCTION RANDd(M) !------------------------------------------------ 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 FUNCTION RANDd