!-------------------------------------------------------------------- ! Read Catshort.DAT ! Construct halos averaged properties !-------------------------------------------------------------------- Module SetUp PARAMETER (Nmax =3000000 ) PARAMETER (Nshells = 150 ) PARAMETER (hubble = 100.d0 ) PARAMETER (Box = 20.d0 ) PARAMETER ( aMasslimit = 1.d10 ) PARAMETER ( RadLmax = 1.0d0 ) ! isolation creterion PARAMETER ( Vclimit = 0.d0 ) ! Read and use halos above this limit PARAMETER ( Vratio = 1.5 ) ! ratio for isolated halos PARAMETER ( Fisolate = 3.0 ) ! limit Vc for isolated halos PARAMETER ( dR = 10. ) ! buffer width for periodic boundaries PARAMETER ( pi = 3.14159265) Real*4, dimension(Nmax) :: xc,yc,zc, & vxc,vyc,vzc, & vc=0., vrms=0., aMass=0., & Rvir=0.,Nbins=0,Conct=0. Integer, dimension(Nmax) :: Liso =1 Real*4, dimension(0:Nshells) :: rdsh,aMsh,Vdsh,Vcsh Integer, dimension(0:Nshells) :: nsh=0 Integer :: Nh,Nhext,Nclusters Character Line*79,LineZ*45 End Module SetUp !-------------------------------------------------------------------- ! Find isolated halos in a given mass range in Catalog.DAT ! Construct their averaged density and vel profile !-------------------------------------------------------------------- PROGRAM HaloProfiles use SetUP CHARACTER Name*80,sName*2 ! Open(1,file='Catshort.00195.DAT') ! Open(1,file='CatshortB.490.DAT') Open(1,file='Catshort.00986.DAT') Vc1 =600. Vc2 =900. write (*,'(A,$)')' Enter range of circular velocities=>' read (*,*) Vc1,Vc2 write (*,'(A,$)')' Enter 2 characters for file name=>' read (*,*) sName write(Name,'(3a,i4.4,a)')'results.',sName,'.',INT(Vc2),'.dat' Open(2,file=Name) ! Write headers write (2,'(T10,"Hubble=",f8.2,T25,"Box=",f8.2,"Mpc/h", & T35,"Vc limits=",3f8.2,"km/s")') & hubble,Box,Vc1,Vc2,Vclimit write (2,'(T10,"Isolation radius=",f8.2)') & RadLmax Lines = 27 ! skip header of catalog file Do i=1,Lines read(1,'(A)')Line if(i.le.40)write(2,*)Line if(i.le.40)write(*,*)Line EndDo Nread =0 Nh = 1 ! read halo catalog 10 read(1,*,end=100,err=100)xc(Nh),yc(Nh),zc(Nh), & vxc(Nh),vyc(Nh),vzc(Nh),aMass(Nh), & Rvir(Nh),Vrms(Nh),vc(Nh),nn,Rmax,conct(Nh) !,Nbins(Nh) Nread =Nread +1 if(mod(Nh,10000)==0) write(*,'(2i9,3f9.3,i4)') Nread,Nh,vc(Nh),conct(Nh) !,Nbins(Nh) if(vc(Nh).lt.Vclimit)goto 10 ! do not take halo if too small if(Nh.ge.Nmax)goto 100 !write(*,*) Nh,Nbins(Nh),vc(Nh) Nh =Nh +1 goto 10 100 close(1) Nh =Nh -1 write (2,'(" Number of halos read=",T30,i8)')Nh write (*,*) '<====== Number of halos read=',Nh write (*,30) (i,xc(i),yc(i),zc(i),aMass(i),vc(i),i=1,1) write (*,30) (i,xc(i),yc(i),zc(i),aMass(i),vc(i),i=Nh,Nh) 30 format(' halo=',i8,' x=',3g11.3,' M=',g11.3,' Vc=',g11.3) CALL ExpandSample ! CALL Isolated3D(Vc1,Vc2) CALL GetProfiles stop end PROGRAM HaloProfiles !-------------------------------------------------------------------- ! make sample larger by replicating ! some halos using periodical ! boundary conditions ! dR is the width of the buffer in Mpc/h ! Nhext = new number of halos SUBROUTINE ExpandSample !-------------------------------------------------------------------- use SetUp Logical inx,iny,inz Xleft =-dR ! new boundaries for the catalog Xright =Box+dR Yleft =-dR Yright =Box+dR Zleft =-dR Zright =Box+dR write (*,*) Xleft ,Xright Nhext = Nh Do ip =1,Nh ! go over all halos x =xc(ip) y =yc(ip) z =zc(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 Nhext =Nhext +1 if(Nhext.gt.Nmax) STOP 'too many halos' xc(Nhext) = x+dx yc(Nhext) = y+dy zc(Nhext) = z+dz vxc(Nhext) = vxc(ip) vyc(Nhext) = vyc(ip) vzc(Nhext) = vzc(ip) aMass(Nhext) = aMass(ip) vc(Nhext) = vc(ip) Rvir(Nhext) = Rvir(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=',Nhext write (2,'(" Nhalos in extended sample=",T30,i8)')Nhext RETURN END SUBROUTINE ExpandSample !-------------------------------------------------------------------- ! select isolated halos ! range of circular velocities Vc1=",g12.4," rvir=",f7.1," =",f7.1)') & Nx,Vc1,Vc2, Fisolate,aMav/Nx, rvirav/Nx, vcav/Nx RETURN End SUBROUTINE Isolated3D !-------------------------------------------------------------------- ! Cosnstruct averaged profiles of isolated halos ! SUBROUTINE GetProfiles !-------------------------------------------------------------------- use SetUp real*4, DIMENSION(Nshells) :: Ss,Slope,dv2,beta rdsh =0. ; aMsh=0. ; Vdsh =0. ; Vcsh =0. ; nsh =0 dv = 0.1 Nx =0 ! count halos. Find number of bins Nb aMav =0. rvirav =0. vcav =0. vMax =0. x =12.540; y=9.401; z =9.181 Do ip =1,Nh dd = sqrt((xc(ip)-x)**2+(yc(ip)-y)**2+(zc(ip)-z)**2) If(Liso(ip).eq.1.and.dd> 0.20)Then Liso(ip)=0. EndIf Enddo Do ip =1,Nh If(Liso(ip).eq.1)Then Nx =Nx +1 aMav =aMav +aMass(ip) rvirav=rvirav +Rvir(ip) vcav =vcav +vc(ip) vMax =max(vMax,vc(ip)) EndIf EndDo aMav =aMav /Nx rvirav =rvirav/Nx vcav =vcav /Nx write(*,*) ' Max Vc =',vMax,Nx Do ip =1,Nh ! get profiles If(Liso(ip).eq.1)Then ii =min(INT(log10(vc(ip))/dv),Nshells) Vcsh(ii) = Vcsh(ii) +vc(ip) aMsh(ii) = aMsh(ii) +aMass(ip) nsh(ii) = nsh(ii) +1 EndIf EndDo ! print results write(2,'(3x,a,T25,a,T40,a,T55,a,T70,a,T85,a,T100,2a)') & 'Vcirc','Nhalos','f(Vcirc)','Mass' iFlag =0 Do j=1,Nshells nn = max(nsh(j),1) nup = 0 do k=j,Nshells nup =nup+nsh(k) enddo aMsh(j) = aMsh(j)/nn Vcsh(j) = Vcsh(j)/nn if(nsh(j).ne.0)iFlag=1 If(iFlag==1)Then If(nn<3)Vcsh(j) =10.**((j+0.5)*dv) vv =10.**((j+1)*dv) - 10.**((j)*dv) if(Vcsh(j)