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 long output file = ' READ (*,'(A)') FileASCII OPEN(17,FILE=FileASCII,STATUS='UNKNOWN') WRITE (*,'(A,$)') 'Enter Name for short output file = ' READ (*,'(A)') FileASCII OPEN(18,FILE=FileASCII,STATUS='UNKNOWN') WRITE (*,'(A,$)') 'Enter Name of input long file = ' READ (*,'(A)') FileASCII OPEN(20,FILE=FileASCII,STATUS='UNKNOWN') WRITE (*,'(A,$)') 'Enter Name of input short file = ' READ (*,'(A)') FileASCII OPEN(21,FILE=FileASCII,STATUS='UNKNOWN') Box =INPUT(' Enter box size in comoving Mpc/h =') dBox = 5. ! 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) c CALL PrintSelected(Nh,Nhex) CALL AverageSelected(Nh,Nhex) 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 write (*,*) ' X position limits=',X0-dR,X0+dR write (*,*) ' Y position limits=',Y0-dR,Y0+dR write (*,*) ' Z position limits=',Z0-dR,Z0+dR nn = 0 aMav =0. rav =0. Nlg =0 ! number of large particles 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) 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) Nlg = Nlg + nlarge(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) Vadpr(j) = Vdpr(j,i)**2 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) Vadpr(j) = Vadpr(j) +Vdpr(j,i)**2 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 Vadpr(j) = sqrt(Vadpr(j)/nn) EndDo aMav = aMav/nn rav = rav/nn Nlg = Nlg/nn write(17,'(" Nselected=",i6," Mass=",g11.3," Rvir=", & f8.1," NlargePart=",i7," Conc=",2f6.1," Visol=",f6.2, & " Risolate=",f6.2)') & nn,aMav,rav,Nlg,Cc-dCc,Cc+dCc,aMiso,Riso write(17,50) Xp(ifirst),Yp(ifirst),Zp(ifirst), & VXp(ifirst),VYp(ifirst),VZp(ifirst), & aM(ifirst),Radh(ifirst),Vrms(ifirst),Vcirc(ifirst), & Npart(ifirst),Rmax(ifirst),Conc(ifirst),Nlines(ifirst) 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)+Npr(j,ifirst)+1)/2, & 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) x = Xp(i) y = Yp(i) z = Zp(i) inside = (amm.gt.aMass-dM.and. amm.lt.aMass+dM ) 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) If(inside)Then CALL CheckIsolated(i,Nh,Nhex,isolated) c isolated =.true. If(isolated)Then nn = nn +1 write(*,40) nn,i,Xp(i),Yp(i),Zp(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) write(18,55) nn,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) Do j=1,Nlines(i) write(17,60) Rpr(j,i),Npr(j,i), & bM(j,i),Dpr(j,i),Vdpr(j,i), & Vcpr(j,i),RHOpr(j,i), & Npr(j,i),0,0,Vrpr(j,i) EndDo EndIf EndIf EndDo 40 Format(2i7,3g12.5,g12.4,f9.2) 50 Format(F8.3,2F8.3,3F8.1,g11.3,f8.2,2f7.1,I9,f8.1,f8.2,i4) 55 Format(i7,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 If(nlarge(i).gt.Npart(i)/3)Then isolated = .false. write (*,*) ' Too many nLarge=',nlarge(i), Npart(i) return EndIf isolated = .true. c return !!!!!!!!!!!!!!!!! amm = aM(i) !Vcirc(i) ! Rhalo = Radh(i) x = Xp(i) y = Yp(i) z = Zp(i) Do j=1,Nhex If(i.ne.j.and.amm.lt.aM(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. c write (*,'(2i7,10g11.3)')i,j,sqrt(dd),diso,Rhalo, c & Radh(j),Vcirc(i),Vcirc(j) c write (*,'(14x,3f9.3," M=",g11.3)')x,y,z,aM(i) c write (*,'(14x,3f9.3," M=",g11.3)')Xp(j),Yp(j),Zp(j),aM(j) return EndIf EndIf EndDo Return End c-------------------------------------------------------------------- SUBROUTINE ReadCatalog(Nh) c-------------------------------------------------------------------- INCLUDE 'PMselect.h' Character*120 Line c Lines = 31 Lines = 27 c Lines = 29 Do i=1,Lines read(20,'(A)')Line write(*,'(A)')Line read(21,'(A)')Line c write(*,'(A)')Line write(17,'(A)')Line if(i.le.13.)write(18,'(A)')Line EndDo Nh = 1 nnn = 0 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=105,err=105) Xx,Yy,Zz, & VXp(Nh),VYp(Nh),VZp(Nh), & aM(Nh),Radh(Nh),Vrms(Nh),Vcirc(Nh), & Npart(Nh),Rmax(Nh),Conc(Nh),Nlines(Nh), & error,serN,nlarge(Nh) If(abs(Xx-Xp(Nh))+abs(Yy-Yp(Nh))+abs(Zz-Zp(Nh)).gt.0.1)Then write (*,*) ' mismatch between short and long files' write (*,*) ' Nh=',Nh write (*,'(9G11.3)') Xx,Yy,Zz write (*,'(9G11.3)') Xp(Nh),Yp(Nh),Zp(Nh) stop EndIf 50 format(g11.4,I9,g11.3,g10.3,2f7.1,g10.3,i8,i7,i5,g11.3) 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) c If(Nh.ge.20384) c & write (*,'(" header=",i7,g11.3,i6,3f8.3)') c & Nh,aM(Nh),Nlines(Nh),Xp(Nh),Yp(Nh),Zp(Nh) If(Nlines(Nh).gt.Nrmax)Then write (*,*) ' Error: halo=',Nh,' Nlines=',Nlines(Nh) write (*,'(3f8.3)') Xp(Nh),Yp(Nh),Zp(Nh) write (*,'(i7,2f8.3)') Npart(Nh),Rmax(Nh),Conc(Nh) STOP ' wrong number of profile lines' EndIf If(Nlines(Nh).le.0)Then write (*,*) ' Error: Nh=',Nh,' Nlines=',Nlines(Nh) write (*,*) Xp(Nh),Yp(Nh),Zp(Nh) write (*,*) Npart(Nh),Rmax(Nh),Conc(Nh) STOP ' wrong number of profile lines' EndIf Do i=1,Nlines(Nh) read (20,50) 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(aM(Nh).gt.2.e12.and.aM(Nh).le.4.e12)nnn=nnn+1 if(Nh.ge.Np)goto100 Nh =Nh +1 goto 10 105 write (*,*) ' exit from short file' 100 Nh =Nh -1 write (*,'(" Number of halos read=",T30,i8)') Nh write (*,'(" Number of halos read=",T30,i8)') nnn 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)Then write (*,*) ' Current halo=',ip, ' Max=',Np write (*,*) ' Current box =',i,j,k write (*,*) ' Next halo =',Next STOP 'too many halos' EndIf 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 Nhex = Next write (*,*) ' Nhalos in extended sample=',Nhex RETURN END C---------------------------------- Read in variables REAL FUNCTION INPUT(text) C------------------------------------------------ Character text*(*) write (*,'(A,$)')text read (*,*) x INPUT =x Return End