!---------------------------------------------------------------------------------------- ! 3d properies of dark matter halos ! with different isolation conditions. ! Input: Halo Catalog (long list) ! List of distinct halos ! Output: statistics of dm profiles ! Parameters: Isolation parameters: Riso ! ! MODULE setHalos PUBLIC Real, PARAMETER :: & dLog =0.07, & ! bin size in log(r) dMlog =0.15, & dVlog =0.040 Integer, PARAMETER :: Nbin = 40, Nbins =27 Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vcirc Integer, ALLOCATABLE, DIMENSION(:) :: Halo,Lisolate,jIsolate Integer :: Ncenter,Nmain,Ntotal,nnx Real :: Box, Risomax, Vcmin,Vcmax,Vfrac, dBuffer, & Vc_av,aM_av,Rv_av Real, DIMENSION (Nbin) :: Nin,Rad_out,Rad_bin,Mr,Dens, & VcircB,Vrms,Vrad,Vrad2,Jx,Jy,Jz,Jtot 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 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/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 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 !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadHalos(iFlag,Vclimit) !-------------------------------------------------------------- Character :: Line*120, Txt*8='Coordina', Txt2*8,txt1*4 Ncenter = 0; Rvir_a = 0.; aM_a = 0. Vc_a = 0.; aL_a = 0.; x_a = 0. vir_a = 0. Do i=1,11 Read(3,'(a)')Line !If(i<10)write(2,'(a)')Line !write(*,'(a)')Line EndDo Do ! kk=1,1000 read (3,*,iostat=iStat) j,X0, Y0 , Z0, VvX,VvY,VvZ, & aM,rr,vrmss,vcrc,aMunb,vrms_u,VirRat,x_off,aJ,aLambda !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 Ncenter = Ncenter +1 aM_a = aM_a + aMunb Vc_a = Vc_a + vcrc aL_a = aL_a + aLambda x_a = x_a + x_off vir_a = vir_a + VirRat 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 Halo(Ncenter) = j EndIf EndIf Do i =1,Nbins read(3,*,iostat =iStat)n,r1,r2,am1,d1,vc1,vs1,vr1,vr2,ajx1,ajy1,ajz1,aj EndDo Enddo write(*,'(a,i6,a,2f8.1)') ' Read Halos =',Ncenter,' Vc limit =',Vclimit write(*,'(a,g13.4)') & ' Mass =',aM_a/Ncenter, & ' Vcirc =',Vc_a/Ncenter, & ' Lambda =',aL_a/Ncenter, & ' X_offset =', x_a/Ncenter, & ' 2-Epot/Ekin =',vir_a/Ncenter rewind(3) End SUBROUTINE ReadHalos !--------------------------------------------------------- ! Read halos catalog, get final statistics ! SUBROUTINE GetHalos !-------------------------------------------------------------- Character :: Line*120, Txt*8='Coordina', Txt2*8,txt1*4 Logical :: take ! invert table of isolated halos Do i=1,Nmain jIsolate(Halo(i)) = Lisolate(i) EndDo Ncenter = 0; Rvir_a = 0.; aM_a = 0. Vc_a = 0.; aL_a = 0.; x_a = 0. vir_a = 0. Nin = 0; Rad_out =0.; Rad_bin = 0.; Mr =0. Dens = 0.; VcircB = 0.; Vrms = 0.; Vrad = 0.; Vrad2 = 0. Jx = 0.; Jy = 0.; Jz = 0.; Jtot = 0. Do i=1,11 Read(3,'(a)')Line If(i<10)write(2,'(a)')Line EndDo Do ! kk=1,1000 read (3,*,iostat=iStat) j,X0, Y0 , Z0, VvX,VvY,VvZ, & aM,rr,vrmss,vcrc,aMunb,vrms_u,VirRat,x_off,aJ,aLambda !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 > Vlow .and. vcrc < Vup.and.abs(VirRat)<0.5)Then If(vcrc > Vcmin .and. vcrc < Vcmax.and.jIsolate(j)==1)Then take = .true. Ncenter = Ncenter +1 aM_a = aM_a + aMunb Vc_a = Vc_a + vcrc aL_a = aL_a + aLambda x_a = x_a + x_off vir_a = vir_a + VirRat !write (32,'(2i9,3f10.4,g12.4,2f9.2)') kk,Ncenter, X0, Y0 , Z0,aM,vcrc,aLambda Else take = .false. EndIf Do i =1,Nbins read(3,*,iostat =iStat)n,r1,r2,am1,d1,vc1,vs1,vr1,vr2,ajx1,ajy1,ajz1,aj If(take)Then Nin(i) = Nin(i) + n Rad_out(i) = Rad_out(i) + r1 Rad_bin(i) = Rad_bin(i) + r2 Mr(i) = Mr(i) + am1 Dens(i) = Dens(i) + d1 VcircB(i) = VcircB(i) + vc1 Vrms(i) = Vrms(i) + vs1**2 Vrad(i) = Vrad(i) + vr1 Vrad2(i) = Vrad2(i) + vr2 **2 Jtot(i) = Jtot(i) + aj**2 EndIf EndDo Enddo Do j=0,2,2 write(j,'(a,i6,a,2f8.1)') ' Read Halos =',Ncenter,' Vc(km/s) =',Vcmin,Vcmax write(j,'(2(a,g13.4))') & ' Risolate/Rvir=',Risomax, & ' Mass =',aM_a/Ncenter, & ' Vcirc =',Vc_a/Ncenter, & ' Lambda =',aL_a/Ncenter, & ' X_offset =', x_a/Ncenter, & ' 2-Epot/Ekin =',vir_a/Ncenter EndDo write(2,'(T4,a,T12,a,T25,a,T36,a,T44,a,T59,a,T67,a,T76,a,T85,a,T96,a)') & 'Npart/halo','R_out/Rvir','R_av/Rvir','M(0)& write(2,'(f9.1,2f11.4,1p,2g12.4,0p,8f9.3)')(Nin(i))/float(n),Rad_out(i)/n,Rad_bin(i)/n, & Mr(i)/n,Dens(i)/n,VcircB(i)/n, sqrt(Vrms(i)/n), & Vrad(i)/n,sqrt(Vrad2(i)/n),sqrt(Jtot(i)/n) EndDo End SUBROUTINE GetHalos ! ------------------------ 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='CatalogDistinctB.416.DAT',Status='OLD')! 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.0 Vclimit = 80. AEXPN =1.0 ISTEP = 416 write(file1,'(a,i4.4,a,i4.4,a)') 'HalosDistB.',INT((Vcmin+Vcmax)/2.),'.',416,'.DAT' Open( 2,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),Vcirc(Ncenter)) ALLOCATE (Halo(Ncenter),jIsolate(Ncenter)) write (*,*) 'Allocated Centers' CALL ReadHalos(1,vcLimit) ! read halos 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 GetHalos ! get statistics of halos End PROGRAM HaloSatellites