!---------------------------------------------------------------------------------------- ! Match FOF halos with BDM ! ! Input: Halo Catalog from Stefan ! Catshort catalog ! Output: list of matching halos ! ! ! MODULE setHalos PUBLIC Real, PARAMETER :: & dR =1.0, & ! fraction of virial radius for matching Risomax = 1., & ! radius in Rvir units to distinct/isolated halos Vfrac = 1. ! relative velocity for distinct/isolated halos Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc, & Xf,Yf,Zf, & VXf,VYf,Vzf, & Amf,Vrmf,Rsph Integer, ALLOCATABLE, DIMENSION(:) :: Halo,Lisolate,Match Integer :: Ncenter,Nmain,Ntotal,nnx,Nfof Real :: Box, Vcmin,Vcmax, dBuffer, & Vc_av,aM_av,Rv_av 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) Halo(Nsample) =Halo(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 distinct halos ! Input: Ncenter - numer of halos, {Xc ...} = halo data ! ! Output: Lisolate(i) = 1 for distinct halo ! SUBROUTINE HaloIsolated !-------------------------------------------------------------- !$OMP PARALLEL DO DEFAULT(SHARED) Do i=1,Ntotal Lisolate(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/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)Lisolate(ic) =0 ! not isolated EndIf ! Vcirc End Do ! jc EndIf ! mass EndDo ! ic !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (ic) Do ic=Nmain+1,Ntotal Lisolate(ic) = Lisolate(Halo(ic)) EndDo 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(*,*) 'Average Mvir =',aM_av write(*,*) 'Average Rvir =',Rv_av Vvir = 6.582e-5*sqrt(aM_av/Rv_av) write(*,*) 'Averge Vvir =',Vvir Return End SUBROUTINE HaloIsolated !-------------------------------------------------------------- ! Find match of FOF with BDM ! Input: Nfof, Ntotal ! ! Output: write info to files ! SUBROUTINE MatchHalos !-------------------------------------------------------------- Integer, Dimension(1000) :: iMatch Real, Dimension(1000) :: dMatch !$OMP PARALLEL DO DEFAULT(SHARED) Do i=1,Ntotal Match(i) = 0 ! no match initialization EndDo Do ic=1,Nfof ! If(mod(ic,50000)==0)write(*,*)' Distinct: ',ic x = Xf(ic); y = Yf(ic); z = Zf(ic) R0 = (Rsph(ic)*dR)**2 Vcc = Vrmf(ic) Amass = aMf(ic) nMatch = 0 iMatch = 0 dMatch = 0. !write(*,'(a,3x,3f8.3,2x,1P,3g12.3)') ' ------- X=',x,y,z,sqrt(R0),Vcc,Amass Do jc =1,Ntotal If(Lisolate(jc) == 1)Then ! compare with large halos dx = Xc(jc)-x dy = Yc(jc)-y dz = Zc(jc)-z dd = dx**2 +dy**2 +dz**2 If(dd< R0)Then ! found match If(Match(jc) == 1)Then write(*,'(a,2i8,3x,6f10.4,1P,2g13.4)') ' already matched:',ic,jc, & x,y,z,Xc(jc),Yc(jc),Zc(jc),Amass,aMc(jc) Else nMatch = nMatch +1 iMatch(nMatch) = jc dMatch(nMatch) = sqrt(dd)*1000. EndIf ! match EndIf ! dd < R0 EndIf ! Lisolate End Do ! jc If(nMatch == 0)Then write(20,'(i8,3x,3f10.4,1P,g13.4)') ic,x,y,z,Amass Else jj = 0 xmax = 0. Do j = 1,nMatch if(aMc(iMatch(j))>xMax)Then jj = j ; xMax = aMc(iMatch(j)) EndIf EndDo if(jj==0)Stop ' error in max halo' iHalo = iMatch(jj) Match(iHalo) = 1 write(2,'(2i8,3x,7f10.4,1P,2g13.4)') ic,iHalo, & x,y,z,Xc(iHalo),Yc(iHalo),Zc(iHalo),dMatch(jj),Amass,aMc(iHalo) EndIf EndDo ! ic write(*,*) 'Done Matching' nnx=0 Vc_av =0. aM_av =0. Rv_av = 0. Do i=1,Nmain If(Lisolate(i).eq.1.and.Match(i)==0)Then nnx=nnx+1 Vc_av = Vc_av + Vcirc(i) aM_av =aM_av + aMc(i) Rv_av =Rv_av + Rmc(i) write(30,'(i8,3x,3f10.4,1P,g13.4)') i, & Xc(i),Yc(i),Zc(i),aMc(i) EndIf EndDo Vc_av = Vc_av/nnx aM_av =aM_av/nnx Rv_av =Rv_av/nnx write(*,*) 'Number of BDM missed halos=',nnx write(*,*) 'Average Vcirc =',Vc_av write(*,*) 'Average Mvir =',aM_av write(*,*) 'Average Rvir =',Rv_av Vvir = 6.582e-5*sqrt(aM_av/Rv_av) write(*,*) 'Averge Vvir =',Vvir Return End SUBROUTINE MatchHalos !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadHalos(iFlag,Vclimit) !-------------------------------------------------------------- Character :: Line*120, Txt*8='Coordina', Txt2*8,txt1*4 Ncenter = 0; aM_a = 0. Vc_a = 0. Do i=1,25 Read(3,'(a)')Line !If(i<10)write(2,'(a)')Line write(*,'(a)')Line EndDo Do ! kk=1,50000 read (3,*,iostat=iStat) X0, Y0 , Z0, VvX,VvY,VvZ, & aM,rr,vrmss,vcrc !write (*,'(2i9,3f10.4,g12.4,2f9.2)') kk,Ncenter, X0, Y0 , Z0,aM,vcrc 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 aM_a = aM_a + aM Vc_a = Vc_a + vcrc 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. Vrms(Ncenter) = vrmss Vcirc(Ncenter) = vcrc Halo(Ncenter) = Ncenter EndIf EndIf Enddo write(*,'(a,i6,a,2f8.1)') ' Read BDM Halos =',Ncenter,' Vc limit =',Vclimit write(*,'(a,g13.4)') & ' Mass =',aM_a/Ncenter, & ' Vcirc =',Vc_a/Ncenter rewind(3) End SUBROUTINE ReadHalos !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadFOF(iFlag) !-------------------------------------------------------------- Character :: Line*120, Txt*8='Coordina', Txt2*8,txt1*4 Nfof = 0; Rvir_a = 0.; aM_a = 0. Do i=1,3 Read(4,'(a)')Line !If(i<10)write(2,'(a)')Line !write(*,'(a)')Line EndDo Do kk=1,10000 read (4,*,iostat=iStat) j, X0, Y0 , Z0, VvX,VvY,VvZ, & k,aM,s1,vrmss,rr,delta,spin !write (*,'(2i9,3f10.4,g12.4,2f9.2)') kk,Ncenter, X0, Y0 , Z0,aM,vcrc,aLambda 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 Nfof = Nfof +1 aM_a = aM_a + aM If(iFlag == 1)Then Xf(Nfof) = X0 Yf(Nfof) = Y0 Zf(Nfof) = Z0 VXf(Nfof) = VvX VYf(Nfof) = VvY VZf(Nfof) = VvZ Amf(Nfof) = aM Vrmf(Nfof) = vrmss Rsph(Nfof) = rr EndIf ! EndIf Enddo write(*,'(a,i10,a,2f8.1)') ' Read FOF Halos =',Nfof write(*,'(a,g13.4)') & ' Mass =',aM_a/Nfof rewind(4) End SUBROUTINE ReadFOF ! ------------------------ 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 MatchFOF use SetHalos Character*120 :: file1,file2,catname,listname,txt*8 !Open( 3,file='CATALOGS/CatshortA.045.DAT',Status='OLD')! short list of halos !Open( 4,file='FOF/cluster_cat_045.l0.20.dat',Status='OLD')! short list of halos Open( 3,file='CATALOGS/CatshortA.416.DAT',Status='OLD')! short list of halos Open( 4,file='FOF/cluster_cat_416.l0.17.dat',Status='OLD')! short list of halos Box =250. dBuffer = 15. Open( 2,file='CatshortMatch.416.DAT')! short list of halos Open( 20,file='CatshortNoMatchFOF.416.DAT')! short list of halos Open( 30,file='CatshortNoMatchBDM.416.DAT')! 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),Vcirc(Ncenter)) ALLOCATE (Vrms(Ncenter),Halo(Ncenter),Match(Ncenter)) write (*,*) 'Allocated Centers' CALL ReadHalos(1,vcLimit) ! read halos CALL ReadFOF(0) ALLOCATE (Xf(Nfof),Yf(Nfof),Zf(Nfof)) ALLOCATE (VXf(Nfof),VYf(Nfof),VZf(Nfof)) ALLOCATE (Amf(Nfof),Vrmf(Nfof),Rsph(Nfof)) CALL ReadFOF(1) 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 Call MatchHalos ! get statistics of halos End PROGRAM MatchFOF