C ______________________________ Selects a halo from a large Catalog.DAT C C C Scales internal PM coordinates and velocities C to coordinates in Mpc/h and km/s C In PM the coordinates are in the range 1 - (NGRID+1) C velocities are P = a_expansion*V_pec/(x_0H_0) C where x_0 = comoving cell_size=Box/Ngrid C H_0 = Hubble at z=0 C C NROW = number of particles in 1D C NGRID= number of cells in 1D INCLUDE 'PMselect.h' REAL INPUT Character*50 FileASCII Logical Inside C................................................................... C Read data and open files WRITE (*,'(A,$)') 'Enter Name for output file = ' READ (*,'(A)') FileASCII OPEN(17,FILE=FileASCII,STATUS='UNKNOWN') WRITE (*,'(A,$)') 'Enter Name of input file = ' READ (*,'(A)') FileASCII OPEN(20,FILE=FileASCII,STATUS='UNKNOWN') WRITE (*,'(A,$)') 'Enter Name of catalog with', & ' corrected concentrations file = ' READ (*,'(A)') FileASCII OPEN(21,FILE=FileASCII,STATUS='UNKNOWN') Box =INPUT(' Enter box size in comoving Mpc/h =') dBox = 10. ! expand box by dBox Mpc/h write (*,'(A,$)')' Enter selection Limits: ', + 'Mass dM X Y Z dR C dC =>' read (*,*) aMass, dM, X0, Y0, Z0, dR, Cc, dCc CALL ReadCatalog(Nh,Nhex) CALL ExpandSample(Nh,Nhex,dBox,Box) write(17,30)aMass-dM,aMass+dM,Cc-dCC,Cc+dCC,aMiso,Riso c CALL PrintSelected(Nh,Nhex) CALL AverageSelected(Nh,Nhex) 30 format(' Mass range=',2g11.4,' Concentrations=',2g11.3, & /' Vcirc_neighbor/Vcirc < ',f8.2,' Rneighbor/Rvir> ',f8.2) END c-------------------------------------------------------------------- SUBROUTINE AverageSelected(Nh,Nhex) c-------------------------------------------------------------------- INCLUDE 'PMselect.h' Logical inside,isolated COMMON /AVER/ Rapr(Nrmax),Napr(Nrmax), + Vacpr(Nrmax),RHOapr(Nrmax),Varpr(Nrmax), + bMa(Nrmax),Dapr(Nrmax),Vadpr(Nrmax), & write (*,*) ' Concentration limits=',Cc-dCc,Cc+dCc nn = 0 aMav =0. rav =0. Do i=1,Nh amm = aM(i) Cnc = Conc(i) x = Xp(i) y = Yp(i) z = Zp(i) inside = (amm.gt.aMass-dM.and. amm.lt.aMass+dM ) inside = (Cnc.gt.Cc-dCc.and.Cnc.lt.Cc+dCc.and.inside) inside = (x.gt.X0-dR.and. x.lt.X0+dR .and.inside) inside = (y.gt.Y0-dR.and. y.lt.Y0+dR .and.inside) inside = (z.gt.Z0-dR.and. z.lt.Z0+dR .and.inside) inside = (errFit(i).lt.errLimit.and.inside) If(inside)Then CALL CheckIsolated(i,Nh,Nhex,isolated) c isolated =.true. If(isolated)Then write(*,50) Xp(i),Yp(i),Zp(i), & VXp(i),VYp(i),VZp(i), & aM(i),Radh(i),Vrms(i),Vcirc(i), & Npart(i),Rmax(i),Conc(i),Nlines(i) nn = nn +1 aMav = aMav + aM(i) rav = rav + Radh(i) If(nn.eq.1)Then ifirst = i Do j=1,Nlines(i) Rapr(j) = Rpr(j,i)/ Radh(i) Vacpr(j) = Vcpr(j,i) RHOapr(j)= RHOpr(j,i) bMa(j) = RHOpr(j,i)**2 Varpr(j) = Vrpr(j,i) EndDo else ilast = i Do j=1,Nlines(i) Rapr(j) = Rapr(j) +Rpr(j,i)/ Radh(i) Vacpr(j) = Vacpr(j) +Vcpr(j,i) RHOapr(j)= RHOapr(j)+RHOpr(j,i) bMa(j) = bMa(j) +RHOpr(j,i)**2 Varpr(j) = Varpr(j) +Vrpr(j,i) EndDo EndIf EndIf EndIf EndDo If(nn.eq.0)Then write (*,*) ' no halos were selected' return EndIf write (*,*) nn,' halos were selected' Do j=1,Nlines(ifirst) Rapr(j) = Rapr(j)/nn Vacpr(j) = Vacpr(j)/nn RHOapr(j)= RHOapr(j)/nn bMa(j) = sqrt(bMa(j)/nn -RHOapr(j)**2) Varpr(j) = Varpr(j)/nn EndDo aMav = aMav/nn rav = rav/nn write(17,'(" Nselected=",i6," Mass=",g11.3," Rvir=",g11.3)') & nn,aMav,rav write(17,50) Xp(ilast),Yp(ilast),Zp(ilast), & VXp(ilast),VYp(ilast),VZp(ilast), & aM(ilast),Radh(ilast),Vrms(ilast),Vcirc(ilast), & Npart(ilast),Rmax(ilast),Conc(ilast),Nlines(ilast) Do j=1,Nlines(ilast) write(17,60) Rapr(j),Npr(j,ilast), & bM(j,ilast),Dpr(j,ilast),Vdpr(j,ilast), & Vacpr(j),RHOapr(j)*Dpr(1,1)/RHOpr(1,1), & Npr(j,ilast),0,0,Varpr(j),bMa(j)*Dpr(1,1)/RHOpr(1,1) EndDo 40 Format(2i7,3g11.3,g12.4,f9.2) 50 Format(F8.3,2F8.3,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4) 60 format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,2g11.3) Return End c-------------------------------------------------------------------- SUBROUTINE PrintSelected(Nh,Nhex) c-------------------------------------------------------------------- INCLUDE 'PMselect.h' Logical inside,isolated nn = 0 Do i=1,Nh amm = aM(i) Cnc = Conc(i) x = Xp(i) y = Yp(i) z = Zp(i) inside = (amm.gt.aMass-dM.and. amm.lt.aMass+dM ) inside = (Cnc.gt.Cc-dCc.and.Cnc.lt.Cc+dCc.and.inside) inside = (x.gt.X0-dR.and. x.lt.X0+dR .and.inside) inside = (y.gt.Y0-dR.and. y.lt.Y0+dR .and.inside) inside = (z.gt.Z0-dR.and. z.lt.Z0+dR .and.inside) inside = (errFit(i).lt.errLimit.and.inside) If(inside)Then CALL CheckIsolated(i,Nh,Nhex,isolated) c isolated =.true. If(isolated)Then nn = nn +1 write(*,40) nn,i,Conc(i),aM(i),Vcirc(i) write(17,50) Xp(i),Yp(i),Zp(i), & VXp(i),VYp(i),VZp(i), & aM(i),Radh(i),Vrms(i),Vcirc(i), & Npart(i),Rmax(i),Conc(i),Nlines(i) c Do j=1,Nlines(i) c write(17,60) Rpr(j,i)/Radh(i),Npr(j,i), c & bM(j,i),Dpr(j,i),Vdpr(j,i), c & Vcpr(j,i),RHOpr(j,i)*Dpr(1,i)/RHOpr(1,i), c & Npr(j,i),0,0,Vrpr(j,i) c EndDo EndIf EndIf EndDo write (*,*) nn,' halos were selected' 40 Format(2i7,3g11.3,g12.4,f9.2) 50 Format(F8.3,2F8.3,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4) 60 format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,g11.3) Return End c-------------------------------------------------------------------- SUBROUTINE CheckIsolated(i,Nh,Nhex,isolated) c-------------------------------------------------------------------- INCLUDE 'PMselect.h' Logical isolated isolated = .true. amm = Vcirc(i) ! aM(i) c amm = aM(i) Rhalo = Radh(i) x = Xp(i) y = Yp(i) z = Zp(i) Do j=1,Nhex If(i.ne.j.and.amm.lt.Vcirc(j)*aMiso)Then dd = (Xp(j)-x)**2+(Yp(j)-y)**2+(Zp(j)-z)**2 diso = Riso*max(Rhalo,Radh(j))/1.e3 ! put it in Mpc c write (*,'(2i7,10g11.3)')i,j,sqrt(dd),diso,Rhalo, c & Radh(j),Vcirc(i),Vcirc(j) If(dd.lt.diso**2)Then isolated = .false. return EndIf EndIf EndDo If(Vrms(i).gt.1.4*Vcirc(i))write (*,*)' Unbound halo:', & i,Vrms(i),Vcirc(i) Return End c-------------------------------------------------------------------- SUBROUTINE ReadCatalog(Nh) c-------------------------------------------------------------------- INCLUDE 'PMselect.h' Character*120 Line Lines = 27 Do i=1,Lines read(20,'(A)')Line read(21,'(A)')Line write(*,'(A)')Line write(17,'(A)')Line EndDo Nh = 1 10 read(20,*,end=100,err=100) Xp(Nh),Yp(Nh),Zp(Nh), & VXp(Nh),VYp(Nh),VZp(Nh), & aM(Nh),Radh(Nh),Vrms(Nh),Vcirc(Nh), & Npart(Nh),Rmax(Nh),Conc(Nh),Nlines(Nh) read(21,*,end=100,err=100) Xp0,Yp0,Zp0, & VXp0,VYp0,VZp0, & aM0,Radh0,Vrms0,Vcirc0, & Npart0,Rmax0,Conc0,Nlines0,errFit(Nh) If(Nh.lt.10.or.Nh/5000*5000.eq.Nh) & write (*,'(" header=",i7,g11.3,i6,3g12.4)') & Nh,aM(Nh),Nlines(Nh),Xp(Nh),Yp(Nh),Zp(Nh) If(Nh.lt.10.or.Nh/5000*5000.eq.Nh) & write (*,'(" header=",i7,g11.3,i6,3g12.4)') & Nh,aM0,Nlines0,Xp0,Yp0,Zp0 If(abs(Xp(Nh)-Xp0).gt.1.e-2.or.abs(Yp(Nh)-Yp0).gt.1.e-2 & .or.abs(Zp(Nh)-Zp0).gt.1.e-2)stop 'Mismatch between catalogs' If(abs(aM(Nh)-aM0).gt.1.e-2*aM0.or.Nlines(Nh).ne.Nlines0) & stop 'Mismatch between catalogs' Conc(Nh) =Conc0 If(Nlines(Nh).gt.Nrmax) STOP ' wrong number of profile lines' If(Nlines(Nh).le.0) STOP ' wrong number of profile lines' Do i=1,Nlines(Nh) read (20,*) Rpr(i,Nh),Npr(i,Nh), & bM(i,Nh),Dpr(i,Nh),Vdpr(i,Nh), & Vcpr(i,Nh),RHOpr(i,Nh), & nsm,nmed,nlg,Vrpr(i,Nh) EndDo if(Nh.ge.Np)goto100 Nh =Nh +1 goto 10 100 Nh =Nh -1 write (*,'(" Number of halos read=",T30,i8)') Nh Return End c-------------------------------------------------------------------- c make sample larger by replicating c some halos using periodical c boundary conditions c dR is the width of the buffer in Mpc/h c Next = new number of halos SUBROUTINE ExpandSample(Nh,Nhex,dBox,Box) c-------------------------------------------------------------------- INCLUDE 'PMselect.h' Logical inx,iny,inz Xleft =-dBox ! new boundaries for the catalog Xright =Box+dBox Yleft =-dBox Yright =Box+dBox Zleft =-dBox Zright =Box+dBox write (*,*) Xleft ,Xright Next = Nh Do ip =1,Nh ! go over all halos x =Xp(ip) y =Yp(ip) z =Zp(ip) Do k=-1,1 ! go over all 27 periodical images dz = k*Box inz =(z+dz.gt.Zleft.and.z+dz.lt.Zright) Do j=-1,1 dy = j*Box iny =(y+dy.gt.Yleft.and.y+dy.lt.Yright) Do i=-1,1 dx =i*Box If(i**2+j**2+k**2.ne.0)Then inx =(x+dx.gt.Xleft.and.x+dx.lt.Xright) If(inx.and.iny.and.inz)Then ! add this halo to the list Next =Next +1 if(Next.gt.Np) STOP 'too many halos' Xp(Next) =x+dx Yp(Next) =y+dy Zp(Next) =z+dz VXp(Next) =VXp(ip) VYp(Next) =VYp(ip) VZp(Next) =VZp(ip) aM(Next) =aM(ip) Radh(Next) =Radh(ip) Vrms(Next) =Vrms(ip) Vcirc(Next)=Vcirc(ip) Conc(Next) =Conc(ip) EndIf ! end inside box EndIf ! end exclude central box EndDo ! end i EndDo ! end j EndDo ! end k EndDo ! end ip write (*,*) ' Nhalos in extended sample=',Next Nhex =Next RETURN END C---------------------------------- Read in variables REAL FUNCTION INPUT(text) C------------------------------------------------ Character text*(*) write (*,'(A,$)')text read (*,*) x INPUT =x Return End