!---------------------------------------------------------------------------------------- ! Find properties of large distinct halos ! ! Input: Halo Catalog (short list) ! ! MODULE setHalos PUBLIC Real, PARAMETER :: & dLogV = 0.075, & ! bin size in log(Vc) dLogR = 0.075, & ! bin size in log(r) dMlog =0.15, & Risolate = 1. ! Radius of isolation in Rvir units Integer, ALLOCATABLE :: LDistinct(:),Lst(:),Label(:,:,:) Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc Integer :: Ncenter,Nmx,Nbx,Nmy,Nby,Nmz,Nbz Real :: Box, Cell , Rmax ,Vlim1, Vlim2 Contains !-------------------------------------------------------------- ! Make list of distinct halos ! Input: Ncenter - numer of halos, {Xc ...} = halo data ! ! Output: LDistinct(i) = 1 for distinct halo ! SUBROUTINE HaloDistinct !-------------------------------------------------------------- !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Ncenter LDistinct(i) =1 EndDo write(*,*) 'Box=',Box !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ic,jc,x,y,z,R0,Vcc,dx,ux,dy,uy,dz,uz, & !$OMP N0x,N1x,N0y,N1y,N0z,N1z, i,i1,j,j1,k,k1) Do ic=1,Ncenter !If(mod(ic,100000)==0)write(*,*)' Distinct: ',ic x = Xc(ic); y = Yc(ic); z = Zc(ic) R0 =Rmc(ic) Vcc =Vcirc(ic) N0x = (x-Rmax)/Cell ; N1x = (x+Rmax)/Cell ! limits for linker list N0y = (y-Rmax)/Cell ; N1y = (y+Rmax)/Cell N0z = (z-Rmax)/Cell ; N1z = (z+Rmax)/Cell Do i =N0x,N1x i1 =mod(i+50*(Nbx+1),Nbx+1) ! periodical boundaries Do j =N0y,N1y j1 =mod(j+50*(Nby+1),Nby+1) ! periodical boundaries Do k =N0z,N1z k1 =mod(k+50*(Nbz+1),Nbz+1) ! periodical boundaries jc = Label(i1,j1,k1) Do while (jc.ne.0) ! loop over all particles in cell If(jc /= ic.and.Vcirc(jc)>=Vcc)Then ux = Xc(jc)-x uy = Yc(jc)-y uz = Zc(jc)-z ! periodical boundary conditions dx = min(abs(ux),abs(ux+Box),abs(uy-Box)) dy = min(abs(uy),abs(uy+Box),abs(uy-Box)) dz = min(abs(uz),abs(uz+Box),abs(uz-Box)) dd =dx**2 +dy**2 +dz**2 If(dd< (Risolate*Max(Rmc(jc),R0))**2) Then If(Vcirc(jc).gt.Vcc.or.jc>ic)LDistinct(ic) =0 ! subhalo If(Vcc.gt.0.95*Vcirc(jc))LDistinct(ic) =-1 ! Fake EndIf ! dd R0 jc = Lst(jc) End Do End Do ! k End Do ! j End Do ! i EndDo ! ic write(*,*) 'Done LDistinct' nn=0 nf =0 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) REDUCTION(+:nn,nf) Do i=1,Ncenter If(LDistinct(i).eq.1)nn=nn+1 If(LDistinct(i).eq.-1)nf=nf+1 !if(Vcirc(i)>1600.)write(*,'(i8,i3,f8.2,3x,3f9.3)')i,LDistinct(i),Vcirc(i),Xc(i),Yc(i),Zc(i) EndDo write(*,*) 'Number of distinct halos =',nn,' Fake=',nf Do i=1,Ncenter Write(2,'(i10,i3,1P,8g12.4)')i,LDistinct(i),aMc(i),Vcirc(i),Xc(i),Yc(i),Zc(i),Rmc(i) EndDo Return End SUBROUTINE HaloDistinct !-------------------------------------------------------------- ! Mass function of distinct halos ! SUBROUTINE MassDistinct !-------------------------------------------------------------- Integer, Dimension(-100:100) :: & mHalos,iHalos Real, Dimension(-100:100) :: massbin mHalos = 0 massbin = 0. nn=0 ! count all halos Do i=1,Ncenter If(LDistinct(i) == 1)Then nn =nn +1 ind =INT(log10(Amc(i))/dMlog+100)-100 If(ind<100)Then mHalos(ind) = mHalos(ind) +1 massbin(ind) = massbin(ind) +Amc(i) EndIf EndIf EndDo massbin = massbin/(max(mHalos,1)) istart =-100 ! find first and last non-empty bin Do i =-100,100 If(mHalos(i)>0)Then istart =i ; exit EndIf EndDo ifinish =100 Do i =100,-100,-1 If(mHalos(i)>0)Then ifinish =i ; exit EndIf EndDo iHalos = mHalos Do i =ifinish-1,istart,-1 mHalos(i) = mHalos(i) +mHalos(i+1) EndDo write(*,*) ' Mass(Msunh) N(>M) Mass MdN/dM' Do i=istart,ifinish ! print results If(mHalos(i)>0)Then aMm =10.**(i*dMlog) dM = 10.**((i+1)*dMlog)-aMm write(*,'(3x,1P,G11.4,0P,i9,2x,1P,2g11.4)')aMm, mHalos(i),massbin(i),massbin(i)*iHalos(i)/dM EndIf EndDo End SUBROUTINE MassDistinct !-------------------------------------------------------------- ! Find satellite halos ! SUBROUTINE Satellites !-------------------------------------------------------------- Integer, Dimension(-100:100) :: iSat,iHalos,iFake,ibin Real, Dimension(-100:100) :: vbin,rbin iSat = 0 iHalos = 0 vbin = 0. iFake = 0 rbin = 0. ibin = 0 ! count satellites !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ic,jc,x,y,z,R0,aM0,Vcc,dx,ux,dy,uy,dz,uz,dd,ind, & !$OMP N0x,N1x,N0y,N1y,N0z,N1z, i,i1,j,j1,k,k1 ) Do ic=1,Ncenter If(LDistinct(ic).eq.1)Then ! distinct halo x = Xc(ic); y = Yc(ic); z = Zc(ic) R0 =Rmc(ic) aM0 =Amc(ic) Vcc =Vcirc(ic) N0x = (x-Rmax)/Cell ; N1x = (x+Rmax)/Cell ! limits for linker list N0y = (y-Rmax)/Cell ; N1y = (y+Rmax)/Cell N0z = (z-Rmax)/Cell ; N1z = (z+Rmax)/Cell Do i =N0x,N1x i1 =mod(i+50*(Nbx+1),Nbx+1) ! periodical boundaries Do j =N0y,N1y j1 =mod(j+50*(Nby+1),Nby+1) ! periodical boundaries Do k =N0z,N1z k1 =mod(k+50*(Nbz+1),Nbz+1) ! periodical boundaries jc = Label(i1,j1,k1) Do while (jc.ne.0) ! loop over all particles in cell If(jc /= ic)Then ux = Xc(jc)-x uy = Yc(jc)-y uz = Zc(jc)-z ! periodical boundary conditions dx = min(abs(ux),abs(ux+Box),abs(ux-Box)) dy = min(abs(uy),abs(uy+Box),abs(uy-Box)) dz = min(abs(uz),abs(uz+Box),abs(uz-Box)) dd =sqrt(dx**2 +dy**2 +dz**2)/R0 If(ddVlim1.and.Vcc<=Vlim2)Then ind =MIN(INT(log10(dd)/dlogR+100)-100,100) !$OMP atomic rbin(ind) = rbin(ind) +dd !$OMP atomic ibin(ind) = ibin(ind) + 1 EndIf EndIf ! dd< R0 EndIf ! jc /= ic jc = Lst(jc) End Do ! jc End Do ! k End Do ! j End Do ! i EndIf ! LDistinct EndDo ! ic nn=0 ! count all halos Do i=1,Ncenter ind =INT(log10(Vcirc(i))/dlogV+100)-100 If(ind<100)Then If(LDistinct(i)>=0)Then iHalos(ind) = iHalos(ind) +1 vbin(ind) = vbin(ind) +Vcirc(i) Else iFake(ind) =iFake(ind) +1 EndIf EndIf EndDo vbin = vbin/(max(iHalos,1)) !------------------ Velocity function istart =-100 ! find first and last non-empty bin Do i =-100,100 If(iHalos(i)>0)Then istart =i ; exit EndIf EndDo ifinish =100 Do i =100,-100,-1 If(iHalos(i)>0)Then ifinish =i ; exit EndIf EndDo Do i =ifinish-1,istart,-1 iHalos(i) = iHalos(i) +iHalos(i+1) iSat(i) = iSat(i) +iSat(i+1) iFake(i) = iFake(i) +iFake(i+1) EndDo write(*,*) ' Vcirc(km/s) All(>V) Sat(>V) Fake(>V) Sat/All' Do i=istart,ifinish ! print results If(iHalos(i)>0)Then v =10.**(i*dlogV) write(*,'(3x,f9.2,3x,3i7,f8.3)')v, iHalos(i), iSat(i),iFake(i),float(iSat(i))/iHalos(i) EndIf EndDo !------------------ Radial Profile istart =-100 ! find first and last non-empty bin Do i =-100,100 If(ibin(i)>0)Then istart =i ; exit EndIf EndDo ifinish =100 Do i =100,-100,-1 If(ibin(i)>0)Then ifinish =i ; exit EndIf EndDo write(*,*) ' (R/Rvir) dNsat dNsat/dV/Nsat ' Do i=istart,ifinish ! print results r0 =10.**(i*dlogR) r1 =10.**((i+1)*dlogR) Vol = 4.188*(r1**3-r0**3) rr = (r1+r0)/2. r = rr If(ibin(i)>1)r =rbin(i)/ibin(i) write(*,'(3x,f9.4,3x,i7,1P,g12.4)')r,ibin(i),ibin(i)/Vol EndDo Return End SUBROUTINE Satellites !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadHalos(iFlag,vcLimit) !-------------------------------------------------------------- Character :: Line*80, Txt*6='Number', Txt2*28,txt1*4 Rdmax =0. Do i=1,2 Read(3,'(a)')Line If(iFlag==0)write(2,'(a)') Line EndDo Read(3,'(a4,i4)')txt1,ii If(iFlag==0) write(2,'(a4,i4)') txt1,ii If(iFlag==1) write(2,'(a,i8,a,f8.2)')' Halos selected=',Ncenters,' Vc limit=',vcLimit indx =0 Do i=1,30 ! find second appearence of 'Number' read(3,*,err=20)Txt2 !write(*,*)TRIM(Txt2) If(TRIM(Txt2)==Txt)then If(indx ==1)Exit indx=1 EndIf EndDo If(iFlag==0)Ncenter =0 If(iFlag==1)Then !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Ncenter Xc(i) =0. ; Yc(i) =0. ; Zc(i) =0. Vxc(i) =0. ; VYc(i) =0. ; VZc(i) =0. Amc(i) =0. ; Rmc(i) = 0.; Vcirc(i) =0. Vrms(i) =0. EndDo EndIf Ncenter =0 Do read (3,45,iostat=iStat) X0, Y0 , Z0, VvX,VvY,VvZ, & aM,rr,vrmss,vcrc,Npart,Rmax,Conc !,Nb ! write (*,'(i9,3f10.4,f8.2,i10,2f9.2,i6)') Ncenter, X0, Y0 , Z0,vcrc,Npart,Rmax,Conc,Nb 45 Format(F9.4,2F9.4,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4) If(iStat.ne.0)Exit If(vcrc> vcLimit)Then Ncenter = Ncenter +1 If(iFlag==1)Then Xc(Ncenter) = X0 Yc(Ncenter) = Y0 Zc(Ncenter) = Z0 VXc(Ncenter) = VvX VYc(Ncenter) = VvY VZc(Ncenter) = VvZ Amc(Ncenter) = aM Rmc(Ncenter) = rr/1000. Vcirc(Ncenter) = vcrc ! *sqrt(1.+(20./vcrc)**2) Vrms(Ncenter) = vrmss Rdmax = MAX(rr/1000.,Rdmax) EndIf EndIf ! test mass Enddo 40 Format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,g11.3) If(iFlag==0)write(*,*) 'Read Halos=',Ncenter If(iFlag==1)write(*,'(a,i8,a,f8.3,a)') ' Read Halos=',Ncenter, & ' Max Radius=',Rdmax,'Mpch' rewind(3) Return 10 write(*,*) ' Could not read radial bins' stop 20 write(*,*) ' Could not read text R(kpc/h)' stop 30 write(*,*) ' Unexpected End of file' stop End SUBROUTINE ReadHalos !------------------------------------------------------------------------- ! Make linker list of Halos SUBROUTINE ListHalos !------------------------------------------------------------------------ !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) Do i=1,Ncenter Lst(i)= -1 EndDo !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,j,k) Do k=Nmz,Nbz Do j=Nmy,Nby Do i=Nmx,Nbx Label(i,j,k)=0 EndDo EndDo EndDo Do jp=1,Ncenter i=INT(Xc(jp)/Cell) j=INT(Yc(jp)/Cell) k=INT(Zc(jp)/Cell) i=MIN(MAX(Nmx,i),Nbx) j=MIN(MAX(Nmy,j),Nby) k=MIN(MAX(Nmz,k),Nbz) Lst(jp) =Label(i,j,k) Label(i,j,k) =jp EndDo write(*,*) ' Done Linker List' Return End SUBROUTINE ListHalos ! ------------------------ real function seconds () ! ------------------------ ! ! purpose: returns elapsed time in seconds Integer, SAVE :: first=0,rate=0 Real , SAVE :: t0=0. !real*8 dummy !real tarray(2) If(first==0)Then CALL SYSTEM_CLOCK(i,rate) first =1 t0 =float(i)/float(rate) seconds = 0. Else ! seconds = etime(tarray) CALL SYSTEM_CLOCK(i) seconds = float(i)/float(rate)-t0 EndIf write(*,*) ' time=',seconds,i,t0 return end function seconds end module SetHalos !-------------------------------------------------------------------------------------- ! ! PROGRAM HaloSatellites use SetHalos Character*120 :: file1,file2,catname,listname,txt*8 Open( 3,file='CatshortA.416.DAT',Status='OLD')! short list of halos Open( 2,file='CatdistinctZ.416.DAT')! short list of halos Box =250. Cell = 2.3 Rmax = 2.2 write(*,'(a,$)')' Enter limit on circular velocity =' read(*,*)vcLimit write(*,*) write(*,'(a,$)')' Enter Vcirc range for radial profile =' read(*,*) Vlim1,Vlim2 write(*,*) tt =seconds() CALL ReadHalos(0,vcLimit) ! count halos ALLOCATE (Xc(Ncenter),Yc(Ncenter),Zc(Ncenter)) ALLOCATE (VXc(Ncenter),VYc(Ncenter),VZc(Ncenter)) ALLOCATE (Amc(Ncenter),Rmc(Ncenter),Vrms(Ncenter),Vcirc(Ncenter)) ALLOCATE (Lst(Ncenter)) Nmx =0 ; Nbx = Box/Cell +1 Nmy =0 ; Nby = Box/Cell +1 Nmz =0 ; Nbz = Box/Cell +1 ALLOCATE (Label(Nmx:Nbx,Nmy:Nby,Nmz:Nbz)) tt =seconds() write (*,*) 'Allocated Centers' CALL ReadHalos(1,vcLimit) tt =seconds() ! read halos CALL ListHalos Close(3) tt =seconds() ALLOCATE (LDistinct(Ncenter)) Call HaloDistinct ! make list of distinct halos close(2) tt =seconds() write (*,*) 'Found Distinct halos' Call Satellites ! make list of satellites write (*,*) 'Found Satellites' tt =seconds() Call MassDistinct ! make list of satellites write (*,*) 'Found Mass distribution' tt =seconds() End PROGRAM HaloSatellites