!---------------------------------------------------------------------------------------- ! Find properties of FOF halos ! ! Input: Halo Catalog (short list) ! ! Output: Distribution functions ! Parameters: ! ! ! MODULE setHalos PUBLIC Real, PARAMETER :: & dLog = 0.1, & ! bin size in log(r) dMlog =0.15 Integer, ALLOCATABLE :: LDistinct(:),Lst(:),Label(:,:,:) Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms Integer :: Ncenter,Nmain Real :: Box Contains !-------------------------------------------------------------- ! 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 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 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(2,*) ' Mass(Msunh) N(>M) Mass dN MdN/dM' Do i=istart,ifinish ! print results If(mHalos(i)>0)Then aMm =10.**(i*dMlog) dM = 10.**((i+1)*dMlog)-aMm write(2,'(3x,1P,G11.4,0P,i9,2x,1P,g11.4,0P,i9,2x,1P,g11.4)')aMm, mHalos(i),massbin(i),iHalos(i), & massbin(i)*iHalos(i)/dM/Box**3 EndIf EndDo End SUBROUTINE MassDistinct !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadHalos(iFlag) !-------------------------------------------------------------- Character :: Line*80, Txt*6='Number', Txt3*8='Coordina',Txt2*8,txt1*4 Nmain =0 Rdmax =0. Do i=1,3 Read(3,'(a)')Line If(iFlag==0)write(*,'(a)') Line If(iFlag==1)write(2,'(a)') Line EndDo Do read (3,*,iostat=iStat) i,X0, Y0 , Z0, VvX,VvY,VvZ, & nn,aM,sig,vrmss,rr,delt,alambda !read (3,*,iostat=iStat) ii,X0, Y0 , Z0, & ! aM0,aM,vnow,vcrc !rr = aM0**0.33333/4.82e1 ! virial radius in kpch for Om=0.27 delta=360.8 !write (*,'(i9,3f10.4,2g12.3,3f9.2)') Nmain, X0, Y0 , Z0,aM !,aM,rr,vnow,vcrc If(iStat.ne.0)Exit Nmain = Nmain +1 If(iFlag==1)Then if(Nmain>Ncenter)Stop ' Not enough space. Increase Ncenter' Xc(Nmain) = X0 Yc(Nmain) = Y0 Zc(Nmain) = Z0 VXc(Nmain) = VvX VYc(Nmain) = VvY VZc(Nmain) = VvZ Amc(Nmain) = aM Rmc(Nmain) = rr Vrms(Nmain) = vrmss EndIf 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=',Nmain If(iFlag==1)write(*,'(a,i8,a,f8.3,a)') ' Read Halos=',Nmain, & ' 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 end module SetHalos !-------------------------------------------------------------------------------------- ! ! PROGRAM HaloFOF use SetHalos Character*120 :: file1,file2,catname,listname,txt*8 Open( 3,file='cluster_cat_045.l0.20.dat',Status='OLD')! short list of halos Open( 2,file='Halo.massFOF_045.l0.020.dat') ! output Box =250. CALL ReadHalos(0) ! count halos Ncenter = Nmain write(*,*) ' Maximum of halos Ncenter =',Ncenter ALLOCATE (Xc(Ncenter),Yc(Ncenter),Zc(Ncenter)) ALLOCATE (VXc(Ncenter),VYc(Ncenter),VZc(Ncenter)) ALLOCATE (Amc(Ncenter),Rmc(Ncenter),Vrms(Ncenter)) !$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. ; Vrms(i) =0. EndDo write (*,*) 'Allocated Centers' CALL ReadHalos(1) ! read halos Close(3) Call MassDistinct ! make list of satellites write (*,*) 'Found Mass distribution' End PROGRAM HaloFOF Subroutine R_mrgref (XVALT, IRNGT,N) ! Ranks array XVALT into index array IRNGT, using merge-sort ! __________________________________________________________ ! This version is not optimized for performance, and is thus ! not as difficult to read as some other ones. ! Michel Olagnon - April 2000 ! __________________________________________________________ ! _________________________________________________________ Real, Dimension (N), Intent (In) :: XVALT Integer, Dimension (N), Intent (Out) :: IRNGT ! __________________________________________________________ ! Integer, Dimension (:), Allocatable :: JWRKT Integer :: LMTNA, LMTNC Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB ! NVAL = Min (SIZE(XVALT), SIZE(IRNGT)) If (NVAL <= 0) Then Return End If ! ! Fill-in the index array, creating ordered couples ! Do IIND = 2, NVAL, 2 If (XVALT(IIND-1) < XVALT(IIND)) Then IRNGT (IIND-1) = IIND - 1 IRNGT (IIND) = IIND Else IRNGT (IIND-1) = IIND IRNGT (IIND) = IIND - 1 End If End Do If (Modulo (NVAL, 2) /= 0) Then IRNGT (NVAL) = NVAL End If ! ! We will now have ordered subsets A - B - A - B - ... ! and merge A and B couples into C - C - ... ! Allocate (JWRKT(1:NVAL)) LMTNC = 2 LMTNA = 2 ! ! Iteration. Each time, the length of the ordered subsets ! is doubled. ! Do If (LMTNA >= NVAL) Exit IWRKF = 0 LMTNC = 2 * LMTNC IWRK = 0 ! ! Loop on merges of A and B into C ! Do IINDA = IWRKF IWRKD = IWRKF + 1 IWRKF = IINDA + LMTNC JINDA = IINDA + LMTNA If (IWRKF >= NVAL) Then If (JINDA >= NVAL) Exit IWRKF = NVAL End If IINDB = JINDA ! ! Shortcut for the case when the max of A is smaller ! than the min of B (no need to do anything) ! If (XVALT(IRNGT(JINDA)) <= XVALT(IRNGT(JINDA+1))) Then IWRK = IWRKF Cycle End If ! ! One steps in the C subset, that we create in the final rank array ! Do If (IWRK >= IWRKF) Then ! ! Make a copy of the rank array for next iteration ! IRNGT (IWRKD:IWRKF) = JWRKT (IWRKD:IWRKF) Exit End If ! IWRK = IWRK + 1 ! ! We still have unprocessed values in both A and B ! If (IINDA < JINDA) Then If (IINDB < IWRKF) Then If (XVALT(IRNGT(IINDA+1)) > XVALT(IRNGT(IINDB+1))) & & Then IINDB = IINDB + 1 JWRKT (IWRK) = IRNGT (IINDB) Else IINDA = IINDA + 1 JWRKT (IWRK) = IRNGT (IINDA) End If Else ! ! Only A still with unprocessed values ! IINDA = IINDA + 1 JWRKT (IWRK) = IRNGT (IINDA) End If Else ! ! Only B still with unprocessed values ! IRNGT (IWRKD:IINDB) = JWRKT (IWRKD:IINDB) IWRK = IWRKF Exit End If ! End Do End Do ! ! The Cs become As and Bs ! LMTNA = 2 * LMTNA End Do ! ! Clean up ! Deallocate (JWRKT) Return ! End Subroutine R_mrgref