c-------------------------------------------------------------------- c Satellites in N-body simulations c Make dV-dR diagram c-------------------------------------------------------------------- include 'subhalos.h' INTEGER listfiles(4) Character Name*80 DATA dR/2.d0/ ! Mpc/h: width of buffer for volume replication DATA listfiles/1,30,31,32/ c Open(1,file='hlist_0.9996.20.dat') c Open(1,file='CatshortA.DAT') write (*,*) ' Enter 3 digits for snapshot ' read (*,*) moment write(Name,'(a,i3.3,a)')'Subhalos.',moment,'.dat' Open(1,file=Name) write(Name,'(a,i3.3,a)')'Catshort.',moment,'.DAT' Open(10,file=Name) CALL ReadParticles(dR) 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) Do i=1,1 write (listfiles(i),'(T10,"Hubble=",f8.2,T28,"Box=",f8.2,"Mpc/h", & /T10,"Selection: Mass limit=",g11.3,"Msun/h", & T55,"Vc limit=",f8.2,"km/s", & T77,"width of buffer=",f8.2,"Mpc/h")') & hubble,Box, aMasslimit,Vclimit, & dR EndDo CALL ListPrimaries CALL Subhalos(1) stop end c-------------------------------------------------------------------- c read particlesNx,Ny,Nz and c replicate in buffer of width dR Mpc SUBROUTINE ReadParticles(dR) c-------------------------------------------------------------------- include 'subhalos.h' Logical inx,iny,inz Lines =32 Do i=1,Lines read(10,'(A)')Line If(i.le.3)write(1,*)Line EndDo Nh = 1 10 read(10,*,end=100,err=100)xc(Nh),yc(Nh),zc(Nh), & vxc(Nh),vyc(Nh),vzc(Nh), & aMass(Nh),rhalo(Nh),vrms(Nh),vc(Nh) if(Nh.ge.Nmax)goto100 Nh =Nh +1 goto 10 100 Nh =Nh -1 NpBox = Nh write (1,'(" Number of halos read=",T30,i8)') Nh write (*,*) '<====== Number of particles read=',Nh write (*,30) (i,xc(i),yc(i),zc(i),i=1,2) write (*,30) (i,xc(i),yc(i),zc(i),i=Nh-3,Nh) 30 format(' halo=',i8,' x=',3g11.3) c now relpicate periodically Xleft =-dR ! new boundaries for the catalog Xright =Box+dR Yleft =-dR Yright =Box+dR Zleft =-dR Zright =Box+dR write (*,*) Xleft ,Xright Ne = 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 Ne =Ne +1 if(Ne.gt.Nmax) STOP 'too many particles' xc(Ne) =x+dx yc(Ne) =y+dy zc(Ne) =z+dz vxc(Ne) =vxc(ip) vyc(Ne) =vyc(ip) vzc(Ne) =vzc(ip) vc(Ne) =vc(ip) vrms(Ne)=vrms(ip) aMass(Ne)=aMass(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=',Ne write (1,'(" Nhals in extended sample=",T30,i8)')Ne Next =Ne ! Np is the number of particles in extened set close (10) RETURN END c-------------------------------------------------------------------- c find 'clusters': halos with large Vcirc c SUBROUTINE ListPrimaries c-------------------------------------------------------------------- include 'subhalos.h' Nprimaries = 0 Do ip =1,Next ! look for primaries in extended halo list If(vc(ip).gt.Vcl)Then Nprimaries = Nprimaries + 1 ListCl(Nprimaries) = ip EndIf EndDo write (*,'(" N primaries=",i6," with Vc >=",2f8.3)')Nprimaries,Vcl RETURN End c-------------------------------------------------------------------- c SUBROUTINE Subhalos(iHalo) c-------------------------------------------------------------------- include 'subhalos.h' Logical inx,iny,inz DIMENSION Nprof(0:Nshells),rbin(0:Nshells), & Nvd(0:Nshells) dlogR = 0.1 ! in units of Mpc/h dVc = 2.5 ! in km/s If(iHalo.gt.Nprimaries)Then write(*,'(a,2i5)')' Too large index for halo:',iHalo,Nprimaries return EndIf Nx = 0 N2 = 0 N30= 0 N40= 0 xx = xc(ListCl(iHalo)) yy = yc(ListCl(iHalo)) zz = zc(ListCl(iHalo)) rvir = rhalo(ListCl(iHalo)) vcprime =vc(ListCl(iHalo)) write (1,'(T7,a,T30,a,T55,a,T72,a,T89,a,T101,a,T113,a)') & 'Number','Coord(Mpc/h)','Vcirc(kms)','Vrms3D','M(sun/h)', & 'Rvir(kpc/h)' write (1,'(a,i5,7g14.4)')' Halo=',ListCl(iHalo),xx,yy,zz, & vcprime,vrms(ListCl(iHalo)), & aMass(ListCl(iHalo)),rvir Do i=0,Nshells Nprof(i) =0 Nvd(i) =0 rbin(i) =0. EndDo Do ip =1,Next x =xc(ip) y =yc(ip) z =zc(ip) vv =vc(ip) aM =aMass(ip) If(vv.gt.Vclimit.and.am.gt.aMasslimit)Then dd = sqrt((x-xx)**2+(y-yy)**2+(z-zz)**2)*1000. indx = MIN(INT(LOG10(dd)/dlogR),Nshells) If(dd.le.rvir)Nx = Nx +1 Nprof(indx) = Nprof(indx) +1 rbin(indx) = rbin(indx) +dd c If(dd.le.rvir.and.vv.gt.40.) c & write (*,'(3f8.4,3g12.4)')x,y,z,vv,aM,dd If(dd.lt.30.)Then N30 =N30+1 write (*,'(3f8.4,3g12.4)')x,y,z,vv,aM,dd EndIf If(dd.lt.40.)N40 =N40+1 If(dd.le.2.*rvir)Then N2 =N2+1 iv =MIN(INT(vv/dVc),Nshells) Nvd(iv) = Nvd(iv)+1 EndIf EndIf EndDo ! ip loop v0 = 0. write (1,'(a,i4,a,i4)')' Halos inside Rvir=',Nx, & ' Halos inside 2Rvir=',N2 write (1,'(" R(kpc/h) R/Rvir dN/dV N")') write (*,*) ' Inside 30kpc n=',N30 write (*,*) ' Inside 40kpc n=',N40 Do i=1,Nshells rout = 10.**(dlogR*i) rin = 10.**(dlogR*(i-1)) v1 = 4.1888*rout**3 v0 = 4.1888*rin**3 dv = v1-v0 Nbin = Nprof(i-1) dd = Nbin/dv rad= rbin(i-1)/max(Nbin,1) if(Nbin.le.3)rad=(rout+rin)/2. if(rout.lt.1.e3.and.Nbin.gt.0.and.rin.gt.20.) & write (1,'(2f8.3,g12.3,i5,3g12.3)') & rad,rad/rvir, & dd,Nprof(i-1),rin,rout EndDo write (1,'(" V_c(km/s) N(>Vc) N(Vc) dV=",f8.3)')dVc Nx =0 Do i=Nshells-1,0,-1 Vout = dVc*(i+1) Vin = dVc*i If(Vin.lt. 150.)Then Nx = Nx + Nvd(i) write (1,'(f8.3,2i5)') & Vin,Nx,Nvd(i) EndIf EndDo RETURN End