C ______________________________ Selects and Scales PM particles C Fof MMART -- initial conditions C December 1998 A. Klypin (aklypin@nmsu.edu) C C Scales internal PM coordinates and velocities C to coordinates in Mpc/h and km/s C C NROW = number of particles in 1D C NGRID= number of cells in 1D INCLUDE 'PMparameters.h' REAL INPUT Character*70 FileASCII,Line,FileCatalog Logical Inside Nseed = 121071 C................................................................... C Read data and open files WRITE (*,'(A,$)') 'Enter File Name for ASCII data = ' READ (*,'(A)') FileASCII OPEN(17,FILE=FileASCII,STATUS='UNKNOWN') WRITE (*,'(A,$)') 'Enter File Name for Halo Catalog = ' READ (*,'(A)') FileCatalog OPEN(18,FILE=FileCatalog,STATUS='UNKNOWN') CALL RDTAPE CALL SetWeight write (*,*) ' RDTAPE is done' WRITE (17,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv 100 FORMAT(1X,'Header=>',A45,/ + 1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ + 1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,/ + 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I6,/ + 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/ + 1x,' Omega_curvature=',F7.3) Box =INPUT(' Enter box size in comoving Mpc/h =') c write (*,*) c + ' Enter Limits for box (XYZ_min/max)', c + ' or Center and Radius (XYZ R 0 0)' c Xmin =INPUT(' X_min (Mpc/h)=') c Xmax =INPUT(' X_max (Mpc/h)=') c Ymin =INPUT(' Y_min (Mpc/h)=') c Ymax =INPUT(' Y_max (Mpc/h)=') c Zmin =INPUT(' Z_min (Mpc/h)=') c Zmax =INPUT(' Z_max (Mpc/h)=') Fraction =INPUT(' Enter fraction of particles to select=') c Fraction =1.0 BoxV = Box*100. ! Box size in km/s ScaleV = BoxV/AEXPN/NGRID ! scale factor for Velocities ScaleC = Box/NGRID ! scale factor for Coordinates Do i=1,16 ! skip header lines read(18,'(A)') Line write(*,'(A)') Line EndDo Radius =1. 200 read(18,*,err=300,end=300) Xcentr,Ycentr,Zcentr, + Vvx0,Vvy0,Vvz0,Am,Radius,Vrms,Vcirc c 200 read(18,*,err=300,end=300) Xcentr,Ycentr,Zcentr, c + Am,Vcir,Radius Radius =Radius/1000. ! Scale radius to Mpc/h ! Radius =Radius * 1.5 ! increase radius beyond the virial radius Xmin =Xcentr -Radius Xmax =Xcentr +Radius Ymin =Ycentr -Radius Ymax =Ycentr +Radius Zmin =Zcentr -Radius Zmax =Zcentr +Radius rad2 = Radius**2 xxmax =-1.e+9 xxmin = 1.e+9 vmax =-1.e+9 vmin = 1.e+9 xrand=0. Icount =0 inside = .true. ifile =1 If(Nspecies.eq.0)Then ! old case: 1 complete set Npages =NROW ! Number of pages N_in_last=NPAGE ! Number of particles in last page N_line =NROW N2 =N_line**2 Else ! multiple masses N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -NPAGE*(Npages-1) N_line=INT((lspecies(1))**0.33333+0.5) N2 =N_line**2 EndIf write (*,*) ' Pages=',Npages,' Species=',Nspecies write (*,*) ' Nparticles=',N_particles,' N_in_last=',N_in_last write (*,*) ' Particles_1d=',N_line DO IROW=1, Npages ! Loop over particle pages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last c write (*,*)' Move page=',IROW,' file=',ifile,' N=',In_page iL = NPAGE*(IROW-1) CALL GETROW(IROW,ifile) ! read in a page of particles DO IN=1, In_page ! Loop over particles X =ScaleC* (XPAR(IN)-1.) ! scale coordinates Y =ScaleC* (YPAR(IN)-1.) Z =ScaleC* (ZPAR(IN)-1.) If(Z.gt.Zmin .and. Z.le.Zmax)Then ! take particles only If(Y.gt.Ymin .and. Y.le.Ymax)Then ! inside the box If(X.gt.Xmin .and. X.le.Xmax)Then dd =(X-Xcentr)**2+(Y-Ycentr)**2+(Z-Zcentr)**2 inside = dd.lt.rad2.and.RANDd(Nseed).lt.fraction IF(inside)Then ! get the particle Ipart =IN+iL ! current particle number WPAR =iWeight(Ipart) ! particles weight Vxs=ScaleV* VX(IN) ! scale velocities Vys=ScaleV* VY(IN) Vzs=ScaleV* VZ(IN) Icount =Icount +1 ! count particles xxmax =MAX(xxmax,X,Y,Z) xxmin =MIN(xxmin,X,Y,Z) vmax =MAX(vmax,Vxs,Vys,Vzs) vmin =MIN(vmin,Vxs,Vys,Vzs) If(Nspecies.le.1)Then ! find lagrangian coords c Jpart =N_particles -Ipart+1 ! reverse the order of particles Jpart =Ipart ! direct order of particles k1 =INT((Jpart-1)/N2)+1 jj =Jpart-(k1-1)*N2 j1 =INT((jj-1)/N_line)+1 i1 =jj-(j1-1)*N_line write (17,250) IN,X,Y,Z,Vxs,Vys,Vzs,WPAR,i1,j1,k1 & ,Jpart,IROW ! write to a file Else write (17,250) IN,X,Y,Z,Vxs,Vys,Vzs,WPAR ! write to a file EndIf EndIf 250 Format(i8,3F9.4,3F9.2,f8.1,3i4,3i9) EndIf EndIf EndIf ENDDO ENDDO write (*,*) ' Scaled Coordinates were written to:', & FileASCII write (*,*) ' Number of particles =',Icount, & ' Fraction=',Fraction write (*,*) ' Limits:', Xmin,Xmax write (*,*) ' ', Ymin,Ymax write (*,*) ' ', Zmin,Zmax write (*,*) ' Min/Max of coordinates(Mpc/h)=',xxmin,xxmax write (*,*) ' Min/Max of velocities (km/s) =',vmin,vmax goto 200 300 write (*,*) ' <======= the End ========> ' END