!-------------------------------------------------------------------- ! Read Catshort.DAT ! Construct halos velocity/Mass function !-------------------------------------------------------------------- Module SetUp PARAMETER (Nmax =520000 ) PARAMETER (Nshells = 500 ) PARAMETER (hubble = 73.d0 ) PARAMETER (Box = 160.d0 ) PARAMETER ( aMasslimit = 1.d9 ) PARAMETER ( Vclimit = 50.d0 ) ! Read and use halos above this limit PARAMETER ( dR = 2.0 ) ! buffer width for periodic boundaries PARAMETER ( pi = 3.14159265) Real*4, dimension(Nmax) :: xc,yc,zc, & vxc,vyc,vzc, & vc=0., vrms=0., aMass=0., & Rvir=0.,Nbins=0,Conct=0. Integer, dimension(Nmax) :: Liso =1, Ltrue =0 Real, ALLOCATABLE, DIMENSION(:) :: Xmax,Ymax,Zmax Real*4, dimension(-Nshells:Nshells) :: rdsh,aMsh,Vdsh,Vcsh Integer, dimension(-Nshells:Nshells) :: nsh=0 Integer :: Nh,Nhext,Nclusters,Nmaxima Character Line*79,LineZ*45 End Module SetUp !-------------------------------------------------------------------- ! Find isolated halos in a given mass range in Catalog.DAT ! Construct their averaged density and vel profile !-------------------------------------------------------------------- PROGRAM HaloFunction use SetUP CHARACTER Name*80 INTEGER sName ! Open(1,file='Catshort.00195.DAT') ! Open(1,file='CatshortA_00486.DAT') write (*,'(A,$)')' Enter 3 characters for file name=>' read (*,*) sName write(Name,'(a,i3.3,a)')'results.',sName,'.dat' Open(2,file=Name) write(Name,'(a,i3.3,a)')'Catshort.',sName,'.DAT' Open(1,file=Name) ! Write headers write (2,'(T10,"Hubble=",f8.2,T25,"Box=",f8.2,"Mpc/h", & T35,"Vc limits=",3f8.2,"km/s")') & hubble,Box,Vc1,Vc2,Vclimit write (2,'(T10,"Isolation radius=",f8.2)') & RadLmax Lines = 27 ! skip header of catalog file Do i=1,Lines read(1,'(A)')Line if(i.le.40)write(2,*)Line if(i.le.40)write(*,*)Line EndDo Nread =0 Nh = 1 ! read halo catalog 10 read(1,*,end=100,err=100)xc(Nh),yc(Nh),zc(Nh), & vxc(Nh),vyc(Nh),vzc(Nh),aMass(Nh), & Rvir(Nh),Vrms(Nh),vc(Nh),nn,Rmax,conct(Nh) !,Nbins(Nh) Nread =Nread +1 if(mod(Nh,20000)==0) write(*,'(2i9,3f9.3,i4)') Nread,Nh,vc(Nh),conct(Nh) !,Nbins(Nh) if(vc(Nh).lt.Vclimit)goto 10 ! do not take halo if too small if(Nh.ge.Nmax)goto 100 !write(*,*) Nh,Nbins(Nh),vc(Nh) Nh =Nh +1 goto 10 100 close(1) Nh =Nh -1 write (2,'(" Number of halos read=",T30,i8)')Nh write (*,*) '<====== Number of halos read=',Nh write (*,30) (i,xc(i),yc(i),zc(i),aMass(i),vc(i),i=1,1) write (*,30) (i,xc(i),yc(i),zc(i),aMass(i),vc(i),i=Nh,Nh) 30 format(' halo=',i8,' x=',3g11.3,' M=',g11.3,' Vc=',g11.3) CALL ExpandSample CALL Isolated3D Ltrue = 1 ! make all halos real CALL GetDistribution ! CALL GetMaxima ! CALL SelectHaloMax ! fiind Ltrue ! CALL GetDistribution !Vcmin =380. ; Vcmax = 420. Vcmin =500. ; Vcmax = 700. CALL GetSubhalos(Vcmin,Vcmax) stop end PROGRAM HaloFunction !-------------------------------------------------------------------- ! make sample larger by replicating ! some halos using periodical ! boundary conditions ! dR is the width of the buffer in Mpc/h ! Nhext = new number of halos SUBROUTINE ExpandSample !-------------------------------------------------------------------- use SetUp Logical inx,iny,inz Xleft =-dR ! new boundaries for the catalog Xright =Box+dR Yleft =-dR Yright =Box+dR Zleft =-dR Zright =Box+dR write (*,*) Xleft ,Xright Nhext = Nh Do ip =1,Nh ! go over all halos x =xc(ip) y =yc(ip) z =zc(ip) Do k=-1,1 ! go over all 27 periodical images dz = k*Box inz =(z+dz.gt.Zleft.and.z+dz.lt.Zright) Do j=-1,1 dy = j*Box iny =(y+dy.gt.Yleft.and.y+dy.lt.Yright) Do i=-1,1 dx =i*Box If(i**2+j**2+k**2.ne.0)Then inx =(x+dx.gt.Xleft.and.x+dx.lt.Xright) If(inx.and.iny.and.inz)Then ! add this halo to the list Nhext =Nhext +1 if(Nhext.gt.Nmax) STOP 'too many halos' xc(Nhext) = x+dx yc(Nhext) = y+dy zc(Nhext) = z+dz vxc(Nhext) = vxc(ip) vyc(Nhext) = vyc(ip) vzc(Nhext) = vzc(ip) aMass(Nhext) = aMass(ip) vc(Nhext) = vc(ip) Rvir(Nhext) = Rvir(ip) EndIf ! end inside box EndIf ! end exclude central box EndDo ! end i EndDo ! end j EndDo ! end k EndDo ! end ip write (*,*) ' Nhalos in extended sample=',Nhext write (2,'(" Nhalos in extended sample=",T30,i8)')Nhext RETURN END SUBROUTINE ExpandSample !-------------------------------------------------------------------- ! select isolated halos ! Liso = 1 if halo is isolated SUBROUTINE Isolated3D !-------------------------------------------------------------------- use SetUp aMav =0. rvirav=0. vcav =0. Do ip =1,Nh ! Initialize counters Liso(ip) = 1 EndDo Nx =0 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ip,jp,x,y,z,vv,rr,dd) Do ip =1,Nh ! Test against other halos if(mod(ip,10000).eq.0)write(*,*) ' isolated=',ip x =xc(ip) y =yc(ip) z =zc(ip) vv=vc(ip) rr=Rvir(ip) Do jp =1,Nhext if(jp.ne.ip)Then ! distance dd = ((xc(jp)-x)**2+(yc(jp)-y)**2+(zc(jp)-z)**2) If(dd.lt.(Rvir(jp)/1.e3)**2 & .and.vc(jp).gt.vv)Liso(ip) =0 EndIf ! jp=ip EndDo ! jp loop EndDo ! ip loop Do ip =1,Nh If(Liso(ip).eq.1)Then Nx =Nx +1 aMav =aMav +aMass(ip) rvirav=rvirav +Rvir(ip) vcav =vcav +vc(ip) EndIf EndDo write(*,'(" Nisolated=",3i7)') Nx write(2,'(" Nisolated halos =",T20,i7, & 3x,"=",g12.4," rvir=",f7.1," =",f7.1)') & Nx,aMav/Nx, rvirav/Nx, vcav/Nx RETURN End SUBROUTINE Isolated3D !-------------------------------------------------------------------- ! select halos with close distance to maxima ! Ltrue = 1 if halo is close to maximum SUBROUTINE SelectHaloMax !-------------------------------------------------------------------- use SetUp aMav =0. rvirav=0. vcav =0. Ltrue =0 Nx =0 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ip,jp,x,y,z,vv,rr,dd) Do ip =1,Nh ! Test against other halos if(mod(ip,10000).eq.0)write(*,*) ' isolated=',ip !If(Liso(ip) == 1)Then x =xc(ip) y =yc(ip) z =zc(ip) vv=vc(ip) rr=Rvir(ip) !Liso(ip) = -2 ! flip isolation: Liso=1 only if close to a Maximum Do jp =1,Nmaxima ! distance dd = ((Xmax(jp)-x)**2+(Ymax(jp)-y)**2+(Zmax(jp)-z)**2) If(dd.lt.(Rvir(jp)/1.e3)**2)Ltrue(ip) =1 EndDo ! jp loop !EndIf EndDo ! ip loop Do ip =1,Nh If(Liso(ip)==1.and.Ltrue(ip)==1)Then Nx =Nx +1 aMav =aMav +aMass(ip) rvirav=rvirav +Rvir(ip) vcav =vcav +vc(ip) EndIf EndDo write(*,'(" Nisolatedand Maxima =",3i7)') Nx write(2,'(" Nisolated halos =",T20,i4, & 3x,"=",g12.4," rvir=",f7.1," =",f7.1)') & Nx,aMav/Nx, rvirav/Nx, vcav/Nx Do ip =1,Nh ! Test against other halos If(Liso(ip) == 1 .and. Ltrue(ip) == 1 .and. abs(aMass(ip)-1.e14)<3.e13)Then write(*,20)' isolated:',xc(ip),yc(ip),zc(ip), & vxc(ip),vyc(ip),vzc(ip),aMass(ip), & Rvir(ip),Vrms(ip),vc(ip) EndIf EndDo 20 format(a,3f9.4,3x,3f9.3,g12.4,3f8.2) RETURN End SUBROUTINE SelectHaloMax !-------------------------------------------------------------------- ! Construct averaged profiles of isolated halos ! SUBROUTINE GetDistribution !-------------------------------------------------------------------- use SetUp real*4, DIMENSION(-Nshells:Nshells) :: rad rdsh =0. ; aMsh=0. ; Vdsh =0. ; Vcsh =0. ; nsh =0 dv = 0.1 Nx =0 ! count halos. Find number of bins Nb aMav =0. rvirav =0. vcav =0. vMax =0. aMax = 0. Do ip =1,Nh If(Liso(ip)==1.and. Ltrue(ip)==1)Then Nx =Nx +1 aMav =aMav +aMass(ip) rvirav=rvirav +Rvir(ip) vcav =vcav +vc(ip) vMax =max(vMax,vc(ip)) aMax =max(aMax,aMass(ip)) EndIf EndDo aMav =aMav /Nx rvirav =rvirav/Nx vcav =vcav /Nx write(*,*) ' Max Vc =',vMax,Nx write(*,*) ' Max M =',aMax Do ip =1,Nh ! get profiles If(Liso(ip).eq.1)Then ii =min(INT(log10(aMass(ip))/dv),Nshells) Vcsh(ii) = Vcsh(ii) +vc(ip) aMsh(ii) = aMsh(ii) +aMass(ip) nsh(ii) = nsh(ii) +1 EndIf EndDo ! print results write(2,'(3x,a,T15,a,T25,a,T35,a,T70,a,T85,a,T100,2a)') & 'Vcirc','Nhalos','MdN/dM/V','Mass' iFlag =0 Do j=1,Nshells nn = max(nsh(j),1) nup = 0 do k=j,Nshells nup =nup+nsh(k) enddo aMsh(j) = aMsh(j)/nn Vcsh(j) = Vcsh(j)/nn if(nsh(j).ne.0)iFlag=1 If(iFlag==1)Then If(nn<3)Vcsh(j) =10.**((j+0.5)*dv) vv =10.**((j+1)*dv) - 10.**((j)*dv) if(Vcsh(j) Nmax)Stop ' Error: too few Nmax for the list of maxima' Xmax(Nmaxima) =x Ymax(Nmaxima) =y Zmax(Nmaxima) =z EndDo write(*,*) ' Number of maxima =',Nmaxima RETURN End SUBROUTINE GetMaxima !-------------------------------------------------------------------- ! Construct subhalo Vel function ! SUBROUTINE GetSubhalos(Vcmin,Vcmax) !-------------------------------------------------------------------- use SetUp real*4, DIMENSION(-Nshells:Nshells) :: rads,dns rdsh =0. ; aMsh=0. ; Vdsh =0. ; Vcsh =0. ; nsh =0 rads =0.; dns =0. dv = 0.025; dlogr =0.1 Rout = 2. ! outer radius for satellites profile in rvir units Nx =0 ! count halos. Find number of bins Nb aMav =0. rvirav =0. vcav =0. vMax =0. aMax = 0. Do ip =1,Nh If(Liso(ip).eq.1.and.vc(ip)>Vcmin.and.vc(ip)Vc)','Mass/Mparent' Do j=0,Nshells nn = max(nsh(j),1) aMsh(j) = aMsh(j)/nn Vcsh(j) = Vcsh(j)/nn If(nsh(j) == 0)Vcsh(j)=dv*(j+0.5) nsum = 0 i = j Do while (dv*i< 1.) nsum =nsum +nsh(i) i= i +1 EndDo if(Vcsh(j)< 1.) & write(2,'(f8.3,3x,i6,2x,4g13.4,3x,i8)') Vcsh(j),nsh(j), & float(nsh(j))/Nx, float(nsum)/Nx,& aMsh(j),dv*j EndDo ! print results for density profile write(0,'(3x,a,T15,a,T25,a,T35,a,T50,a,T85,a,T100,2a)') & 'R/Rvir','dN/dVol',' n' iFlag =0 Do j=-Nshells,Nshells nn = max(dns(j),1.) If(dns(j) > 0.5)iFlag =1 If(iFlag ==1)Then r0 =10.**(j*dlogr) r1 =10.**((j+1)*dlogr) dVol = 4.188*(r1**3-r0**3) radius = max(rads(j)/nn,(r0+r1)/2.) If(radius< Rout)& write(0,'(f8.3,3x,g12.4,3x,i6,2x,2g13.4,3x,i8)') radius,dns(j)/dVol,INT(dns(j)),r0,r1,j EndIf EndDo RETURN End SUBROUTINE GetSubhalos