!---------------------------------------------------------------------------------------- ! 3d properies of satellites of Isolated halos ! ! Input: Halo Catalog (short list) ! List of distinct halos ! Output: statistics of subhalos ! Parameters: ! ! MODULE setHalos PUBLIC Real, PARAMETER :: & dLog =0.07, & ! bin size in log(r) dMlog =0.15, & dVlog =0.040 Integer, PARAMETER :: Nbin = 40 Integer, ALLOCATABLE :: LDistinct(:),Lisolate(:) Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc Integer :: Ncenter,Nmain,Ntotal,nnx Real :: Box, Risomax, Vcmin,Vcmax,Vfrac, dBuffer, & Vc_av,aM_av,Rv_av Real, Dimension(-Nbin:Nbin) :: Vc_bin Integer, Dimension(-Nbin:Nbin) :: iVc_bin Contains !-------------------------------------------------------------- ! add buffer of periodically replicated halos ! Width of the buffer is dBuffer ! Number of particles in the expanded sample is Ntotal SUBROUTINE Expand !-------------------------------------------------------------- Logical :: inside,insidex,insidey,insidez Nsample = Nmain !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ic,i,j,k,x,y,z,x1,y1,z1, inside,insidex,insidey,insidez) Do ic=1,Nmain x = Xc(ic); y = Yc(ic); z = Zc(ic) do k=-1,1 z1 = z+k*Box If(z1.le.Box+dBuffer.and.z1.ge.-dBuffer)Then insidez =.true. Else insidez =.false. EndIf do j=-1,1 y1 = y+j*Box If(y1.le.Box+dBuffer.and.y1.ge.-dBuffer)Then insidey =.true. Else insidey =.false. EndIf do i=-1,1 x1 = x+i*Box If(x1.le.Box+dBuffer.and.x1.ge.-dBuffer)Then insidex =.true. Else insidex =.false. EndIf inside = insidex.and.insidey.and.insidez If(abs(i) +abs(j) +abs(k) .eq. 0)inside=.false. If(inside)Then !$OMP critical Nsample = Nsample +1 If(Nsample.le.Ncenter)Then Xc(Nsample) = x1 Yc(Nsample) = y1 Zc(Nsample) = z1 VXc(Nsample) = VXc(ic) VYc(Nsample) = VYc(ic) VZc(Nsample) = VZc(ic) Amc(Nsample) =Amc(ic) Rmc(Nsample) =Rmc(ic) Vcirc(Nsample) =Vcirc(ic) Vrms(Nsample) =Vrms(ic) LDistinct(Nsample) = LDistinct(ic) EndIf !$OMP end critical Endif EndDo EndDo EndDo EndDo write(*,*) 'Buffer width=',dBuffer,' New number of halos=',Nsample write(*,*) ' Maximum number of halos=',Ncenter If(Nsample > Ncenter) Stop ' Error: Too many centers. Increase fExpand ' Ntotal = Nsample End SUBROUTINE Expand !-------------------------------------------------------------- ! 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) ! Do i=1,Ntotal ! Lisolate(i) = 1 ! LDistinct(i) = 1 ! EndDo !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (i,Vcc) & !$OMP REDUCTION(+:Nsample) Do i=1,Nmain Vcc = Vcirc(i) !If(Vcc >Vcmin .and. Vcc < Vcmax.and.aMc(i)< 2.2e12.and.aMc(i)>1.1e12)Then If(Vcc >Vcmin .and. Vcc < Vcmax)Then Nsample =Nsample +1 Lisolate(i) =1 Else Lisolate(i) =0 EndIf EndDo !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ic,jc,x,y,z,R0,R1,Amass,Vcc,dx,dy,dz,dd) Do ic=1,Nmain !If(mod(ic,50000)==0)write(*,*)' Distinct: ',ic x = Xc(ic); y = Yc(ic); z = Zc(ic) R0 = Rmc(ic) Vcc = Vcirc(ic) Amass = aMc(ic) If(Lisolate(ic) == 1)Then !write(*,'(a,3x,3f8.3,2x,1P,3g12.3)') ' ------- X=',x,y,z,sqrt(R0),Vcc,Amass Do jc =1,Ntotal !If(jc /= ic.and.Vcirc(jc)>=Vcc/1.5)Then ! compare with large halos If(jc /= ic.and.Vcirc(jc)>=Vcc/Vfrac)Then ! compare with large halos dx = Xc(jc)-x dy = Yc(jc)-y dz = Zc(jc)-z R1 = (max(R0,Rmc(jc))*Risomax)**2 dd = dx**2 +dy**2 +dz**2 !If(dd< R1 .and.(Vcirc(jc)>=Vcc/Vfrac.or.dd>R0**2))Lisolate(ic) =0 ! not isolated If(dd< R1)Lisolate(ic) =0 ! not isolated EndIf ! Vcirc End Do ! jc EndIf ! mass EndDo ! ic write(*,*) 'Done Distinct' nnx=0 Vc_av =0. aM_av =0. Rv_av = 0. Do i=1,Nmain If(Lisolate(i).eq.1)Then nnx=nnx+1 Vc_av = Vc_av + Vcirc(i) aM_av =aM_av + aMc(i) Rv_av =Rv_av + Rmc(i) EndIf EndDo Vc_av = Vc_av/nnx aM_av =aM_av/nnx Rv_av =Rv_av/nnx write(*,*) 'Number of isolated halos=',nnx write(*,*) 'Average Vcirc =',Vc_av write(*,*) 'Averge Mvir =',aM_av write(*,*) 'Averge Rvir =',Rv_av Vvir = 6.582e-5*sqrt(aM_av/Rv_av) write(*,*) 'Averge Vvir =',Vvir !Do i =1,Ntotal ! write(30,'(i8,3f10.4,g12.4,i3)')i,Xc(i),Yc(i),Zc(i),aMc(i),Lisolate(i) !EndDo Return End SUBROUTINE HaloIsolated !-------------------------------------------------------------- ! Find satellites and ! get their properties SUBROUTINE Satellites !-------------------------------------------------------------- Real :: Vrad(-Nbin:Nbin), Vtan(-Nbin:Nbin), & Vrad2(-Nbin:Nbin), rad(-Nbin:Nbin) Integer :: Nin(-Nbin:Nbin) Vrad = 0. Vrad2 = 0. Vtan = 0. Nin = 0 rad = 0 Rmax0 = dBuffer**2 Nisolate =0 Vc_bin = 0. ; iVc_bin = 0 ! count satellites !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ic,x,y,z,R0,aM0,Vcc,Vvir,jc,dx,dy,dz,dd, & !$OMP dVx,dVy,dVz,Vv,Vr,Vt,indx,V2,Vt2) & !$OMP REDUCTION(+:Nisolate) Do ic=1,Nmain If(Lisolate(ic).eq.1)Then ! isolated halo Nisolate = Nisolate + 1 x = Xc(ic); y = Yc(ic); z = Zc(ic) R0 =Rmc(ic) aM0 =Amc(ic) Vcc =Vcirc(ic) Vvir = 6.582e-5*sqrt(aM0/R0) !write (30,'(2i9,3f10.4,g12.4,2f9.2)') ic,Nisolate, X, Y , Z,aM0,Vcc Do jc=1,Ntotal If(jc /= ic )Then dx = Xc(jc)-x dy = Yc(jc)-y dz = Zc(jc)-z dd = dx**2 +dy**2 +dz**2 If(dd VdN/dV dN N(>V)/halo N(>V)total' indx = 0 Do i =-Nbin,Nbin If(iVc_bin(i) /= 0.and.indx == 0)indx = 1 V = 10.**((i)*dVlog) V1 = 10.**((i+1)*dVlog) If(indx==1 .and.V<=1.)Then Vv = Vc_bin(i) /max(iVc_bin(i),1) aN = iVc_bin(i)/float(Nisolate)/(V1-V)*(V1+V)/2. If(Vv<1.e-5)Vv = V write(1,'(2f9.4,f10.4,i8,f8.3,i8)')V,Vv,aN,iVc_bin(i),Nin(i)/float(Nisolate),Nin(i) EndIf EndDo Return End SUBROUTINE Satellites !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadHalos(iFlag,vcLimit) !-------------------------------------------------------------- Character :: Line*80, Txt*6='Number', Txt2*8,txt1*4 Ncenter =0 Rdmax =0. jSecond = 0 Do i=1,2 Read(3,'(a)')Line If(iFlag==0)write(*,'(a)') Line EndDo Read(3,'(a4,i4)')txt1 Read(3,'(a4,i4)')txt1 If(iFlag == 1)Then jj =0 Do i=1,2 Read(4,'(a)')Line write(*,'(a)') Line EndDo Read(4,'(a4,i4)')txt1,ii Do read(4,*,err=20)Txt2 !write(*,*)TRIM(Txt2) If(TRIM(Txt2)==Txt)Then jj =jj +1 If(jj==2)Exit EndIf EndDo EndIf Do !read (3,45,iostat=iStat) X0, Y0 , Z0, VvX,VvY,VvZ, & ! aM,rr,vrmss,vcrc,Npart,Rmax,Conc !,Nb read (3,*,iostat=iStat) j, X0, Y0 , Z0, & aM,aMmax,Vcrc,VcrcMax,iSub,iSubever ! 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 Amc(Ncenter) = aM !max Vcirc(Ncenter) = Vcrc !Max Do read (4,45,iostat=iStat) X1, Y1 , Z1, VvX,VvY,VvZ, & aM1,rr1,vrmss1,vcrc1,Npart1,Rmax1,Conc1 If(iStat.ne.0)stop ' error reading second catalog' jSecond = jSecond +1 if(jSecond == j)Then VXc(Ncenter) = VvX VYc(Ncenter) = VvY VZc(Ncenter) = VvZ Rmc(Ncenter) = rr1/1000. Vrms(Ncenter) = vrmss1 Rdmax = MAX(rr1/1000.,Rdmax) If(abs(X1-X0)>1.e-3.or.abs(Y1-Y0)>1.e-3.or.abs(Z1-Z0)>1.e-3)Then write(*,*) ' ----- Error in matching halo catalogs: j, jSecond=',j,jSecond Stop EndIf exit EndIf EndDo 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( 4,file='CatshortA.416.DAT',Status='OLD')! short list of halos Open( 3,file='CatshortB.416.DAT',Status='OLD')! short list of halos !Open( 2,file='Catdistinct.416.DAT')! short list of halos Box =250. Risomax = 3.0 dBuffer = 15. Vcmin = 195. Vcmax = 205. write(*,'(a,$)') ' Enter Min and Max for halo Vcirc:' read(*,*) Vcmin,Vcmax Vfrac = 1. Vclimit = 80. AEXPN =1.0 ISTEP = 416 write(file1,'(a,i4.4,a,i4.4,a)') 'BolshoiDistA.',INT((Vcmin+Vcmax)/2.),'.',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 Nmain = Ncenter fExpand = (1.+2.*dBuffer/Box)**3*1.1 ! estimate how much more space ! will be needed Ncenter = Nmain*fExpand ! make more centers for periodical conditions write(*,*) ' fExpand =',fExpand write(*,*) ' Ncenter =',Ncenter 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) Ncenter = Nmain*fExpand ! restore Ncenter Call Expand ! make more halos by adding halos ! due to periodical conditions ALLOCATE (Lisolate(Ncenter)) Call HaloIsolated ! make list of distinct halos close(2) Vvir = 6.582e-5*sqrt(aM_av/Rv_av) write (*,*) 'Found Isolated halos' write(1,'(a,f8.4,a,i4,a,i9,a,f8.2,/a,i8,a,f8.3,a,p1g12.3,0p,a,f8.2/,a,f8.3,a,2f8.3,a,f8.3,a,f6.2,/a,a)') & ' a=',AEXPN,' Step=',ISTEP, & ' Number of halos =',Nmain, & ' Velocity limit =',Vclimit, & ' Number of isolated halos=',nnx, & ' Average Vcirc =',Vc_av, & ' Average Mvir =',aM_av, & ' Average Rvir =',Rv_av, & ' Averge Vvir =',Vvir, & ' Isolation: Halo Vcirc=',Vcmin,Vcmax, & ' Vother/Vhalo>',1./Vfrac, & ' Dist Isolate/Rvir=',Risomax, & ' N R/Rvir_bin Vrad/Vvir Vtanrms/Vvir', & ' Vrad_rms/Vvir V_total/Vvir beta Number-density' Risomax =5.0 Call Satellites ! get statistics of satellites End PROGRAM HaloSatellites