!---------------------------------------------------------------------------------------- ! Find properties of large distinct halos ! Do not remove unbound particles ! Input: Halo Catalog (short list) ! Dark matter particles ! Output: Halo catalog ! Parameters: ! outer radius in units of virial radius ! number of particles for potential energy ! MODULE setHalos PUBLIC Real, PARAMETER :: & dLog = 0.075, & ! bin size in log(r) dMlog =0.15 Integer, ALLOCATABLE :: LDistinct(:), LDistinctX(:),LDistinctY(:),LDistinctZ(:) Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc Integer :: Ncenter,Nmx,Nbx,Nmy,Nby,Nmz,Nbz Real :: Box, Risomax, Vmax,Ammin,Ammax,Vfrac Contains !-------------------------------------------------------------- ! Make list of isolated halos ! Input: Ncenter - numer of halos, {Xc ...} = halo data ! ! Output: LDistinct(i) = 1 for isolated halo ! SUBROUTINE HaloIsolated !-------------------------------------------------------------- Nsample = 0 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,Amass) & !$OMP REDUCTION(+:Nsample) Do i=1,Ncenter Amass = aMc(i) If(Amass >Ammin .and. Amass < Ammax.and.LDistinct(i)>=0)Then Nsample =Nsample +1 LDistinctX(i) =1 LDistinctY(i) =1 LDistinctZ(i) =1 Else LDistinctX(i) =0 LDistinctY(i) =0 LDistinctZ(i) =0 EndIf EndDo write(*,*) 'Box=',Box,' Halos in the mass range=',Nsample !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ic,jc,x,y,z,R0,R1,Amass,Vcc,dx,ux,dy,uy,dz,uz, & !$OMP signn,dd,dV) Do ic=1,Ncenter !If(mod(ic,1000)==0)write(*,*)' Distinct: ',ic x = Xc(ic); y = Yc(ic); z = Zc(ic) R0 = Rmc(ic) ! square of projected max radius Vcc = Vcirc(ic) Amass = aMc(ic) If(Amass >Ammin .and. Amass < Ammax)Then !write(*,'(a,3x,3f8.3,2x,1P,3g12.3)') ' ------- X=',x,y,z,sqrt(R0),Vcc,Amass Do jc =1,Ncenter If(jc /= ic.and.Vcirc(jc)>=Vcc/Vfrac)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)) R1 = (max(R0,Rmc(jc))*Risomax)**2 ! along X-direction signn =1. If(ux>Box/2..or.(ux>-Box/2..and.ux<0.))signn =-1. dd = dy**2 +dz**2 dV = abs(signn*(Vxc(jc) -Vxc(ic)) + 100.*dx) !If(dd< R0)write(*,'(a,3x,4f8.3,2x,1P,8g12.3)') ' ',ux,dx,Xc(jc),x,sqrt(dd),dV,Vxc(jc),Vxc(ic) If(dd< R1 .and. dV < Vmax) LDistinctX(ic) =0 ! not Isolated ! along Y-direction signn =1. If(uy>Box/2..or.(uy>-Box/2..and.uy<0.))signn =-1. dd = dx**2 +dz**2 dV = abs(signn*(Vyc(jc) -Vyc(ic)) + 100.*dy) If(dd< R1 .and. dV < Vmax) LDistinctY(ic) =0 ! not Isolated ! along Z-direction signn =1. If(uz>Box/2..or.(uz>-Box/2..and.uz<0.))signn =-1. dd = dx**2 +dy**2 dV = abs(signn*(Vzc(jc) -Vzc(ic)) + 100.*dz) If(dd< R1 .and. dV < Vmax) LDistinctZ(ic) =0 ! not Isolated EndIf ! Vcirc End Do ! jc EndIf ! mass EndDo ! ic write(*,*) 'Done LDistinct' nnx=0 nny=0 nnz=0 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i) REDUCTION(+:nnx,nny,nnz) Do i=1,Ncenter If(LDistinctX(i).eq.1)nnx=nnx+1 If(LDistinctY(i).eq.1)nny=nny+1 If(LDistinctZ(i).eq.1)nnz=nnz+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 =',nnx,nny,nnz Return End SUBROUTINE HaloIsolated !-------------------------------------------------------------- ! Find ad write dV -dR for satellites ! SUBROUTINE Satellites !-------------------------------------------------------------- Rmax0 = Risomax**2 ! count satellites !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ic,jc,x,y,z,R0,aM0,Vcc,dx,ux,dy,uy,dz,uz,dd,signn,dV) Do ic=1,Ncenter x = Xc(ic); y = Yc(ic); z = Zc(ic) R0 =Rmc(ic) aM0 =Amc(ic) Vcc =Vcirc(ic) ! ------ X-direction If(LDistinctX(ic).eq.1)Then ! isolated halo Do jc=1,Ncenter If(jc /= ic .and. LDistinct(jc).ne.-1)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)) ! along X-direction signn =1. If(ux>Box/2..or.(ux>-Box/2..and.ux<0.))signn =-1. dd = dy**2 +dz**2 dV = abs(signn*(Vxc(jc) -Vxc(ic)) + 100.*dx) If(ddBox/2..or.(uy>-Box/2..and.uy<0.))signn =-1. dd = dx**2 +dz**2 dV = abs(signn*(Vyc(jc) -Vyc(ic)) + 100.*dy) If(ddBox/2..or.(uz>-Box/2..and.uz<0.))signn =-1. dd = dx**2 +dy**2 dV = abs(signn*(Vzc(jc) -Vzc(ic)) + 100.*dz) If(dd 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 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) If(iFlag ==1)Then Do i=1,Ncenter read(2,*) j,LDistinct(i) if(j/=i)Stop ' error in Distinct halos' EndDo EndIf 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 ! ------------------------ 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 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='Catdistinct.416.DAT')! short list of halos Box =250. Risomax = 1.5 Vmax = 1500. Ammin = 2.0e12 Ammax = 3.0e12 write(*,'(a,$)') ' Enter Min and Max for halo masses:' read(*,*) Ammin,Ammax Vfrac = 2. Vclimit = 60. AEXPN =1.0 ISTEP = 416 write(file1,'(a,i4.4,a,i4.4,a)') 'BolshoiSat.',INT((Ammin+Ammax)/2./1.e11),'.',416,'.DAT' Open( 1,file=file1)! short list of halos ! write(*,'(a,$)')' Enter limit on circular velocity =' ! read(*,*)vcLimit 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 (LDistinct(Ncenter)) write (*,*) 'Allocated Centers' CALL ReadHalos(1,vcLimit) ! read halos Close(3) ALLOCATE (LDistinctX(Ncenter),LDistinctY(Ncenter),LDistinctZ(Ncenter)) Call HaloIsolated ! make list of distinct halos close(2) write (*,*) 'Found Isolated halos' write(1,'(a,f8.4,a,i4,/a,1P,2g11.3,0P,a,f6.3,/10x,a,f6.3,a,f6.1,/2a)') & ' a=',AEXPN,' Step=',ISTEP, & ' Isolation: Halo Masses=',Ammin,Ammax, & ' Vother/Vhalo>',1./Vfrac, & ' Dist Isolate/Rvir=',Risomax,' dV(km/s)<',Vmax, & ' dR(kpch) dV(km/s) MassSat(Msunh) Vcirc(km/s)', & ' MassPrime RvirPrime(kpc)' Risomax =1.5 Vmax =1500. Call Satellites ! make list of satellites ! write (*,*) 'Found Satellites' End PROGRAM HaloSatellites