C------------------------------------------------- C A.Klypin 25 Nov.1997 (NMSU) C Bound-Density-Maxima code: C Find maxima of mass in spheres of radius Rsearch1 C for particles in a periodic cube of size Box C Coordinates: 0-NGRID, Box =NGRID C All variables and constants are scaled to the Hubble constant of your C model. Only FINAL values are rescled to H=100km/s/Mpc C More detailed description of the algorithm and the code can be C found in file BDM_describe C-------------------------------------------------- C Arrays Label and Lst provide quick access to neibhor particles C Label(i,j,k) = number of the first particle in c the cell (i,j,k), if Label .ne. 0 C = 0 - empty cell C Cell = size of a cell for Label C Lst(i) = number of the next particle in the cell, C where i-th was found. C If Lst=0 - last point in the cell was found. C INCLUDE 'PMparameters.h' INCLUDE 'PMlists2.h' Real INPUT,Mhalo Real Overdens(Nrad),Mass(Nrad),MnMass(Nrad) COMMON /OVER/ Ovlim,Ovcnt DIMENSION Radsc(Nrad) C ------------------------------------------ Read data If(iHelp.eq.1)Call Help CALL RDTAPE If(iDebug.eq.1)write (*,*) ' RDTAPE is done' CALL SetWeights If(extras(100).gt.0.)Then Box =extras(100) Else Box =INPUT(' Enter box size in comoving Mpc/h =') EndIf WRITE (*,120) + AEXPN,ISTEP,NROW,NGRID write (*,*) ' Number of mass Species =',Nspecies write (*,*) ' Total number of particles=',lspecies(Nspecies) write (*,*) ' Number of small particles=',lspecies(1) write (*,*) ' Box size (Mpc/h) = ',Box 120 FORMAT(1X, + 5x,'Expansion Parameter= ',F8.3,' Step= ',i5,/ + 5x,'Nrow= ',I4,' Ngrid= ',I4) hsmall= hubble ! Hubble/100. Box = Box /hsmall ! Scale to real Mpc Xscale= Box/Ngrid ! Scale for comoving coordinates Vscale= 100.*hsmall*Xscale/AEXPN ! Scale for velocities Dscale= 2.746e+11*hsmall**2*(Box/NROW)**3 *weightSmall ! mass scale aMassOne= Om0*Dscale*hsmall write (*,'(6x,"Mass of smallest particle",/9x, + "in units M_sun/h is =",3x,g10.3)') + aMassOne write(*,'(6x,"first species=",2g11.3)')wspecies(1),weightSmall If(Nspecies.gt.1)Then write (*,'(6x,"Mass of largest particle",/9x, + "in units M_sun/h is =",3x,g10.3)') + Om0*Dscale*wspecies(Nspecies)*hsmall/(NGRID**3/FLOAT(NROW)**3) EndIf If(iHelp.eq.1)write(*,140) 140 format(10x,'Give average overdensity inside halo radius',/ +10x,'= 178 for CDM, Omg_0=1; =340 for LCDM Omg_0=0.3') Ovlim = OverdenVir() ! set virial overdensity Ovcnt = Ovlim ! Min. Center Overdensity for Halos ! scale factor to get number of virial particles aNvir = 4.1888*Ovlim*(NROW/float(NGRID))**3/weightSmall Amassl = AmassOne*5. ! Minimum halo mass in Msun/h If(iHelp.eq.1)write(*,150) 150 format(10x,'Code searches for maxima of mass inside spheres',/ +10x,'with radii ranging from small radius (second question) to',/ +10x,'large radius (first question). The radii may be equal.',/ +10x,'The search radii should be larger than the resolution ',/ +10x,'of the run. Code will fail if there is too few particles ',/ +10x,'inside the spheres.',/ +10x,'The smaller radius also is used for the search of ',/ +10x,'more accurate position of density peak',/ +10x,' Never give the value of the small radius which is ',/ +10x,'smaller than 1.5-2 highest spatial resolution') Rsearch1=INPUT('Enter comoving search radius(Mpc/h) => ') Rsearch2=INPUT('Enter smaller radius(Mpc/h) of final halos=> ') Rminhalo=0.0001 ! min.radius for halos(Mpc/h) Fract = 1. ! fraction of DM particles (1/4,1/2,1) Fract =Min(Max(Fract,0.),1.0) If(iHelp.eq.1)write(*,160) 160 format(10x,'Unbounding particles.',/ + 10x,' =1 for removing unbound particles; =0 or >10', + ' for no unbounding') Toohot =INPUT('Enter rejection velocity limit (V/Vescape)=> ') If(iHelp.eq.1)write(*,170) 170 format(10x,'Halos with close masses, close velocities and with',/ +11x,'distances less than some value (Mpc/h) are considered ',/ +10x,'as the same halo. Only the largest is kept.') dRdubl =INPUT('Distance to check for Velocity duplicates => ') fVdubl =0.15 ! Define duplicates if (v1-v2)/Vrms < if(iHelp.eq.1)Then write (25,*) Ovlim,' Overdensity Threshold for Halos' write (25,*) Rsearch1,' comoving search radius(Mpc/h)' write (25,*) Rsearch2,' smaller radius(Mpc/h) of final halos' write (25,*) Toohot,' rejection velocity limit (V/Vescape)' write (25,*) dRdubl, + ' Distance to check for Velocity duplicates' EndIf C ------------------------------------------ Open files Open( 2,file='Catshort.DAT',Status='UNKNOWN')! short list of halos Open(20,file='Catalog.DAT',Status='UNKNOWN') ! catalog of halos C ------------------------------------------ Read data CALL ReadPnt(N,Fract,xshift,yshift,zshift) ! particles If(iDebug.eq.1)write (*,*)' done Reading Particles' C CALL ReadMax(N,Fract) ! particles c Rsearch1=Rsearch1/AEXPN ! scale to proper distances c Rsearch2=Rsearch2/AEXPN c dRdubl =dRdubl /AEXPN Boxx =Box ! store box size in Mpc WRITE (2,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN,NROW,NGRID,Om0,Oml0,hubble WRITE (20,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN,NROW,NGRID,Om0,Oml0,hubble 100 FORMAT(1X,'Header=>',A45,/ + ' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ + ' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3, + ' Nrow=',I4,' Ngrid=',I4,/' Omega_0=',f6.3, + ' Omega_L=',f6.3,' h=',f6.3) C ------------------------------------------ Constants --------- Omcold= Om0 Omhot = 0. Amassl= Amassl/hsmall ! Scale to real masses Radius= Rsearch1/hsmall ! Comoving Search radius in real Mpc Rminhalo=Rminhalo/hsmall dRdubl= dRdubl/hsmall C Dscale is a factor used to convert particles to mass C = mass of a particle for Omega=1 model C Xscale is a factor to convert PM coordinates into Mpc C = Box / Ngrid = length of a PM cell in Mpc C Vscale is a factor to convert PM pomenta to V (km/sec) C = the Hubble velocity at the distance C of 1 PM cell C the factor AEXPN is needed to go from momenta to pec.velocity If(iDebug.eq.1)write (*,*) ' Vscale=',Vscale, ' X=',Xscale, + ' Box(Mpc)=',Box Radius=Radius/Xscale ! comoving radius in cell units RadSr1 =Rsearch1/hsmall/Xscale ! max Min search radii RadSr2 =Rsearch2/hsmall/Xscale Box =NGRID ! all internal variables are in Grid units Riter =0.0001 ! Error for iteration of radius in cell units Rmax =Radius ! Treat maxima Identical if dR',f5.1,'V_escape were removed', . ' Cell (in grids) for linker list=',f8.3) write (*,*) 'Pairs made' CALL SmallRem(ANlimit,Ncentr) ! Remove small halos write (*,*) 'Small Rem made. Ncentr=',Ncentr C Final touch: if Rsearch1 and Rsearch2 are different C then improve position of halos by shrinking C the search Radius to Rsearch2 in Nit steps Rsearch2 = Rsearch2/3. ! reduce the small radius c RadSr2 =Rsearch2/hsmall/Xscale ! for more accurate search of center IF(ABS(Rsearch1-Rsearch2).gt.0.001*Rsearch1)Then Nit =3 hNit =(Rsearch1/Rsearch2)**(1./Nit) If(iDebug.eq.1)write (*,*) ' Nit =',Nit,' Factor:',hNit Do ii =1,Nit n1 = 0 r1 =0. Do i=1,Ncentr If(Rmc(i).gt.RadSr2)Then n1 = n1 +1 Rmc(i) = max(Rmc(i)/hNit,RadSr2) r1 = r1 +Rmc(i) Endif Enddo write (*,*)'Iterate:', & ' Mean Comoving Radius for halos(Mpc/h)=', & r1/(max(n1,1))*Xscale*hsmall,' n=',n1 CALL Pairs(Box,N,Ncentr,Riter) EndDo Radius=Rsearch1 EndIf dRdubl =dRdubl/Xscale ! scale to cell units CALL CATALOG(Box,Ncentr,dRdubl) ! Remove duplicates If(Ncentr.gt.Ncc)Then write (*,*) ' Too many centers:',Ncentr,' max (Ncc)=',Ncc STOP ' Too many centers. Increase Ncc parameter' Endif Celloptimal = (Rad(Nrad)/FracSearch) Cell = max(Celloptimal,float(NGRID)/float(Nlinker)) Call List(Box,N,Ncp) ! make linker List, Label write (*,20) N,Ncentr,Box*Xscale*hsmall,Ovlim,Cell write (2,20) N,Ncentr,Box*Xscale*hsmall,Ovlim,Cell write (20,20)N,Ncentr,Box*Xscale*hsmall,Ovlim,Cell 20 Format(5x,'Statistics for ',i9,' points. Maxima=',i8,/ . 5x,'Box(Mpc)=',f6.1, . 5x,'Overdensity Limit=',f6.1, . ' Cell for linker list(grids)=',f8.3) t0 =seconds() write (*,*) ' time before Spairs(sec)=',t0 CALL SPairs(Box,N,Ncentr) ! Accumulate statistics t1 =seconds() write (*,*) ' time for Spairs (sec)=',t1-t0 Amassl = AmassOne*9. ! Minimum halo mass in Msun/h C----------------------------- Auxiliary Settings nfalse =0 Do ir=1,Nrad ! all scaled to hubble, not h^(-1) Radsc(ir) =Rad(ir) *Xscale ! scale radii to comoving Mpc MnMass(ir)=Om0*4.188*2.746e+11*hubble**2*Radsc(ir)**3 c mass at the mean density EndDo write (*,12) (Radsc(i)*hsmall*1000.,i=1,Nrad) write (2,12) (Radsc(i)*hsmall*1000.,i=1,Nrad) write (20,12) (Radsc(i)*hsmall*1000.,i=1,Nrad) 12 Format(5x,'Radii of shells in comoving kpc/h:',/20(4x,10f8.2/)) C---------------------------- Process Halos : Unbind, Dublicates CALL CLEAN(Ncentr,Xscale,Vscale,Dscale, + Fract,Radsc,MnMass,Toohot,Rminhalo) write (*,*) ' CLEAN: Ncenters=',Ncentr t2 =seconds() write (*,*) ' time for CLEAN(sec)=',t2-t1 IF(fVdubl .gt.0. .and. dRdubl.gt.0.)Then ! Remove VelDublicates CALL RemVelDublA(Box,Ncentr,dRdubl,fVdubl,Xscale,Vscale) write (*,25) Ncentr,dRdubl*Xscale*hsmall ,fVdubl write (2,25) Ncentr,dRdubl*Xscale*hsmall ,fVdubl write (20,25)Ncentr,dRdubl*Xscale*hsmall ,fVdubl 25 Format(5x,'Number of centers=',i6, & 2x,'Velocity Duplicates: dR(Mpc/h)<',f6.3, & 2x,'dV/Vrms<',f6.3) ENDIF C---------------------------------- Find largest halo c NpMax =0 c inpMax =0 c Do i=1,Ncentr c write (*,*) ' halo=',i,NbcS(1,i),NbcS(2,i),NbcS(3,i) c If(NbcS(2,i).gt. NpMax)Then c inpMax =i c NpMax =NbcS(2,i) c EndIf c EndDo c If(inpMax.eq.0.or.NpMax.eq.0)Write (*,*) ' Errrror' c Do i=1,Ncentr cc If(i.ne.inpMax) iRadius(i) =0 c If(Vrm(i).lt. 100.) iRadius(i) =0 c EndDo c iRadius(inpMax) =Nrad C------------------------------- Print Results iCount = 0 write ( 2,13) write (20,13) 13 Format(6x,'Coordinates Mpc',9x,'Velocity km/s',7x, & 'Mass Radius Vrms(3D) Vcirc Npart(