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 'PMlists.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' 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 *iWeight(1) ! 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),iWeight(1) 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/iWeight(1) 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' CALL CATALOG(Box,Ncentr) ! Remove duplicates write (*,*) 'CATALOG made' 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/2. ! reduce the small radius ! 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)=Rmc(i)/hNit 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 t3 =seconds() write (*,*) ' time after Pairs=',t3 Radius=Rsearch1 EndIf 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) If(Celloptimal.gt.Cell)Then Cell = Celloptimal Call List(Box,Ncp) ! make linker List, Label EndIf 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) CALL SPairs(Box,N,Ncentr) ! Accumulate statistics t3 =seconds() write (*,*) ' timeafter Spairs=',t3 Amassl = AmassOne*5. ! 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) IF(fVdubl .gt.0. .and. dRdubl.gt.0.)Then ! Remove VelDublicates dRdubl =dRdubl/Xscale ! scale to cell units CALL RemVelDubl(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(i,1),NbcS(i,2),NbcS(i,3) c If(NbcS(i,2).gt. NpMax)Then c inpMax =i c NpMax =NbcS(i,2) 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( ') RETURN END C-------------------------------------------------------------- C Make linker lists of particles in each cell SUBROUTINE List(Box,Ncp) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Do i=1,Ncp Lst(i)=-1 EndDo Do k=Nm,Nb Do j=Nm,Nb Do i=Nm,Nb Label(i,j,k)=0 EndDo EndDo EndDo Do jp=1,Ncp i=INT(Xp(jp)/Cell) j=INT(Yp(jp)/Cell) k=INT(Zp(jp)/Cell) i=MIN(MAX(Nm,i),Nb) j=MIN(MAX(Nm,j),Nb) k=MIN(MAX(Nm,k),Nb) Lst(jp) =Label(i,j,k) Label(i,j,k) =jp EndDo Return End C-------------------------------------------------------------- C Find all neibhours for a center C Xc,Yc,Zc - center; a0 -its weight C SUBROUTINE Neib(Xnew,Ynew,Znew,Amnew,Xc,Yc,Zc,Radius) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' C Initiate counters Xnew =0. Ynew =0. Znew =0. Amnew=0. Radius2=Radius**2 c limits for Label i1=INT((Xc-Radius)/Cell) j1=INT((Yc-Radius)/Cell) k1=INT((Zc-Radius)/Cell) i1=MIN(MAX(Nm,i1),Nb) j1=MIN(MAX(Nm,j1),Nb) k1=MIN(MAX(Nm,k1),Nb) i2=INT((Xc+Radius)/Cell) j2=INT((Yc+Radius)/Cell) k2=INT((Zc+Radius)/Cell) i2=MIN(MAX(Nm,i2),Nb) j2=MIN(MAX(Nm,j2),Nb) k2=MIN(MAX(Nm,k2),Nb) C Look for neibhours nn =0 Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) 10 If(jp.ne.0)Then dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2 If(dd.lt.Radius2)Then nn =nn +1 ww =iWeight(jp) Xnew =Xnew +Xp(jp)*ww Ynew =Ynew +Yp(jp)*ww Znew =Znew +Zp(jp)*ww Amnew=Amnew+ww EndIf jp =Lst(jp) GoTo10 EndIf EndDo EndDo EndDo If(Amnew.GT.0.)Then Xnew =Xnew /Amnew Ynew =Ynew /Amnew Znew =Znew /Amnew EndIf Return End C-------------------------------------------------------------- C Find all neibhours for a center C Xc,Yc,Zc - center; C SUBROUTINE SNeib(Xc,Yc,Zc,Vx1,Vy1,Vz1,Nob1,NobS,NobM, & NobG,Rob1,Vxx ) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' DIMENSION Vx1(Nrad),Vy1(Nrad),Vz1(Nrad),Vxx(Nrad),Nob1(Nrad), . Rob1(Nrad),NobM(Nrad),NobG(Nrad),NobS(Nrad) iFlag = 0 Radmax =Rad(Nrad)/FracSearch d0 = Rad(1) 50 dmax = Radmax**2 Ninside = aNvir*Radmax**3 iw1 = 1 iw8 = 8 iw64 = 64 Do i=1,Nrad Nob1(i) =0 NobS(i) =0 NobM(i) =0 NobG(i) =0 Rob1(i) =0. Vx1(i) =0. Vy1(i) =0. Vz1(i) =0. Vxx(i) =0. EndDo c limits for Label i1=INT((Xc-Radmax)/Cell) j1=INT((Yc-Radmax)/Cell) k1=INT((Zc-Radmax)/Cell) i1=MIN(MAX(Nm,i1),Nb) j1=MIN(MAX(Nm,j1),Nb) k1=MIN(MAX(Nm,k1),Nb) i2=INT((Xc+Radmax)/Cell) j2=INT((Yc+Radmax)/Cell) k2=INT((Zc+Radmax)/Cell) i2=MIN(MAX(Nm,i2),Nb) j2=MIN(MAX(Nm,j2),Nb) k2=MIN(MAX(Nm,k2),Nb) nps =0 C Look for neibhours Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) ! Dark matter 10 If(jp.ne.0)Then dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2 If(dd.lt.dmax)Then d1 =sqrt(dd) ir =INT(max(0.,4.*(sqrt(d1/d0)-1.) +1)) +1 ir =1 ir = min(max(1,ir),Nrad) ww =iWeight(jp) iww = INT(ww/iWeight(1)+0.05) Nob1(ir)= Nob1(ir)+ iww c If(iww.eq.iw1)NobS(ir)= NobS(ir)+ 1 c If(iww.eq.iw8)NobM(ir)= NobM(ir)+ 1 c If(iww.eq.iw64)NobG(ir)= NobG(ir)+ 1 c Rob1(ir)= Rob1(ir)+ d1*iww c Vxx(ir) = Vxx(ir) + c + iww*(VXp(jp)**2+VYp(jp)**2+VZp(jp)**2) Vx1(ir) = Vx1(ir) +iww*VXp(jp) Vy1(ir) = Vy1(ir) +iww*VYp(jp) Vz1(ir) = Vz1(ir) +iww*VZp(jp) EndIf jp =Lst(jp) GoTo10 EndIf EndDo EndDo EndDo Do ir =2,Nrad Nob1(ir)= Nob1(ir) + Nob1(ir-1) c NobS(ir)= NobS(ir) + NobS(ir-1) c NobM(ir)= NobM(ir) + NobM(ir-1) c NobG(ir)= NobG(ir) + NobG(ir-1) c Rob1(ir)= Rob1(ir) + Rob1(ir-1) c Vxx(ir) = Vxx(ir) + Vxx(ir-1) Vx1(ir) = Vx1(ir) + Vx1(ir-1) Vy1(ir) = Vy1(ir) + Vy1(ir-1) Vz1(ir) = Vz1(ir) + Vz1(ir-1) EndDo Nobj = Nob1(6) Vx14 = Vx1(6)/Nobj Vy14 = Vy1(6)/Nobj Vz14 = Vz1(6)/Nobj Vmn = Vx14**2+Vy14**2+Vz14**2 Do ir=1,Nrad c Vmn =0. If(Nob1(ir).GT.0)Then Nobj = Nob1(ir) Vx1(ir) = Vx1(ir)/Nobj Vy1(ir) = Vy1(ir)/Nobj Vz1(ir) = Vz1(ir)/Nobj c Vmn = Vx1(ir)**2+Vy1(ir)**2+Vz1(ir)**2 c Vxx(ir) = SQRT(ABS(Vxx(ir)/Nobj -Vmn)) EndIf EndDo If(Nob1(Nrad).gt.Ninside.and.iFlag.eq.0)Then Radmax =Rad(Nrad) Iflag = 1 goto50 EndIf Return End C-------------------------------------------------------------- C Find all neibhours for a center C Xc,Yc,Zc - center; C SUBROUTINE SNeibOld(Xc,Yc,Zc,Vx1,Vy1,Vz1,Nob1,NobS,NobM, & NobG,Rob1,Vxx ) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' DIMENSION Vx1(Nrad),Vy1(Nrad),Vz1(Nrad),Vxx(Nrad),Nob1(Nrad), . Rob1(Nrad),NobM(Nrad),NobG(Nrad),NobS(Nrad) Radmax =Rad2(Nrad) iw1 = 1 iw8 = 8 iw64 = 64 Do i=1,Nrad Nob1(i) =0 NobS(i) =0 NobM(i) =0 NobG(i) =0 Rob1(i) =0. Vx1(i) =0. Vy1(i) =0. Vz1(i) =0. Vxx(i) =0. EndDo c limits for Label i1=INT((Xc-Radmax)/Cell) j1=INT((Yc-Radmax)/Cell) k1=INT((Zc-Radmax)/Cell) i1=MIN(MAX(Nm,i1),Nb) j1=MIN(MAX(Nm,j1),Nb) k1=MIN(MAX(Nm,k1),Nb) i2=INT((Xc+Radmax)/Cell) j2=INT((Yc+Radmax)/Cell) k2=INT((Zc+Radmax)/Cell) i2=MIN(MAX(Nm,i2),Nb) j2=MIN(MAX(Nm,j2),Nb) k2=MIN(MAX(Nm,k2),Nb) nps =0 C Look for neibhours Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) ! Dark matter 10 If(jp.ne.0)Then dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2 c write (*,'(" jp=",i9," d=",7g11.3)') c & jp,dd,Rad2(Nrad),Xc,Yc,Zc Do ir =1,Nrad If(dd.lt.Rad2(ir))Then ww =iWeight(jp) iww = INT(ww/iWeight(1)+0.05) Nob1(ir)= Nob1(ir)+ iww If(iww.eq.iw1)NobS(ir)= NobS(ir)+ 1 If(iww.eq.iw8)NobM(ir)= NobM(ir)+ 1 If(iww.eq.iw64)NobG(ir)= NobG(ir)+ 1 Rob1(ir)= Rob1(ir)+ sqrt(dd)*iww Vxx(ir) = Vxx(ir) + + iww*(VXp(jp)**2+VYp(jp)**2+VZp(jp)**2) Vx1(ir) = Vx1(ir) +iww*VXp(jp) Vy1(ir) = Vy1(ir) +iww*VYp(jp) Vz1(ir) = Vz1(ir) +iww*VZp(jp) c write (*,'(" old ",i9,2i7," ww=",5g11.3)') c & jp,ir,Nob1(ir),sqrt(dd),Rad(ir) EndIf EndDo jp =Lst(jp) GoTo10 EndIf EndDo EndDo EndDo c write(*,*) '------------ Old Nob1',Xc,Yc,Zc c write (*,'(10i11)')Nob1 c write(*,*) ' Old Rob1' c write (*,'(1P,10G11.3)')Rob1 c write(*,*) ' Old Vxx' c write (*,'(1P,10G11.3)')Vxx Nobj =Nob1(6) Vx14 = Vx1(6)/Nobj Vy14 = Vy1(6)/Nobj Vz14 = Vz1(6)/Nobj Vmn = Vx14**2+Vy14**2+Vz14**2 Do ir=1,Nrad c Vmn =0. If(Nob1(ir).GT.0)Then Nobj = Nob1(ir) Vx1(ir) = Vx1(ir)/Nobj Vy1(ir) = Vy1(ir)/Nobj Vz1(ir) = Vz1(ir)/Nobj c Vmn = Vx1(ir)**2+Vy1(ir)**2+Vz1(ir)**2 Vxx(ir) = SQRT(ABS(Vxx(ir)/Nobj -Vmn)) EndIf EndDo Return End C-------------------------------------------------------------- C Add points around the Box to take into account periodical conditions C It replicates points periodically and keeps those which are at distance C less than Radius from the Box. C N = number of points inside the Box C Ncp = total number of points SUBROUTINE Points(Box,N,Ncp,Radius) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Logical Inside B1 =-Radius B2 =Box+Radius B3 =Box-Radius Ncp =N Do i=1,N ! Add dark matter particles x0 =xp(i) y0 =yp(i) z0 =zp(i) Inside = (x0.gt.Radius).and.(x0.lt.B3) Inside =Inside.and.((y0.gt.Radius).and.(y0.lt.B3)) Inside =Inside.and.((z0.gt.Radius).and.(z0.lt.B3)) If(.not.Inside)THEN ux =VXp(i) uy =VYp(i) uz =VZp(i) ww=iWeight(i) Do k1 = -1,1 zr =z0+k1*Box If(zr.GT.B1.AND.zr.LT.B2)Then kk =k1**2 Do j1 = -1,1 yr =y0+j1*Box If(yr.GT.B1.AND.yr.LT.B2)Then jj = kk +j1**2 Do i1 = -1,1 xr =x0+i1*Box ii =jj +i1**2 If((xr.GT.B1.AND.xr.LT.B2).AND.ii.ne.0) . Then Ncp =Ncp +1 If(Ncp.le.Np)Then xp(Ncp) =xr yp(Ncp) =yr zp(Ncp) =zr VXp(Ncp) =ux VYp(Ncp) =uy VZp(Ncp) =uz iWeight(Ncp) =ww EndIf EndIf EndDo EndIf EndDo EndIf EndDo EndIf EndDo If(Ncp.Gt.Np)Then write(*,*)' Too many points:',Ncp,' Maximum is =', Np STOP EndIf Return End C------------------------------------------------------------- C Find positions of 'Ncentr' spheres, which maximize number of particles C inside 'radius Radius'. C Iterations are repeated until 1) mass stops increasing OR C 2) the sphere stops moving: dR is less than Riter C Looks for all particles inside Radius for periodical set SUBROUTINE Pairs(Box,N,Ncentr,Riter) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' INTEGER Lab(Nc) write (*,*) ' Pairs: start',N,Ncentr Do i=1,Ncentr Lab(i) =1 Amc(i) =0 If(Rmc(i).lt.1.e-8)write (*,*) ' Error in Rmc: zero',Rmc(i),i If(Rmc(i).le.0.5*RadSr2) & write (*,*) ' Error in Rmc: too small ',Rmc(i),i If(Rmc(i).gt.RadSr1)write (*,*) ' Error in Rmc: big ',Rmc(i),i c Rmc(i) =Radius EndDo Niter =0 10 Iter =0 Niter =Niter +1 Iter1 =0 Do i=1,Ncentr If(Lab(i).gt.0.and.Rmc(i).ge.0.999*RadSR2)Iter1 =Iter1 +1 EndDo C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (i,Xc,Yc,Zc,a0,in,Radius) C$OMP+PRIVATE (Xnew,Ynew,Znew,Amnew,Dr) C$OMP+SCHEDULE(DYNAMIC,10000) Do i=1,Ncentr Xc =Xm(i) Yc =Ym(i) Zc =Zm(i) a0 =Amc(i) in =Lab(i) Radius =Rmc(i) If(Radius.lt.0.999*RadSR2)in =-Niter If(in.gt.0)Then ! If not converged, keep iterating Call Neib(Xnew,Ynew,Znew,Amnew,Xc,Yc,Zc,Radius) ! new center of mass Dr =MAX(ABS(Xnew-Xc),ABS(Ynew-Yc),ABS(Znew-Zc)) IF(Dr.LT.Riter.or.Amnew.lt.Amc(i)) . Lab(i)=-Niter ! iterations converged c Iter =Iter +1 IF(Xnew.lt.0.)Xnew =Xnew+Box IF(Ynew.lt.0.)Ynew =Ynew+Box IF(Znew.lt.0.)Znew =Znew+Box IF(Xnew.gt.Box)Xnew =Xnew-Box IF(Ynew.gt.Box)Ynew =Ynew-Box IF(Znew.gt.Box)Znew =Znew-Box Xm(i) = Xnew Ym(i) = Ynew Zm(i) = Znew Amc(i) = Amnew c write (*,'(i6,6g12.5)') i, Xm(i),Ym(i),Zm(i),Amc(i) EndIf EndDo Write (*,*) ' Iteration=',Niter,' Centers left=',Iter,Iter1 If(Iter1.GE.1 .and. Niter.LT.100)Goto10 Return End C-------------------------------- Remove dublicates in Velocity C If distance is less than 'dRdubl' and C if difference in velocitiy is less than a fraction C 'fVdubl' of the velocity dispersion --> C dublicates: keep only one (largest) halo C C Condition of proximity in velocity depends on how close are halos C in space. Halos with velocity difference dv =sqrt(vx_1-vx_2)^2 +(y)+ (z)), C at distance d one from the other, with mass ratio 1/1.5 < M1/M2 < 1.5 C and with internal velocity dispersions Vrms_1 and Vrms_2 C are considered too close ('dublicates') C if dv < fVdubl*max(Vrms_1,Vrms_2)*sqrt(6-5(d/dRdubl)^2) SUBROUTINE RemVelDubl(Box,Ncentr,dRdubl,fVdubl,Xscale,Vscale) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Logical Lab(Nc) Real Mhalo,Mhalo2 dR2 = dRdubl**2 fV2 = fVdubl**2 Vol_2 =Rad(2)**3-Rad(1)**3 nrem= 0 Do i=1,Ncentr Lab(i) =.true. EndDo If(iDebug.eq.1)write (*,50) 50 format(5X,'<======== Remove fake halos ========>', & /17x,'difference in velocities dV is small as compared with', & /17x,'internal rms velocity Vrms', & ' and halos are close in space',/ & 4x,'dR(Mpc)',' dV(km/s)',' Vrms(km/s)', & ' N_part_halo_1',' N_part_halo_2', & ' dVx(km/s) dVy(km/s) dVz(km/s)') Do i=1,Ncentr-1 Mhalo =Amc(i) irad =Max(iRadius(i),2) Rh_1 =Rad(irad) Vrh2 =(Vrmc(i,irad)/Vscale)**2 If(Mhalo.eq.0. .or. irad.eq.0)Lab(i) =.false. If(Lab(i))Then Do j=i+1,Ncentr dd =(Xm(j)-Xm(i))**2 + . (Ym(j)-Ym(i))**2 + . (Zm(j)-Zm(i))**2 Mhalo2 =Amc(j) If(dd.lt.dR2.and.Mhalo2.gt.Mhalo/1.5)Then ! close distance pair If(Mhalo2.lt.Mhalo*1.5)Then ! close mass range jrad=Max(iRadius(j),2) Rh_2 =Rad(jrad) Rhh =MAX(Rh_1,Rh_2,1.e-8) dv =(Wxc(j,jrad)-Wxc(i,irad))**2 + . (Wyc(j,jrad)-Wyc(i,irad))**2 + . (Wzc(j,jrad)-Wzc(i,irad))**2 VV =MAX(Vrh2,(Vrmc(j,jrad)/Vscale)**2) ! max Vrms of the two c If(dv.lt.(6.-5.*dd/dR2)*fV2*VV)Then ! close in velocity If(dv.lt.(2.-dd/(dRdubl*Rhh))*fV2*VV)Then ! close in velocity If(Nbc(i,2).gt.Nbc(j,2))Then ! find if density at the distance idens =i ! between halos does not change more than idens_0 =j Else ! five times -- then this is dublicates idens_0 =i idens =j Endif dist =sqrt(dd) Do irr =1,Nrad If(dist.lt.rad(irr))Then idist =irr goto 300 Endif Enddo 300 irr=max(irr,2) rho_cent =(Nbc(idens_0,2)-Nbc(idens_0,1))/Vol_2 rho_dist =(Nbc(idens,irr)-Nbc(idens,irr-1))/ + (Rad(irr)**3-Rad(irr-1)**3) If(iDebug.eq.1)write (*,31) sqrt(dd)*Xscale, + sqrt(dv)*Vscale,sqrt(VV)*Vscale, . Nbc(i,irad),Nbc(j,jrad), . (Wxc(j,jrad)-Wxc(i,irad))*Vscale, + (Wyc(j,jrad)-Wyc(i,irad))*Vscale, + (Wzc(j,jrad)-Wzc(i,irad))*Vscale, + rho_cent,rho_dist,irr,Nbc(idens_0,2), + Nbc(idens,2),idens_0,idens IF(rho_dist.gt.rho_cent/1000.)Then ! dublicates nrem = nrem +1 31 format(f10.4,f9.3,f11.2,2i14,3f10.3,2g11.3,i4,2i5,2i3) 30 format(f10.4,f9.3,f11.2,2i14,3f10.3,2g11.3,i4) c If(Mhalo.gt.Mhalo2)Then ! choose the biggest If(Nbc(i,2).gt.Nbc(j,2))Then ! choose the densiest Lab(j)=.false. c write (*,30) sqrt(dd)*Xscale,sqrt(dv)*Vscale, c + sqrt(VV)*Vscale, c . Nbc(i,irad),Nbc(j,jrad), c . (Wxc(j,jrad)-Wxc(i,irad))*Vscale, c + (Wyc(j,jrad)-Wyc(i,irad))*Vscale, c + (Wzc(j,jrad)-Wzc(i,irad))*Vscale, c + rho_cent,rho_dist,irr c + Xm(j)*Xscale,Ym(j)*Xscale,Zm(j)*Xscale Else Lab(i)=.false. c write (*,30) sqrt(dd)*Xscale,sqrt(dv)*Vscale, c + sqrt(VV)*Vscale, c . Nbc(i,irad),Nbc(j,jrad), c . (Wxc(j,jrad)-Wxc(i,irad))*Vscale, c + (Wyc(j,jrad)-Wyc(i,irad))*Vscale, c + (Wzc(j,jrad)-Wzc(i,irad))*Vscale, c + rho_cent,rho_dist,irr c + Xm(i)*Xscale,Ym(i)*Xscale,Zm(i)*Xscale EndIf EndIf EndIf EndIf EndIf EndDo EndIf EndDo m =0 meanmass =0 Do i=1,Ncentr If(Lab(i))Then m = m +1 Xm(m) =Xm(i) Ym(m) =Ym(i) Zm(m) =Zm(i) Amc(m) =Amc(i) Rmc(m) =Rmc(i) RmaxV(m) =RmaxV(i) Vrm(m) =Vrm(i) iRadius(m) =iRadius(i) Do ir=1,Nrad ! Store Results Wxc(m,ir) =Wxc(i,ir) Wyc(m,ir) =Wyc(i,ir) Wzc(m,ir) =Wzc(i,ir) Nbc(m,ir) =Nbc(i,ir) NbcS(m,ir) =NbcS(i,ir) NbcM(m,ir) =NbcM(i,ir) NbcG(m,ir) =NbcG(i,ir) Rrc(m,ir) =Rrc(i,ir) Vrmc(m,ir) =Vrmc(i,ir) EndDo meanmass =meanmass +Nbc(m,Nrad) EndIf EndDo xmeanmass = FLOAT(meanmass)/FLOAT(m) write(*,10) m,Ncentr,xmeanmass write(20,10) m,Ncentr,xmeanmass write(*,20) dRdubl*Xscale*hubble,fVdubl 10 format(5x,'Vel.Dublicates: New number of objects=',i6,' Old=',i6, . ' mean N_particles=',g11.3) 20 format(' Distance Difference=',f8.4,/ . ' Vrelative/Vrms < ',f8.4) Ncentr =m Return End C-------------------------------- Update Statistics of pairs C SUBROUTINE SPairs(Box,N,Ncentr) C--------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' DIMENSION Vx1(Nrad),Vy1(Nrad),Vz1(Nrad),Vxx(Nrad),Nob1(Nrad), . Rob1(Nrad),NobM(Nrad),NobG(Nrad),NobS(Nrad) Do i=1,Ncentr Do ir=1,Nrad Wxc(i,ir) = 0. Wyc(i,ir) = 0. Wzc(i,ir) = 0. Nbc(i,ir) = 0 NbcS(i,ir) = 0 NbcM(i,ir) = 0 NbcG(i,ir) = 0 Rrc(i,ir) = 0. Vrmc(i,ir)= 0. EndDo EndDo C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (i,ir,Xc,Yc,Zc,Vx1,Vy1,Vz1,Nob1,NobS,NobM) C$OMP+PRIVATE (NobG,Rob1,Vxx ) C$OMP+SCHEDULE (DYNAMIC,2000) Do i=1,Ncentr if(mod(i,10000).eq.0)write (*,*) ' Sneib halo=',i Xc =Xm(i) Yc =Ym(i) Zc =Zm(i) Call SNeib(Xc,Yc,Zc,Vx1,Vy1,Vz1,Nob1,NobS,NobM, & NobG,Rob1,Vxx ) Do ir=1,Nrad ! Store Results Wxc(i,ir) =Vx1(ir) Wyc(i,ir) =Vy1(ir) Wzc(i,ir) =Vz1(ir) Nbc(i,ir) =Nob1(ir) NbcS(i,ir) =NobS(ir) NbcM(i,ir) =NobM(ir) NbcG(i,ir) =NobG(ir) Rrc(i,ir) =Rob1(ir) Vrmc(i,ir)=Vxx(ir) EndDo EndDo Return End C-------------------------------- Remove small halos C C SUBROUTINE SmallRem(ANlimit,Ncentr) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' New =0 Do i=1,Ncentr If(Amc(i).GT.ANlimit*Rmc(i)**3)Then New =New +1 Xm(New) =Xm(i) Ym(New) =Ym(i) Zm(New) =Zm(i) Amc(New)=Amc(i) Rmc(New)=Rmc(i) EndIf EndDo Ncentr =New write (*,*) ' small centers removed. New Ncentr=',Ncentr Return End C-------------------------------- Get Catalog of Halos C Remove close centers, keep only the largest C Emass - small difference in mass SUBROUTINE CATALOG(Box,Ncentr) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Dimension Lab(Ncc) New =0 If(Ncentr.gt.Ncc)Then write (*,*) ' Error: too small Ncc in CATALOG:',Ncc,Ncentr stop EndIf Do i=1,Ncentr Lab(i) =1 EndDo C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (i,j,Xc,Yc,Zc,a0,Rh,dd) C$OMP+SCHEDULE (DYNAMIC,5000) Do i=1,Ncentr If(Rmc(i).lt.1.e-8)Lab(i)=0 Xc =Xm(i) Yc =Ym(i) Zc =Zm(i) a0 =Amc(i) Rh =Rmc(i) Do j=1,Ncentr If(i.ne.j)Then dd =((Xc-Xm(j))**2+(Yc-Ym(j))**2+(Zc-Zm(j))**2) IF(dd.lt.( (Rh+Rmc(j))/2. )**2 )Then If(a0.GT.Amc(j))Then! take largest of maxima Lab(j)=0 Else If(a0.lt.Amc(j))Then Lab(i)=0 Else ! equal masses -> take max(i,j) Lab(min(i,j)) =0 EndIf EndIf EndIf EndIf EndDo EndDo write (*,*) ' CATALOG: done Lab clean',Ncentr c icount =0 c Do i=1,Ncentr c icount =icount +Lab(i) c Enddo c write (*,*) ' Count=',icount c write (*,'(40i2)') (Lab(i),i=1,Ncentr) c STOP C Remove maxima with Label =0 Do i=1,Ncentr If(Lab(i).eq.1)Then New =New +1 Xm(New) =Xm(i) Ym(New) =Ym(i) Zm(New) =Zm(i) Amc(New)=Amc(i) Rmc(New)=Rmc(i) Vrm(New)=Vrm(i) Wxc(New,1)=Wxc(i,1) Wyc(New,1)=Wyc(i,1) Wzc(New,1)=Wzc(i,1) EndIf EndDo Ncentr =New write (*,*) ' dublicates removed. New Ncentr=',Ncentr Return End C--------------------------------------------------------- C Select particles c as initial seeds of centers SUBROUTINE ReadCent(N,Ncentr,Boxx) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Real INPUT,Iweigh_first,Iweigh_last,Iweigh_min Logical Ldiff c Xccc =INPUT(' Center x=') c write (*,*) Xccc c Yccc =INPUT(' Center y=') c write (*,*) Yccc c Zccc =INPUT(' Center z=') c write (*,*) Zccc c Rccc =INPUT(' Center r=') c write (*,*) Rccc c Xccc =Xccc/Boxx*NGRID c Yccc =Yccc/Boxx*NGRID c Zccc =Zccc/Boxx*NGRID c Rccc =Rccc/Boxx*NGRID Ncentr =0 Iweigh_first =iWeight(1) Iweigh_last =iWeight(N) Iweigh_min = Iweigh_first If(iDebug.eq.1) & write (*,*) ' Weights: first species=',Iweigh_first, & ' last=',Iweigh_last,' Weight minimum center=',Iweigh_min n1 = 0 nhalf = 0 n2 = 0 If(RadSr1.lt.1.1*RadSr2.or.ABS(RadSr2).lt.1.e-8)Then ! there is only one radius Ldiff =.false. Radius = RadSr1 R1 = RadSr1 R2 = RadSr1 Rhalf = RadSr1 Else Ldiff =.true. R1 = 1.1*RadSr2 R2 = 0.9*RadSr1 Rhalf =(RadSr1+RadSr2)/2. ialpha =3 RRs = (RadSr1-RadSr2) If(iDebug.eq.1)Then write(*,*) ' Min Radius=',RadSr2,' Max=',RadSr1 write(*,*) ' Power index =',ialpha EndIf Endif If(iHelp.eq.1)Write(*,130) 130 format(/ +10x,'These parameters define how and how many spheres are ',/ +10x,'used to search for density maxima. Tests can be done ',/ +10x,'with few spheres (answers: 100 0). For final search ',/ +10x,'a large fraction of particles should be used. ',/ +10x,'Combinations (np/10 10) or (np/3 5) are recommended',/ +10x,'Only particles of the smallest mass are used for ',/ +10x,'initial centers of shperes') Nline=INPUT(' Enter Number of particles for initial seeds=> ') jneib=INPUT(' Enter Number of neibhours for a seed => ') If(Nline.eq.0)Return fr = N/FLOAT(Nline) xrand =1.00001 Nstep =MAX(INT(fr+1.),1) fr =1./fr Nseed = 121071 write (*,*) ' Every ',Nstep,' particle is used', & ' as initial center',' Fraction=',fr Nplast =0 Do i=1,N ! take centers if weight is less than Iweigh_min If( iWeight(i).le.Iweigh_min)Then Nplast =Nplast +1 Endif Enddo write (*,*) ' Total pool for particles centers=',Nplast C Do i=1,N,Nstep Do i=1,Nplast If(fr.lt.0.999)xrand =RANDd(Nseed) IF(xrand.lt.fr)Then c IF(ABS(Xp(i)-Xccc).lt.Rccc.and. c . ABS(Yp(i)-Yccc).lt.Rccc.and. c . ABS(Zp(i)-Zccc).lt.Rccc)THEN Ncentr =Ncentr+1 IF(Ldiff)Then ! different radii Rrrr =RadSr2+ & RRs*(MAX(RANDd(Nseed),1.e-8))**ialpha Rmc(Ncentr) = Rrrr If(Rrrr.lt.R1)n1 = n1+1 If(Rrrr.lt.R2.and.Rrrr.ge.R1)nhalf = nhalf+1 If(Rrrr.gt.R2)n2 = n2+1 Else Rmc(Ncentr) =Radius Endif Amc(Ncentr) =1. Xm(Ncentr) =Xp(i) Ym(Ncentr) =Yp(i) Zm(Ncentr) =Zp(i) j=Ncentr c EndIf ! endif for small box test EndIf EndDo write (*,*) ' Ncenters (random fraction)=',Ncentr If(jneib.eq.0)goto 150 ! No more points Imax =MIN(INT(NGRID/Cell),Nb) write (*,*) ' Size of Grid centers=',Imax,Cell,NGRID Do k=0,Imax Do j=0,Imax Do i=0,Imax jp =Label(i,j,k) Ncount =0 jcount =0 10 jcount= jcount+1 If(jp.ne.0.and.jcount.lt.jneib)Then jp=Lst(jp) ! find if there are enough neib goto 10 endif jcount =0 If(jp.ne.0.and.iWeight(jp).le.Iweigh_min)Then Ncount =Ncount +1 c IF(ABS(Xp(jp)-Xccc).lt.Rccc.and. c . ABS(Yp(jp)-Yccc).lt.Rccc.and. c . ABS(Zp(jp)-Zccc).lt.Rccc)THEN Ncentr=Ncentr+1 If(Ncentr.le.Nc)Then IF(Ldiff)Then ! different radii Rrrr =RadSr2+ & RRs*(MAX(RANDd(Nseed),1.e-8))**ialpha Rmc(Ncentr) = Rrrr If(Rrrr.lt.R1)n1 = n1+1 If(Rrrr.lt.R2.and.Rrrr.ge.R1)nhalf = nhalf+1 If(Rrrr.gt.R2)n2 = n2+1 Else Rmc(Ncentr) =Radius Endif Amc(Ncentr) =1. Xm(Ncentr) =Xp(jp) Ym(Ncentr) =Yp(jp) Zm(Ncentr) =Zp(jp) EndIf c EndIf ! endif for small-box test c write (*,*) ' Count=',Ncount,i,j,k,Ncentr If(Ncount.lt.3) Goto 10 EndIf EndDo EndDo EndDo 150 write(*,*) ' Ncentr (total)=',Ncentr,' Max=',Nc If(Ncentr.gt.Nc)Then write (*,*) ' Too many centers' STOP ' Too many centers' Endif If(Ncentr.gt.0)Then write (*,200) RadSR2,R1,n1,R1,R2,nhalf,R2,RadSR1,n2 write (2,200) RadSR2,R1,n1,R1,R2,nhalf,R2,RadSR1,n2 write (20,200) RadSR2,R1,n1,R1,R2,nhalf,R2,RadSR1,n2 200 format(' Number of centers with radii between ',3x,2f7.4, & ' was ',i8, & /18x,' with radii between ',3x,2f7.4,' was ',i8, & /18x,' with radii between ',3x,2f7.4,' was ',i8) Endif Return End C--------------------------------------------------------- C Read list of maxima c Xnew =X-1: coord. now are 0-NGRID SUBROUTINE ReadMax(N,Xscale) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' read(3) + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN, + NROWC,NGRIDC,NREC0,Om0,Oml0,hubble, + Ocurv write(*,200) + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN, + NROWC,NGRIDC,NREC0,Om0,Oml0,hubble, + Ocurv 200 FORMAT( + 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=',I9,/ + 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/ + 1x,' Omega_curvature=',F7.3) N =0 xxm =1000. xxn =-1000. 10 read (3,err=5,end=5) X0, Y0, Z0, VvX,VvY,VvZ N=N+1 If(N.le.Np)Then Xp(N) =X0-1. Yp(N) =Y0-1. Zp(N) =Z0-1. VXp(N)=VvX VYp(N)=VvY VZp(N)=VvZ xxm =min(xxm,Xp(N)) xxn =max(xxn,Xp(N)) EndIf GoTo 10 5 write(*,*) ' Read Particles: X{Min,Max}=',xxm,xxn,' N=',N If(N.GT.Np)Then write (*,*) ' Too many particles: Max=',Np Stop EndIf Return End C--------------------------------------------------------- C Read particles SUBROUTINE ReadPnt(N,Fract,xshift,yshift,zshift) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' N =0 Iz=1 Iy=1 Ff=1. c xshift = 100. c yshift = 100. c zshift = 0. xshift = 0. yshift = 0. zshift = 0. xmin =1.e10 xmax =-1.e10 ymin =1.e10 ymax =-1.e10 zmin =1.e10 zmax =-1.e10 Nsmall =0 Lsmall =lspecies(1) If(Fract.le.0.666)Then Iz=2 Ff=0.5 If(Fract.le.0.333)Then Iy=2 Ff=0.25 Endif EndIf Fract =Ff 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 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) EndIf c write (*,*) ' Pages=',Npages,' Species=',Nspecies c write (*,*) ' N_in_last=',N_in_last DO IROW=1, Npages,Iz ! Loop over particle pages c DO IROW=32,66 In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last iL = NPAGE*(IROW-1) C Loop over particles READ (21,REC=IROW) RECDAT DO IN=1,In_page,Iy N = N +1 Xp(N) =XPAR(IN)-1.+xshift Yp(N) =YPAR(IN)-1.+yshift Zp(N) =ZPAR(IN)-1.+zshift If(Xp(N).gt.NGRID)Xp(N) =Xp(N)-NGRID If(Yp(N).gt.NGRID)Yp(N) =Yp(N)-NGRID If(Zp(N).gt.NGRID)Zp(N) =Zp(N)-NGRID If(N.le.Lsmall)Then Nsmall =Nsmall + 1 xmin =min(Xp(N),xmin) xmax =max(Xp(N),xmax) ymin =min(Yp(N),ymin) ymax =max(Yp(N),ymax) zmin =min(Zp(N),zmin) zmax =max(Zp(N),zmax) EndIf VXp(N)=VX(IN) VYp(N)=VY(IN) VZp(N)=VZ(IN) Ipart =min(IN+iL,Nmaxpart) ! current particle number iWeight(N) =iWeight(Ipart) ! particles weight ENDDO ENDDO write(*,*) ' Small Particles:',Nsmall write(*,*) ' x:',xmin,xmax write(*,*) ' y:',ymin,ymax write(*,*) ' z:',zmin,zmax Return End C--------------------------------------------------- C Remove Unbound particles, find radius and mass C No removal of unbound particles If: C Toohot =< 0 or Toohot =>5 C Vscale =100.*hsmall*Xscale/AEXPN C Xscale=Box/Ngrid C Dscale=2.75e+11*hsmall**2*(Box/NROW)**3 SUBROUTINE CLEAN(Ncentr,Xscale,Vscale,Dscale, + Fract,Radsc,MnMass,Toohot,Rminhalo) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' dimension Mass(Nrad),MnMass(Nrad),Overdens(Nrad),Radsc(Nrad) Real Mass,MnMass,Mhalo C.... Open_MP C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (i,ir,iter,Xc,Yc,Zc,Summass,Mass,Overdens) C$OMP+PRIVATE (irad,Mhalo,Rhalo,Vmaxrot,Rmaxrot,Vcirc,Vcc) C$OMP+PRIVATE (NobG,Rob1,Vxx ) C$OMP+SCHEDULE (DYNAMIC,2000) Do i=1,Ncentr if(mod(i,10000).eq.0)write(*,*)' Clean halo=',i,Ncentr Xc =Xm(i) *Xscale ! Scale all to physical units: comoving Yc =Ym(i) *Xscale Zc =Zm(i) *Xscale iflag =0 ifirst=0 Do ir =1,Nrad Vrmc(i,ir)=Vrmc(i,ir) *Vscale ! Scale to km/s Summass = Om0*Nbc(i,ir)/Fract If(Nbc(i,ir).gt.5.and.iflag.eq.0)Then ifirst = ir ! find first significant bin iflag = 1 EndIf Mass(ir) = Dscale*Summass Overdens(ir)= Mass(ir)/MnMass(ir) EndDo C find size, mass, rad ... for a halo Call Decide(Mass,Radsc,Overdens,Nrad,irad,Mhalo,Rhalo, + Vmaxrot,Rmaxrot,ifirst ) Vcirc=sqrt(Vmaxrot/AEXPN)*6.58e-5 ! =10^-5*sqrt(G 2e33/3.08e24) If(Toohot.gt.0. .and. Toohot.lt.5.)Then Do iter =0,4 ! Remove High Energy particles from halos If(irad .ge.1 .and. Rhalo.gt.Rminhalo)Then Vcc =Vcirc/Vscale ! back to dimensionless units Vmaxrot=Vcc*Toohot*1.1**(4-iter) ! Increase the Rotatioanl Vel. Call RemEng(i,Vmaxrot,Rmaxrot/Xscale) ! Remove Unbound Particles iflag =0 ifirst=0 Do ir =1,Nrad If(Nbc(i,ir).gt.5.and.iflag.eq.0)Then ifirst = ir ! find first significant bin iflag = 1 EndIf Vrmc(i,ir)=Vrmc(i,ir) *Vscale Summass =Om0*Nbc(i,ir)/Fract Mass(ir) = Dscale*Summass Overdens(ir)=Mass(ir)/MnMass(ir) EndDo Call Decide(Mass,Radsc,Overdens,Nrad,irad,Mhalo,Rhalo, + Vmaxrot,Rmaxrot,ifirst) Vcirc=sqrt(Vmaxrot/AEXPN)*6.58e-5 !=10^-5*sqrt(2e33/3.08e24) c If(Vcirc.gt.100.)write (*,'(5x,f7.2,i3,2g11.3)') c & Vcirc,irad,Mhalo,Rhalo EndIf EndDo ! end iter energy EndIf Amc(i) =Mhalo ! store results Rmc(i) = Rhalo RmaxV(i) = Rmaxrot Vrm(i) =Vcirc iRadius(i) =irad c write(*,'(i6,1P,5g11.3)')i,Rmaxrot,Vcirc,Rhalo,Mhalo Enddo ! end i (Ncentr) Return End C------------------------------------------------------------- C Find halo radius and mass SUBROUTINE Decide(Mass,Rad,Overdens,Nrad,irad,Mhalo,Rhalo, + Vmaxrot,Rmaxrot,ifirst ) C-------------------------------------------------------------- COMMON /OVER/ Ovlim,Ovcnt Real Rad(Nrad),Overdens(Nrad),Mass(Nrad),Mhalo DATA slope/0.999/ Vmaxrot =0 Rmaxrot =0. Mhalo=0. Rhalo=0. If(ifirst.eq.0)Then irad =0 ! too few particles RETURN EndIf If(MAX(Overdens(ifirst),Overdens(2)).LT.Ovcnt)Then irad =0 ! max overdensity is less then Ovcnt RETURN EndIf iov =0 Do ir=Nrad,ifirst,-1 If(Overdens(ir).gt.Ovlim .and. iov.eq.0)Then iov =ir EndIf EndDo If(iov.eq.0)Then irad =0 ! max overdensity is less then Ovlim RETURN EndIf iovv =0 Do ir=ifirst,iov irleft =max(ir-1,1) irr1 =min(ir+1,Nrad) irr2 =min(ir+2,Nrad) Overright =(2.*Overdens(irr1)+Overdens(irr2))/3. Overleft =(2.*Overdens(ir)+Overdens(irleft))/3. If(Overright.gt.0.97*Overleft)Then !stop:rising profile If(ir.eq.1)Then ! softer conditions for the first bin If(Overright.gt.1.2*Overleft)Then iovv =ir goto 20 EndIf else iovv =ir goto 20 EndIf EndIf EndDo 20 If(iovv.ne.0)iov =iovv irad =iov If(irad.lt.ifirst)write (*,*)' err: irad=',irad,ifirst 30 irleft =max(irad-2,1) ! check if overdens is steep irright =min(irad+1,Nrad) gradlf =Overdens(irleft) * (Rad(irleft))**slope gradrt =Overdens(irright)* (Rad(irright))**slope If(gradlf.lt.gradrt)Then ! not a steep gradient irad =irad -1 ! reduce radius iovv =irad If(irad.lt.ifirst)Then ! zero radius -> quit irad =0 RETURN Else goto 30 ! iterate EndIf EndIf If(irad.lt.ifirst)write (*,*)' err2: irad=',irad,ifirst If(irad.lt.Nrad)Then ! interpolate to Overlim If(Ovlim.gt.Overdens(irad+1) .and. . Overdens(irad).gt.Ovlim)Then Mhalo=(Mass(irad)*(Ovlim-Overdens(irad+1))+ . Mass(irad+1)*(Overdens(irad)-Ovlim)) / . (Overdens(irad)-Overdens(irad+1)) Rhalo=(Rad(irad)*(Ovlim-Overdens(irad+1))+ . Rad(irad+1)*(Overdens(irad)-Ovlim)) / . (Overdens(irad)-Overdens(irad+1)) Else ! Overdens is > Overlim Mhalo =Mass(irad) Rhalo =Rad(irad) EndIf Else ! Overdens is > Overlim Mhalo =Mass(irad) Rhalo =Rad(irad) EndIf imax = 0 Do ir =1,irad ! find max rotational velocity V2 = Mass(ir)/Rad(ir) If (V2.gt.Vmaxrot)Then Vmaxrot =V2 Rmaxrot =Rad(ir) imax = ir EndIf EndDo ! improve Rmaxrot value by finding ! radius containing F(1)/F(2.15)=0.415 fraction ! of mass at Vmax. Then Rmax=2.15*R imax = max(imax,1) aMrs = 0.415*Mass(imax) c write (*,*) ' == imax=',imax,' mass=',aMrs c write (*,'(1P,10G11.3)') Mass Do ir=1,irad If(Mass(ir).gt.aMrs)Then if(ir.eq.1)Then aM0 =0. r0 =0. Else aM0 =Mass(ir-1) r0 =Rad(ir-1) EndIf Rmaxrot2 =2.15*(r0+(Rad(ir)-r0)*(aMrs-aM0)/ & max(Mass(ir)-aM0,1.e-3)) Rmaxrot =min(Rmaxrot,Rmaxrot2,Rhalo/2.) return EndIf EndDo 50 write (*,'(4i4,1P,10G11.3)') ir,irad,ifirst,iovv, & Mass(ir),Ovlim, Rmaxrot, & Rmaxrot2,Rad(ir),r0 write (*,'(1P,10g11.3)')Overdens Return End C------------------------------------------------ C Remove unbound particles for i-th halo C irad = radius of the halo C Vmaxrot = max rotational Velocity C Rmaxrot = radius at max rot.velocity SUBROUTINE RemEng(i,Vmaxrot,Rmaxrot) C----------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' DIMENSION VXxc(Nrad), VYyc(Nrad), VZzc(Nrad), & Vescape2(Nrad) DIMENSION Vx1(Nrad),Vy1(Nrad),Vz1(Nrad),Vxx(Nrad),Nob1(Nrad), . Rob1(Nrad),NobM(Nrad),NobG(Nrad),NobS(Nrad) Xc =Xm(i) ! store center of mass and velocity Yc =Ym(i) Zc =Zm(i) Do ir =1,Nrad ! store mean velocity VXxc(ir) =Wxc(i,ir) VYyc(ir) =Wyc(i,ir) VZzc(ir) =Wzc(i,ir) ! escape velocity for NFW profile Vescape2(ir) = (2.15*Vmaxrot)**2 * & log(1.0+2.*Rad(ir)/Rmaxrot)/(Rad(ir)/Rmaxrot) Wxc(i,ir) = 0. Wyc(i,ir) = 0. Wzc(i,ir) = 0. Nbc(i,ir) = 0 NbcS(i,ir) = 0 NbcM(i,ir) = 0 NbcG(i,ir) = 0 Rrc(i,ir) = 0. Vrmc(i,ir)= 0. EndDo Call SrejNeib(VXxc,VYyc,VZzc,Xc,Yc,Zc,Vescape2, & Vx1,Vy1,Vz1,Nob1,NobS,NobM,NobG,Rob1,Vxx) Do ir=1,Nrad ! Store Results Wxc(i,ir) =Vx1(ir) Wyc(i,ir) =Vy1(ir) Wzc(i,ir) =Vz1(ir) Nbc(i,ir) =Nob1(ir) NbcS(i,ir) =NobS(ir) NbcM(i,ir) =NobM(ir) NbcG(i,ir) =NobG(ir) Rrc(i,ir) =Rob1(ir) Vrmc(i,ir)=Vxx(ir) EndDo Return End C-------------------------------------------------------------- C Accumilate statistics only C if relative velocity is not large C Xc,Yc,Zc - center; C SUBROUTINE SrejNeib(VXxc,VYyc,VZzc,Xc,Yc,Zc,Vescape2, & Vx1,Vy1,Vz1,Nob1,NobS,NobM,NobG,Rob1,Vxx) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' DIMENSION VXxc(Nrad), VYyc(Nrad), VZzc(Nrad),Vescape2(Nrad) DIMENSION Vx1(Nrad),Vy1(Nrad),Vz1(Nrad),Vxx(Nrad),Nob1(Nrad), . Rob1(Nrad),NobM(Nrad),NobG(Nrad),NobS(Nrad) iFlag = 0 Radmax =Rad(Nrad)/FracSearch d0 = Rad(1) 50 dmax = Radmax**2 Ninside = aNvir*Radmax**3 iw1 = 1 iw8 = 8 iw64 = 64 Do i=1,Nrad Nob1(i) =0 NobS(i) =0 NobM(i) =0 NobG(i) =0 Rob1(i) =0. Vx1(i) =0. Vy1(i) =0. Vz1(i) =0. Vxx(i) =0. EndDo c limits for Label i1=INT((Xc-Radmax)/Cell) j1=INT((Yc-Radmax)/Cell) k1=INT((Zc-Radmax)/Cell) i1=MIN(MAX(Nm,i1),Nb) j1=MIN(MAX(Nm,j1),Nb) k1=MIN(MAX(Nm,k1),Nb) i2=INT((Xc+Radmax)/Cell) j2=INT((Yc+Radmax)/Cell) k2=INT((Zc+Radmax)/Cell) i2=MIN(MAX(Nm,i2),Nb) j2=MIN(MAX(Nm,j2),Nb) k2=MIN(MAX(Nm,k2),Nb) C Look for neibhours Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) ! Dark matter 10 If(jp.ne.0)Then dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2 If(dd.lt.dmax)Then d1 =sqrt(dd) iov =INT( max(0.,4.*(sqrt(d1/d0)-1.) +1)) +1 iov = min(max(1,iov),Nrad) dv = (VXxc(iov)-VXp(jp))**2+ + (VYyc(iov)-VYp(jp))**2+ + (VZzc(iov)-VZp(jp))**2 If(dv .lt. Vescape2(iov))Then ww =iWeight(jp) iww =INT(ww/iWeight(1)+0.005) Nob1(iov)= Nob1(iov)+ iww If(iww.eq.iw1)NobS(iov)= NobS(iov)+ 1 If(iww.eq.iw8)NobM(iov)= NobM(iov)+ 1 If(iww.eq.iw64)NobG(iov)= NobG(iov)+ 1 Rob1(iov)= Rob1(iov)+ d1*iww Vxx(iov) = Vxx(iov) + + iww*(VXp(jp)**2+VYp(jp)**2+VZp(jp)**2) Vx1(iov) = Vx1(iov) +iww*VXp(jp) Vy1(iov) = Vy1(iov) +iww*VYp(jp) Vz1(iov) = Vz1(iov) +iww*VZp(jp) EndIf EndIf jp =Lst(jp) GoTo10 EndIf EndDo EndDo EndDo Do ir =2,Nrad Nob1(ir)= Nob1(ir) + Nob1(ir-1) NobS(ir)= NobS(ir) + NobS(ir-1) NobM(ir)= NobM(ir) + NobM(ir-1) NobG(ir)= NobG(ir) + NobG(ir-1) Rob1(ir)= Rob1(ir) + Rob1(ir-1) Vxx(ir) = Vxx(ir) + Vxx(ir-1) Vx1(ir) = Vx1(ir) + Vx1(ir-1) Vy1(ir) = Vy1(ir) + Vy1(ir-1) Vz1(ir) = Vz1(ir) + Vz1(ir-1) EndDo Do ir=1,Nrad Vmn =0. If(Nob1(ir).GT.0)Then Nobj =Nob1(ir) Vx1(ir) = Vx1(ir)/Nobj Vy1(ir) = Vy1(ir)/Nobj Vz1(ir) = Vz1(ir)/Nobj Vmn = Vx1(ir)**2+Vy1(ir)**2+Vz1(ir)**2 Vxx(ir) = SQRT(ABS(Vxx(ir)/Nobj -Vmn)) EndIf EndDo If(Nob1(Nrad).gt.Ninside.and.iFlag.eq.0)Then Radmax =Rad(Nrad) Iflag = 1 goto50 EndIf Return End C-------------------------------------------------------------- C Accumilate statistics only C if relative velocity is not large C Xc,Yc,Zc - center; C SUBROUTINE SrejNeibOld(VXxc,VYyc,VZzc,Xc,Yc,Zc,Vescape2, & Vx1,Vy1,Vz1,Nob1,NobS,NobM,NobG,Rob1,Vxx) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' DIMENSION VXxc(Nrad), VYyc(Nrad), VZzc(Nrad),Vescape2(Nrad) DIMENSION Vx1(Nrad),Vy1(Nrad),Vz1(Nrad),Vxx(Nrad),Nob1(Nrad), . Rob1(Nrad),NobM(Nrad),NobG(Nrad),NobS(Nrad) Radmax =Rad(Nrad) alph =Rad(1)**2 iw1 = 1 iw8 = 8 iw64 = 64 Do i=1,Nrad Nob1(i) =0 NobS(i) =0 NobM(i) =0 NobG(i) =0 Rob1(i) =0. Vx1(i) =0. Vy1(i) =0. Vz1(i) =0. Vxx(i) =0. EndDo c limits for Label i1=INT((Xc-Radmax)/Cell) j1=INT((Yc-Radmax)/Cell) k1=INT((Zc-Radmax)/Cell) i1=MIN(MAX(Nm,i1),Nb) j1=MIN(MAX(Nm,j1),Nb) k1=MIN(MAX(Nm,k1),Nb) i2=INT((Xc+Radmax)/Cell) j2=INT((Yc+Radmax)/Cell) k2=INT((Zc+Radmax)/Cell) i2=MIN(MAX(Nm,i2),Nb) j2=MIN(MAX(Nm,j2),Nb) k2=MIN(MAX(Nm,k2),Nb) C Look for neibhours Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) ! Dark matter 10 If(jp.ne.0)Then dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2 If(dd.lt.Rad2(Nrad))Then c iov =max(1,INT(4.*sqrt(sqrt(dd/alph)) ) -2) iov =0 Do ir =1,Nrad ! find bin If(dd.le.Rad2(ir))Then iov =ir goto 30 EndIf EndDo 30 If(iov.eq. 0)Then ! error write (*,*) ' Error SrejNeib: iov =',iov,' d=',sqrt(dd) Stop EndIf dv = (VXxc(iov)-VXp(jp))**2+ + (VYyc(iov)-VYp(jp))**2+ + (VZzc(iov)-VZp(jp))**2 If(dv .lt. Vescape2(iov))Then ww =iWeight(jp) iww =INT(ww/iWeight(1)+0.005) Do ir =iov,Nrad ! Nrad Nob1(ir)= Nob1(ir)+ iww If(iww.eq.iw1)NobS(ir)= NobS(ir)+ 1 If(iww.eq.iw8)NobM(ir)= NobM(ir)+ 1 If(iww.eq.iw64)NobG(ir)= NobG(ir)+ 1 Rob1(ir)= Rob1(ir)+ sqrt(dd)*iww c Vrr = - (VXp(jp)*(Xc-Xp(jp))+ ! radial velocity c & VYp(jp)*(Yc-Yp(jp))+ c & VZp(jp)*(Zc-Zp(jp)) )/sqrt(dd) c Vxx(ir) = Vxx(ir) + iww*Vrr Vxx(ir) = Vxx(ir) + + iww*(VXp(jp)**2+VYp(jp)**2+VZp(jp)**2) Vx1(ir) = Vx1(ir) +iww*VXp(jp) Vy1(ir) = Vy1(ir) +iww*VYp(jp) Vz1(ir) = Vz1(ir) +iww*VZp(jp) EndDo EndIf EndIf jp =Lst(jp) GoTo10 EndIf EndDo EndDo EndDo Do ir=1,Nrad Vmn =0. If(Nob1(ir).GT.0)Then Nobj =Nob1(ir) Vx1(ir) = Vx1(ir)/Nobj Vy1(ir) = Vy1(ir)/Nobj Vz1(ir) = Vz1(ir)/Nobj Vmn = Vx1(ir)**2+Vy1(ir)**2+Vz1(ir)**2 Vxx(ir) = SQRT(ABS(Vxx(ir)/Nobj -Vmn)) EndIf EndDo Return End C-------------------------------------------------- C Set Weights of particles for C fast access SUBROUTINE SetWeights C-------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Wpr =NGRID**3/FLOAT(NROW**3) If(Nspecies.eq.0)Then ! old constant weights N_particles =lspecies(1) lspecies(1) = NROW**3 N_particles =lspecies(1) Nspecies = 1 wspecies(1) = 1. Ww = 1. Do i=1,N_particles iWeight(i) =Ww EndDo Else N_particles =min(lspecies(Nspecies),Nmaxpart) jstart =1 If(iDebug.eq.1) + write (*,*) ' N_Particles=',N_particles If(N_particles.le.0)Then write (*,*) ' Wrong number of particles:',N_particles write (*,*) ' Nspecies=',Nspecies,' N=',lspecies STop endif Do j=1,Nspecies jend =min(lspecies(j),Nmaxpart) Do k=jstart ,jend iWeight(k) =wspecies(j)/Wpr EndDo jstart =jend EndDo EndIf If(iDebug.eq.1)Then write (*,*) ' Set Weights for ',N_particles,' particles:', & ' w_first=' ,iWeight(1),' w_last=', iWeight(N_particles) write (*,*) ' Nspecies=',Nspecies,' Weights=',wspecies(1), & wspecies(2), wspecies(3) EndIf RETURN END C----------------------------------------------------------- C Open file with answers for rerunning C Print help and description SUBROUTINE Help INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Open(25,file='InHalos.dat',Status='UNKNOWN') Write (*,100) 100 format(//,55('-'),/ +1x,'You run BDM - Bound Density Maxima - code'/ +5x,'Code reads data from two files: PMcrd.DAT', + ' and PMcrs0.DAT',/ +5x,'It writes output to two files :', + ' Catshort.DAT and Catalog.DAT',/ +5x,'Code will ask you few questions to set parameters', + ' of the search',/ +5x,'Your answers will be written into file InHalos.dat') RETURN END c ------------------------ real function seconds () c ------------------------ c c purpose: returns elapsed time in seconds c uses : xlf90 utility function timef() c (in this case - to be compiled with xlf90 or with -lxlf90) c or dtime (Exemplar) or etime (Power Challenge) c real*8 timef real*8 dummy real tarray(2) dummy = 0.0 c.... for IBM SP dummy = timef () seconds = dummy / 1000 c.... for exemplar c dummy = dtime(tarray) c seconds = dtime(tarray) c.... for hitachi c dummy = dtime(tarray) c seconds = second() c.... for power challenge c dummy = etime(tarray) c seconds = etime(tarray) return end