!--------------------------------------------------------------------------------- ! Read HaloTrace files ! Extract information on halo evolution !-------------------------------------------------------------------------------- MODULE setHalos PUBLIC Integer, parameter :: Ncats =178 ! Maximum number of snapshots to analyze Real, parameter :: Box=250., & Omega0 = 0.27, & ! Omega_matter hubble = 0.70 Real, ALLOCATABLE, DIMENSION (:) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc,VcircMax,AmcMax Integer, ALLOCATABLE, DIMENSION (:) :: iHalo,iSub,iSub0,iPoint Real, DIMENSION(0:100) :: aTrack,aTrSub0,aTrSub1 Integer :: Ncenter,Nall,Ndistinct Real :: da =0.015 TYPE :: hal Real*4 :: aexpn ! expansion parameter Integer*4 :: InCat ! pointer to halo in Catshort Integer*4 :: Parent ! if not 0, pointer to parent halo in Catshort Real :: x,y,z,vx,vy,vz ! coordinates (Mpch units), peculiar velocities(km/s) Real :: aM,Rm,circV ! Mass(Msunh), Vir radius(kpch), max circ Velocity End TYPE hal TYPE :: haloEvolution Integer :: N ! number of moments for which halo was traced Type(hal) :: h(Ncats) End TYPE haloEvolution Type(haloEvolution) :: halo Contains !--------------------------------------------------------- ! Read 'iC' halo catalog ! SUBROUTINE ReadHalos(iC) !-------------------------------------------------------------- Character*120 :: listname, Txt*5, Txt2*50 Type(hal) :: halo2 Real :: buff(Ncats) ! Open file write(listname,'(a,i3.3,a)') 'HaloTree/HaloTrace.',iC,'.dat' Open( iC+100,file=TRIM(listname),Status='UNKNOWN',form='unformatted') read(iC+100) ! first line is empty nHalos = 0 Do ! Loop over all halos in the file read(iC+100,iostat=iStat) ih If(iStat.ne.0)exit Line =0 Do read(iC+100,iostat=iStat) halo2 ! there is a bug in reading these data backspace(iC+100) ! read it twice to get around the bug read(iC+100,iostat=iStat) halo2%aexpn ! read aexpn again If(iStat.ne.0)exit Line =Line +1 halo%N = ih halo%h(Line) = halo2 !write(*,'(i9,i8,f9.4,2i9,3f9.4,3f9.2,1pg12.4,0p2f8.2)')ih,nHalos, halo2 EndDo nHalos = nHalos +1 ! smooth mass history buff(1) = halo%h(1)%aM buff(Line) = (halo%h(Line)%aM+halo%h(Line-1)%aM)/2. Do i=2,Line-1 buff(i) = (halo%h(i-1)%aM+halo%h(i)%aM+halo%h(i+1)%aM)/3. EndDo Do i=1,Line halo%h(i)%aM = buff(i) EndDo ! smooth Vcirc history buff(1) = (halo%h(1)%circV+halo%h(2)%circV)/2. buff(Line) = (halo%h(Line)%circV+halo%h(Line-1)%circV)/2. Do i=2,Line-1 buff(i) = (halo%h(i-1)%circV+halo%h(i)%circV+halo%h(i+1)%circV)/3. EndDo Do i=1,Line halo%h(i)%circV = buff(i) EndDo If(halo%h(1)%aexpn >= 1.)Then ! halo exists at z=0 !write(*,*) ' halo =',nHalos, ' a=', halo%h(1)%aexpn !write(*,'(10i8)') (halo%h(j)%InCat,j=1,Line) iP = iPoint(halo%h(1)%InCat) ! position in selected halo list If(iP >0)Then ! selected halo If(halo%h(1)%Parent /= 0)iSub0(iP) = 1 ! not distinct at z= 0 Nn = 0 aMax = 0. ; vv =0. Do i = 1,Line If(halo%h(i)%Parent /= 0) Nn =Nn +1 aMax = max(halo%h(i)%aM,aMax) vv = max(halo%h(i)%circV,vv) EndDo If(Nn > 1)iSub(iP) = Nn VcircMax(iP) = vv AmcMax(iP) = aMax ! If(iSub0(iP) == 0 ) Then ! halo is distinct at z=0 Ndistinct = Ndistinct + 1 indx = MIN(INT(halo%h(Line)%aexpn/da),100) aTrack(indx) = aTrack(indx) +1 ! last time it was detected If(Nn > 1)aTrSub0(indx) = aTrSub0(indx) +1 ! was a subhalo If(aMax > 1.2*halo%h(1)%aM .and. Nn <= 1)aTrSub1(indx) = aTrSub1(indx) +1 ! mass declined ! EndIf EndIf EndIf EndDo write(*,'(a,i9,a,i4)') ' Read Halos=',nHalos,' File =',iC close(100+iC) Return End SUBROUTINE ReadHalos !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadHalosCat(iFlag,vcLimit) !-------------------------------------------------------------- Character :: Line*80, Txt*6='Number', Txt2*8,txt1*4 Ncenter =0 Rdmax =0. Do i=1,2 Read(3,'(a)')Line If(iFlag==0)write(2,'(a)') Line EndDo Read(3,'(a4,i4)')txt1,ii If(iFlag==0)Write(2,*)txt1,ii indx =0 Do read(3,*,err=20)Txt2 !write(2,*)TRIM(Txt2) If(TRIM(Txt2)==Txt)Then If(indx==1)Exit indx =1 EndIf EndDo Ind = 0 Do ! i=1, 50 read (3,45,iostat=iStat) X0, Y0 , Z0, VvX,VvY,VvZ, & aM,rr,vrmss,vcrc,Npart,Rmax,Conc !,Nb ! 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 Ind = Ind +1 If(vcrc> 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) iHalo(Ncenter) = Ind ! pointer from selected to original iPoint(Ind) = Ncenter ! pointer from original to selected if(iPoint(ind)==0)write(*,*) ' error: ind=',ind,' Ncenter =',Ncenter EndIf EndIf ! test mass Enddo 40 Format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,g11.3) Nall = Ind If(iFlag==0)write(*,*) 'Read Halos=',Ncenter,' Total =',Nall If(iFlag==1)write(*,'(a,i8,a,f8.3,a)') ' Read Halos=',Ncenter, & ' Max Radius=',Rdmax,'Mpch' rewind(3) Return 20 write(*,*) ' error reading halo catalog' stop End SUBROUTINE ReadHalosCat End module setHalos !------------------------------------------------- ! ! PROGRAM HaloCrossIdentity use SetHalos Character*120 :: file1,listname,txt*8 Call CPU_time(t0); call system_clock(iT0,iTc,ib) write(*,*) write(*,'(a,$)')' Enter limit on circular velocity and output file number (0-9)=' read(*,*)vcLimit,iFile Open( 3,file='CATALOGS/CatshortA.416.DAT',Status='OLD')! short list of halos Open( 2,file='CatshortD.416.DAT') ! updated short list of halos write(file1,'(a,i1,a)') 'StatTrack.',iFile,'.DAT' Open( 10,file=file1) ! Table of tracking statistics CALL ReadHalosCat(0,vcLimit) ! count halos write(2,'(a,T20,a,T38,a,5(2x,a))')' Halo','Coords(Mpch)', & 'Mass(z=0)','MassMax','Vcirc(z=0)','Vcirc(max)','Sub(z=0)','SubEver' ALLOCATE (Xc(Ncenter),Yc(Ncenter),Zc(Ncenter)) ALLOCATE (VXc(Ncenter),VYc(Ncenter),VZc(Ncenter)) ALLOCATE (Amc(Ncenter),Rmc(Ncenter),Vrms(Ncenter),Vcirc(Ncenter)) ALLOCATE (AmcMax(Ncenter),VcircMax(Ncenter)) ALLOCATE (iHalo(Ncenter),iSub(Ncenter),iSub0(Ncenter)) ALLOCATE(iPoint(Nall)) iHalo = 0 ; iSub =0 ; iSub0 =0 ; AmcMax =0. ; VcircMax = 0. aTrack = 0.; aTrSub0 = 0.; aTrSub1 = 0. CALL ReadHalosCat(1,vcLimit) ! read z=0 halos AmcMax = Amc VcircMax = Vcirc Ndistinct = 0 Call CPU_time(t1); call system_clock(iT2,iTc,ib) ; write(*,*) ' time for Init =',t1-t0,float(iT2-iT0)/iTc ! --------------------------------------------- Main loop Do iC =1, 100 Call ReadHalos(iC) Call CPU_time(t1); call system_clock(iT2,iTc,ib) ; write(*,*) ' time =',t1-t0,float(iT2-iT0)/iTc EndDo write(*,*) ' Ndistinct =',Ndistinct write(10,'(10("-"),a/,10x,a,f6.2,a,i10)') ' Statistics of halo tracking','Vclimit=',Vclimit,' km/s Ndistinct=',Ndistinct write(10,'(3x,a)')' Aexpn Z dNtracked dN(were_subhalos) dN(were_larger) Ntracked N(were_subhalos) N(were_larger)' g =FLOAT(Ndistinct) iTa =0 ; iT0 =0; iT1 =0 do i=0,100 aa =(i+0.5)*da iTa = iTa +aTrack(i) iT0 = iT0 +aTrSub0(i) iT1 = iT1 +aTrSub1(i) if(aa<=1. .and. iTa>0)write(10,'(2f8.3,1p,6g14.3)') aa,1./aa-1,aTrack(i)/g,aTrSub0(i)/g,aTrSub1(i)/g, & iTa/g, iT0/g, iT1/g Enddo !Do i=1, Ncenter ! write(2,'(i9,3f9.4,1p,2g11.4,0p,2f9.2,2i8)')iHalo(i),Xc(i),Yc(i),Zc(i), & ! Amc(i),AmcMax(i),Vcirc(i),VcircMax(i),iSub0(i),iSub(i) !EndDo stop end Program HaloCrossIdentity