!------------------------------------------------- ! Trace Majot projenitors of halos ! MODULE setHalos PUBLIC Integer, parameter :: Ncats =4, & Nmax = 1400000, & ! max # halos Nbnd = 200 ! max # of particles Integer, DIMENSION(Ncats) :: Nhalos Real, DIMENSION(Ncats) :: Aexpn Real, DIMENSION (Nmax,Ncats) :: Xc,Yc,Zc, & VXc,VYc,Vzc, & Amc,Rmc,Vrms,Vcirc Integer, DIMENSION (Nmax,Ncats) :: nSelect Integer, DIMENSION (Nbnd,Nmax,Ncats) :: nBound Integer, ALLOCATABLE, DIMENSION(:) :: ListMax Contains !--------------------------------------------------------- ! Read halos catalog ! SUBROUTINE ReadHalos(iC) !-------------------------------------------------------------- Character*120 :: Line, Txt*5, Txt2*50 Read(iC+100,'(a)')Line write(*,'(a)') Line read(iC+100,'(a5,f8.3,a)') Txt,aa,Txt2 write(*,'(a5,f8.3,a)')Txt,aa,Txt2 Aexpn(iC) = aa Do i=1,5 Read(iC+100,'(a)')Line write(*,'(a)') Line EndDo Do read(iC+100,*,iostat=iStat) i if(iStat.ne.0)Exit backspace (iC+100) read(iC+100,*,iostat=iStat)j, & Xc(i,iC),Yc(i,iC),Zc(i,iC),VXc(i,iC),VYc(i,iC),VZc(i,iC), & Amc(i,iC),Rmc(i,iC),Vrms(i,iC),Vcirc(i,iC), & nSelect(i,iC) If(iStat.ne.0)Exit If(nSelect(i,iC).ne.0)Then backspace (iC+100) read(iC+100,*,iostat=iStat)j, & Xc(i,iC),Yc(i,iC),Zc(i,iC),VXc(i,iC),VYc(i,iC),VZc(i,iC), & Amc(i,iC),Rmc(i,iC),Vrms(i,iC),Vcirc(i,iC), & nSelect(i,iC),(nBound(k,i,iC),k=1,nSelect(i,iC)) EndIf EndDo 10 Nhalos(iC) =i write(*,*) ' Read Halos=',Nhalos(iC) Return End SUBROUTINE ReadHalos !-------------------------------------------------------------- ! find major progenitor of halo ih of catalog iC ! List of catalogs is an array iCatalog with nCatalogs elements ! Use chain rule: find halo with largest number of ! common particles and use it as ! new center for subsequent search of projenitors SUBROUTINE MajorProjenitor(ih,iC,iCatalog,nCatalogs) !-------------------------------------------------------------- Integer, DIMENSION(nCatalogs) :: iCatalog ihalo = ih ! head of the chain: first halo icat = iC write(1,20) ' Parent=',Aexpn(iC), & Xc(ih,iC),Yc(ih,iC),Zc(ih,iC), & Amc(ih,iC),Vcirc(ih,iC) 10 format(2i6,f8.3,(T40,0p3f9.4,1p5g12.3)) 20 format(12x,a,f8.3,(T40,0p3f9.4,1p5g12.3)) Do i =1,nCatalogs If(i.ne.iC)Then CALL IndPart(ihalo,icat,iCatalog(i),jhalo,nn,0) if(jhalo==0)write(1,'(20x,6i7)') ihalo,icat,iCatalog(i),jhalo,nn If(jhalo.ne.0)Then ihalo = jhalo icat = iCatalog(i) EndIf EndIf EndDo Return End SUBROUTINE MajorProjenitor !--------------------------------------------------------- ! compare particles of halo ih in Catalog iC with catalog iG ! jh = identified halo in iG ! iCross = total number of identical particles SUBROUTINE IndPart(ih,iC,iG,jh,iCross,iFlag) !-------------------------------------------------------------- Nc = nSelect(ih,iC) If(Nc==0)Then write(*,*) ' No selected particles in halo=',ih,iC jh =0 iCross =0 return EndIf ALLOCATE (ListMax(Nhalos(iG))) ListMax = 0 nmx = 0 jmx = 0 !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (j,i,k,Ng,ii,kk) Do j=1,Nhalos(iG) Ng =nSelect(j,iG) If(Ng.ne.0)Then Do i=1,Nc ii =nBound(i,ih,iC) ! particle in first halo Do k=1,Ng kk =nBound(k,j,iG) ! particle in second halo if(kk==ii)ListMax(j) =ListMax(j) +1 EndDo EndDo EndIf EndDo Do j=1,Nhalos(iG) If(ListMax(j).ne.0)Then if(ListMax(j)>nmx)Then nmx =ListMax(j) ; jmx =j EndIf ! write(*,10) ih,j,icr,Xc(ih,iC),Yc(ih,iC),Zc(ih,iC),Amc(ih,iC), & ! Xc(j,iG),Yc(j,iG),Zc(j,iG),Amc(j,iG) EndIf EndDo If(nmx.eq.0)Then ! write(*,10) ih,0,0,Xc(ih,iC),Yc(ih,iC),Zc(ih,iC),Amc(ih,iC) jh =0; iCross =0 Else jh =jmx ; iCross =nmx If(iFlag==0)write(1,10) jmx,nmx,Aexpn(iG), & Xc(jmx,iG),Yc(jmx,iG),Zc(jmx,iG), & Amc(jmx,iG),Vcirc(jmx,iG), & VXc(jmx,iG),VYc(jmx,iG),VZc(jmx,iG) EndIf 10 format(2i6,f8.3,(T40,0p3f9.4,1p2g12.3,0p3f9.2)) DEALLOCATE (ListMax) End SUBROUTINE IndPart end module SetHalos !------------------------------------------------- ! Find halos at different moments, which ! have common particles PROGRAM HaloCrossIdentity use SetHalos Character*120 :: file1,listname,txt*8 Integer, DIMENSION(Ncats) :: iCatalog Character*10, DIMENSION(Ncats) :: iNames ! ------------------------------------------ Read data Open(1,file='HaloCrossList.dat') Do iC=1,Ncats write(*,'(a,$)')' Enter snapshot=' read(*,'(a)')iNames(iC) endDo !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (iC,listname) Do iC=1,Ncats listname='Halolist.'//TRIM(iNames(iC))//'.DAT' Open( iC+100,file=listname,Status='UNKNOWN')! short list of halos Call ReadHalos(iC) ! count halos iCatalog(iC) = iC EndDo Ncatalogs = Ncats nc =0 nm =0 ! Do ih=1,Nhalos(1) ! identify halo in the second catalog ! CALL IndPart(ih,1,2,j,nn,1) ! If(Amc(ih,1)>1.e8)Then ! nm = nm+1 ! if(nn==0)nc =nc+1 ! EndIf ! EndDo !write (*,*) nm,nc Call MajorProjenitor(649733,1,iCatalog,nCatalogs) Call MajorProjenitor(653246,1,iCatalog,nCatalogs) Call MajorProjenitor(414103,1,iCatalog,nCatalogs) Call MajorProjenitor(410918,1,iCatalog,nCatalogs) stop end Program HaloCrossIdentity