C------------------------------------------------- c C INCLUDE 'PMparameters.h' INCLUDE 'PMvlists.h' Real INPUT Np_min =INPUT('Enter Minimum number of particles in a halo => ') V_c_min=INPUT('Enter Minimum circular velocity => ') Radius =INPUT('Enter comoving radius(Mpc/h) for spheres => ') Nsph =INPUT('Enter number of random spheres => ') Fract =INPUT('Enter fraction of DM particles (1/4,1/2,1)=> ') Fract =Min(Max(Fract,0.),1.0) Box =INPUT('Enter Comoving Box size(Mpc/h) => ') write (*,*) C ------------------------------------------ Open files Open( 2,file='Vresults.DAT') ! short list of halos Open( 3,file='Catshort.DAT') C ------------------------------------------ Read data CALL RDTAPE write (*,*) ' RDTAPE is done' write (*,*) ' Nnu=Number of Species=',Nspecies WRITE (2,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 --------- hsmall= hubble ! Hubble/100. Box = Box /hsmall ! Scale to real Mpc Radius= Radius/hsmall ! Comoving Search radius in 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 ! mass scale 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 of 1 PM cell C the factor AEXPN is needed to go from momenta to pec.velocity write (*,*) ' Vscale=',Vscale, ' X=',Xscale,' Box(Mpc)=',Box Radius=Radius/Xscale ! comoving radius in cell units RadSr1 = Radius**2 Cell =2.0 ! Cell-size for Linker List in grid units Box =NGRID ! all internal variables are in Grid units imin =INT(-Radius/Cell)! Make recommendations. Code works imax =INT((Box+Radius)/Cell) ! if the limits are not optimal If(imin.lt.Nm)Write(*,*) . ' Wrong low limit:',Nm,' better be',imin If(imax.gt.Nb)Write(*,*) . ' Wrong upper limit:',Nb,' better be',imax write (*,*) . ' Box(cells)=',Box,' Cell=',Cell,' Rad(cells)=',Radius c DtoNscale = Om0*Dscale/Fract ! scale dark matter particles C C CALL ReadCent(Nsph,Ncentr) ! random centers CALL ReadMax(N,Xscale,Np_min,V_c_min) ! read halos Open( 3,file='Catshort.DAT') Ncentr = N Do i=1,N Rmc(i) =Radius Amc(i) =1. Xm(i) =Xp(i) Ym(i) =Yp(i) Zm(i) =Zp(i) EndDo CALL ReadPnt(N,Fract,Vscale) ! read particles Ndm = N c---------------------------------- DM particles CALL Points(Box,N,Ncp) ! Make more points periodically write (*,*) ' Nparticles in periodic set=',Ncp, . '( Maximim =',Np,')' write (2,*) ' Ndm in periodic set=',Ncp, . '( Maximum =',Np,')' Call List(Box,Ncp) ! make linker List, Label write (*,*) ' Nparticles=',N,' Ncenters=',Ncentr CALL SPairs(Box,N,Ncentr,1) ! Accumulate statistics for particles C---------------------------------- Halos CALL ReadMax(N,Xscale,Np_min,V_c_min) ! read halos Nhalos = N CALL Points(Box,N,Ncp) ! Make more points periodically write (*,*) ' N_halos=',Nhalos,' in periodic set=',Ncp, . '( Maximum =',Np,')',' N_dm min=',Np_min, . ' V_circ_min=',V_c_min write (2,*) ' N_halos=',Nhalos,' in periodic set=',Ncp, . '( Maximum =',Np,')' write (2,*) ' N_dm min=',Np_min, . ' V_circ_min=',V_c_min Call List(Box,Ncp) ! make linker List, Label write (*,*) ' Nparticles=',N,' Ncenters=',Ncentr CALL SPairs(Box,N,Ncentr,2) ! Accumulate statistics for halos C------------------------------- Print Results write ( 2,13) 13 Format(7x,'N_dm Contrast N_halos V_dm rmsV_dm', & ' V_halos rmsV_halos') AMmean =Ndm/FLOAT(NGRID**3)*4.188*Radius**3 Do i=1,Ncentr ! Scale all to physical units Xc =Xm(i) *Xscale Yc =Ym(i) *Xscale Zc =Zm(i) *Xscale aM =Nbc(i,1) Contrast = aM/AMmean c If(Contrast.ge.1.)Then VelocityD = sqrt(max(1.e-10, & Wxc(i,1)**2 +Wyc(i,1)**2 +Wzc(i,1)**2))/ & max(Nbc(i,1),1) VelocityH = sqrt(max(1.e-10, & Wxc(i,2)**2 +Wyc(i,2)**2 +Wzc(i,2)**2))/ & max(Nbc(i,2),1) rmsD=sqrt(max(1.e-10,Vrmc(i,1)/max(Nbc(i,1),1) - & VelocityD**2)) rmsH=sqrt(max(1.e-10,Vrmc(i,2)/max(Nbc(i,2),1) - & VelocityD**2)) rmsH2=sqrt(max(1.e-10,Vrmc(i,2)/max(Nbc(i,2),1) - & VelocityH**2)) rmsH3=sqrt(max(1.e-10,Vrmc(i,2)/max(Nbc(i,2),1) )) c write (*,200) i,Nbc(i,1),Contrast ,Nbc(i,2), c & VelocityD,rmsD, VelocityH, rmsH,rmsH2,rmsH3 write (2,200) i,Nbc(i,1),Contrast ,Nbc(i,2), & VelocityD,rmsD, VelocityH, rmsH,rmsH2,rmsH3, & Xc,Yc,Zc 200 format(i5,i7,f8.3,i5,2(f8.2,f8.2,2x),2f8.2,3f7.3) c Endif Enddo Stop End C-------------------------------------------------------------- C Make linker lists of particles in each cell SUBROUTINE List(Box,Ncp) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMvlists.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; C ic = center, js = 1 (dm), =2 (halos) SUBROUTINE SNeib(ic,js) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMvlists.h' Radmax =Radius 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.RadSr1)Then Nbc(ic,js)= Nbc(ic,js)+ 1 Rrc(ic,js)= Rrc(ic,js)+ sqrt(dd) Vrmc(ic,js) =Vrmc(ic,js) + & VXp(jp)**2+VYp(jp)**2+VZp(jp)**2 Wxc(ic,js) =Wxc(ic,js) +VXp(jp) Wyc(ic,js) =Wyc(ic,js) +VYp(jp) Wzc(ic,js) =Wzc(ic,js) +VZp(jp) EndIf jp =Lst(jp) GoTo10 EndIf EndDo EndDo 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) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMvlists.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) 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 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-------------------------------- Update Statistics of pairs C SUBROUTINE SPairs(Box,N,Ncentr,js) C--------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMvlists.h' If(Ncentr.gt.Ncc)Then write (*,*) ' Number of centers is too large:' write (*,*) ' Ncentr =',Ncentr,' Max=Ncc=',Ncc,' STOP' STOP EndIf Do i=1,Ncentr Wxc(i,js) = 0. Wyc(i,js) = 0. Wzc(i,js) = 0. Nbc(i,js) = 0 Rrc(i,js) = 0. Vrmc(i,js)= 0. EndDo Do ic=1,Ncentr Xc =Xm(ic) Yc =Ym(ic) Zc =Zm(ic) Call SNeib(ic,js) EndDo Return End C--------------------------------------------------------- C Select random particles c as initial seeds of centers SUBROUTINE ReadCent(N,Ncentr) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMvlists.h' Nseed =120731 Do i=1,N Rmc(i) =Radius Amc(i) =1. Xm(i) =NGRID*RANDd(Nseed) Ym(i) =NGRID*RANDd(Nseed) Zm(i) =NGRID*RANDd(Nseed) EndDo Ncentr =N Return End C--------------------------------------------------------- C Read list of maxima c Xnew =X-1: coord. now are 0-NGRID SUBROUTINE ReadMax(Ncentr,Xscale,Np_min,V_c_min) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMvlists.h' Character*79 Text Do i=1,16 read(3,'(A)') Text write(*,'(A)') Text enddo Ncentr =0 xxm =1000. xxn =-1000. 10 read (3,*,err=5,end=5) X0, Y0, Z0, VvX,VvY,VvZ, & amh,rh,wvel,Vcirc,ihh If(ihh.ge.Np_min.and.Vcirc.ge.V_c_min)Then Ncentr=Ncentr+1 If(Ncentr.le.Ncc)Then Xp(Ncentr) =X0/Xscale/hubble Yp(Ncentr) =Y0/Xscale/hubble Zp(Ncentr) =Z0/Xscale/hubble VXp(Ncentr)=VvX VYp(Ncentr)=VvY VZp(Ncentr)=VvZ xxm =min(xxm,Xp(Ncentr)) xxn =max(xxn,Xp(Ncentr)) EndIf EndIf If(Ncentr.lt.Ncc)GoTo 10 5 write(*,*) ' Read Halos: X{Min,Max}=',xxm,xxn,' N=',Ncentr CLOSE (3) If(Ncentr.GT.Ncc)Then write (*,*) ' Too many particles: Max=',Ncc Stop EndIf Return End C--------------------------------------------------------- C Read particles SUBROUTINE ReadPnt(N,Fract,Vscale) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMvlists.h' N =0 Iz=1 Iy=1 Ff=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 C Loop over particles DO 10 IROW=1,NROW,Iz READ (21,REC=IROW) RECDAT DO IN=1,NPAGE,Iy N = N +1 Xp(N) =XPAR(IN)-1. Yp(N) =YPAR(IN)-1. Zp(N) =ZPAR(IN)-1. VXp(N)=VX(IN)*Vscale VYp(N)=VY(IN)*Vscale VZp(N)=VZ(IN)*Vscale ENDDO 10 CONTINUE Return End