c-------------------------------------------------------------------- c Make a density map using adaptive kernel INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' INCLUDE 'PMmap.h' REAL INPUT,moment Character FileASCII*100 Character File1*100,File2*100,Charyz*2,CharDM*2 EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)), & (Box,extras(100)) C................................................................... C Read data and open files write(*,*)' enter 6 digits of the snapshot and two flags I J ' write (*,*) ' I=projection (0-xy,1-xz), J=Disk(0), DM(1)' read(*,*) moment,iZflag,iDMflag If(moment.lt.1.5)Then write(File1,'("FILES/PMcrda",f6.4,".DAT")')moment write(File2,'("FILES/PMcrs0a",f6.4,".DAT")')moment Else imoment =INT(moment+0.1) c read(*,*) imoment write(File1,'("FILES/PMcrda.",i4.4,".DAT")')imoment write(File2,'("FILES/PMcrs0a.",i4.4,".DAT")')imoment endIf write(*,*) File1 write(*,*) File2 CALL RDTAPE(File1,File2) c aMassDisk =3.00e10 ! Nbulge = 0 c Nbulge = 33409 If(iZflag.eq.0)Then Charyz ='xy' Else Charyz ='xz' EndIf If(iDMflag.eq.0)Then CharDM ='SK' Else CharDM ='DM' EndIf write(FileASCII,'("RESULTS/DensityMap",2A2,".",f6.4,".dat")') & Charyz,CharDM,AEXPN OPEN(37,FILE=FileASCII,STATUS='UNKNOWN') write(FileASCII,'("RESULTS/Result.BAR.",f6.4,".dat")') & AEXPN OPEN(43,FILE=FileASCII,STATUS='UNKNOWN') time = (age(AEXPN)-age(0.6))/1.e9 aMassDisk = aMassDisk/hubble c Nbulge = 33334 write (*,*) ' hubble =',hubble,' mass of disk=',aMassDisk write (*,*) ' RDTAPE is done Time=',time,' a=',AEXPN WRITE (37,100) HEADER,AEXPN,AEXP0,AMPLT, + ASTEP,ISTEP,PARTW,EKIN, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv,Ndisk,Nbulge,aMassDisk,time 100 FORMAT(1X,'Header=>',A45,/ + 1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/ + 1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,/ + 1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I8,/ + 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/ + 1x,' Omega_curvature=',F7.3,/ + 1x,' Ndisk/bulge =',2i8,' Mass disk =',G11.3,' Age=',f7.3,'Gyr') c Box =INPUT(' Enter box size in comoving Mpc/h =') c Box =0.5 ! NG3109 c Box = 2 ! 'bulge' BoxV = Box*100. ! Box size in km/s ScaleV = BoxV/AEXPN/NGRID ! scale factor for Velocities (km/s) ScaleC = Box/NGRID*AEXPN/hubble ! scale factor for Coordinates (real Mpc) ScaleH = 100.*hubble*sqrt(Om0/AEXPN**3+Oml0) ! H(a)=km/s/Mpc if(Ndisk.le.0)write (*,*) ' wrong number of disk particles:',Ndisk if(Ndisk.le.0)Stop ScaleM = aMassDisk/Ndisk Ndisk = Ndisk +Nbulge Nbulge = 0 CALL ReadPnts(N,ScaleC,ScaleV,ScaleM) Call ChangePnts(N,Box,ScaleH) Call GetDensity2(Box) stop end C-------------------------------------------------------------- C get density map in projection SUBROUTINE GetDensity2(Box) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' INCLUDE 'PMmap.h' DIMENSION SURFA(Ndens,Ndens),SURFB(Ndens,Ndens) SAVE SURFA,SURFB EQUIVALENCE (Ndisk,extras(90)) Cell = 2.*RsurfD/Nblink ! cell for linker list Sstep = 2.*RsurfD/Ndens ! cell for surface density xoffset = Box/2.-RsurfD yoffset = Box/2.-RsurfD zoffset = Box/2.-RsurfD dStep = 10.e-6 Do i=1,Nsearch RadDens(i) =dStep*i EndDo write (37,'(" i j SurfDensity Filter(kpc) X Y(kpc)")') write (*,*) ' GetDens: Cell,kpc=',Cell*1.e3,dStep*1.e3 write (*,*) ' Offsets =',xoffset*1.e3,yoffset*1.e3 write (*,*) ' Sstep =',Sstep*1.e3 write (*,'(5f8.3)') (RadDens(i)*1.e3,i=1,Nsearch) Do j=1,Ndens Do i=1,Ndens SURFD(i,j)=0. SURFA(i,j)=0. EndDo EndDo nincr =0 If(iDMflag.eq.0)Then CALL ListDisk(Box) else CALL ListDM(Box) EndIf write (*,*) ' Linker List is done' C Find radius of the filter C Place it to SURFD cC$OMP PARALLEL DO DEFAULT(SHARED) cC$OMP+PRIVATE (i,j,xi,yj,indx,Density,nn,Rsmooth) cC$OMP+SCHEDULE (DYNAMIC,50) Do j=1,Ndens yj = yoffset +(j-1)*Sstep indx =0 if(j.eq.100.or.j.eq.101)indx=1 Do i=1,Ndens xi = xoffset +(i-1)*Sstep CALL Neib2(xi,yj,Density,Rsmooth,indx,nincr,Box) SURFD(i,j) = Rsmooth c write (30,'(6g11.3)')xi,yj,Density,Rsmooth EndDo EndDo write (*,*)' Finding smooting radius Rsmooth is done' write (*,*)' Total number of cells =',Ndens**2 write (*,*)' N cells with insufficently large radius=',nincr Do j=1,Ndens ! Smooth map of Rsmooth Do i=1,Ndens SURFA(i,j) = SURFD(i,j) EndDo EndDo C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (i,j) Do j=2,Ndens-1 Do i=2,Ndens-1 SURFA(i,j) = (SURFD(i,j) + & 0.5*(SURFD(i-1,j)+SURFD(i+1,j)+SURFD(i,j-1)+SURFD(i,j+1))+ & 0.25*(SURFD(i-1,j-1)+SURFD(i+1,j-1) & +SURFD(i-1,j+1)+SURFD(i+1,j+1)))/4. EndDo EndDo write (*,*) ' smoothing of filter radius is done' C Find density using variable filter radius C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (i,j,xi,yj,indx,Density,nn,Rsmooth) C$OMP+SCHEDULE (DYNAMIC,50) Do j=1,Ndens yj = yoffset +(j-1)*Sstep indx =0 c if(j.eq.100.or.j.eq.101)indx=1 Do i=1,Ndens xi = xoffset +(i-1)*Sstep CALL Neib3(xi,yj,Density,SURFA(i,j),indx,Box) SURFD(i,j) = Density EndDo EndDo write (*,*) ' Density is done' Do j=1,Ndens ! Smooth map of Density Do i=1,Ndens SURFB(i,j) = SURFD(i,j) EndDo EndDo Do iter =1,5 C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (i,j) Do j=2,Ndens-1 Do i=2,Ndens-1 SURFB(i,j) = (SURFD(i,j) + & 0.5*(SURFD(i-1,j)+SURFD(i+1,j)+SURFD(i,j-1)+SURFD(i,j+1))+ & 0.25*(SURFD(i-1,j-1)+SURFD(i+1,j-1) & +SURFD(i-1,j+1)+SURFD(i+1,j+1)))/4. EndDo EndDo Do j=1,Ndens Do i=1,Ndens SURFD(i,j) = SURFB(i,j) EndDo EndDo EndDo Do j=1,Ndens yj = (yoffset +(j-1)*Sstep-Box/2.)*1.e3 Do i=1,Ndens xi = (xoffset +(i-1)*Sstep-Box/2.)*1.e3 write(37,'(2i4,2g12.5,4f9.4)')i,j,SURFD(i,j),SURFA(i,j)*1.e3, & xi,yj EndDo EndDo Return End C-------------------------------------------------------------- C Make 2D linker lists of particles in each cell SUBROUTINE ListDisk(Box) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' INCLUDE 'PMmap.h' logical inside EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)) Do i=1,Np Lst(i)=-1 EndDo Do k=Nmlink,Nblink Do j=Nmlink,Nblink Do i=Nmlink,Nblink Label(i,j,k)=0 EndDo EndDo EndDo k=1 Do jp=1,Ndisk ! DISK particles i=INT((Xp(jp)-xoffset)/Cell) i=MIN(MAX(Nmlink,i),Nblink) If(iZFlag.eq.0)Then ! x-y projection j=INT((Yp(jp)-yoffset)/Cell) If(abs(Zp(jp)-zoffset-RsurfD).lt.dslice)Then Ncount = Ncount +1 Lst(jp) =Label(i,j,k) Label(i,j,k) =jp EndIf else ! x-z projection j=INT((Zp(jp)-zoffset)/Cell) If(abs(Yp(jp)-yoffset-RsurfD).lt.dslice)Then Ncount = Ncount +1 Lst(jp) =Label(i,j,k) Label(i,j,k) =jp EndIf EndIf EndDo Return End C-------------------------------------------------------------- C Make 2D linker lists of particles in each cell SUBROUTINE ListDM(Box) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' INCLUDE 'PMmap.h' EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)) Do i=1,Np Lst(i)=-1 EndDo Do k=Nmlink,Nblink Do j=Nmlink,Nblink Do i=Nmlink,Nblink Label(i,j,k)=0 EndDo EndDo EndDo Ncount = 0 xright = xoffset + 2.*RsurfD If(iZFlag.eq.0)Then yright = yoffset + 2.*RsurfD Else yright = zoffset + 2.*RsurfD EndIf zleft = Box/2.-dslice zright = Box/2.+dslice k=1 Do jp=1+Ndisk+Nbulge,lspecies(Nspecies) ! DM particles inside = (Xp(jp).gt.xoffset.and.Xp(jp).lt.xright) IF(iZFlag.eq.0)Then ! x-y projection inside=inside.and.(Yp(jp).gt.yoffset.and.Yp(jp).lt.yright) inside=inside.and.(Zp(jp).gt.zleft.and.Zp(jp).lt.zright) Else inside=inside.and.(Zp(jp).gt.zoffset.and.Zp(jp).lt.yright) inside=inside.and.(Yp(jp).gt.zleft.and.Yp(jp).lt.zright) EndIf If(inside)Then i=INT((Xp(jp)-xoffset)/Cell) i=MIN(MAX(Nmlink,i),Nblink) If(iZFlag.eq.0)Then ! x-y projection j=INT((Yp(jp)-yoffset)/Cell) Ncount = Ncount +1 Lst(jp) =Label(i,j,k) Label(i,j,k) =jp else ! x-z projection j=INT((Zp(jp)-zoffset)/Cell) Ncount = Ncount +1 Lst(jp) =Label(i,j,k) Label(i,j,k) =jp EndIf EndIf c write (*,'(" xyz=",3g11.3," dy=",3g11.3)') c & (Xp(jp)-Box/2.)*1.e3,(Yp(jp)-Box/2.)*1.e3,(Zp(jp)-Box/2.)*1.e3, c & (Yp(jp)-yoffset-RsurfD)*1.e3, c & (Xp(jp)-xoffset-RsurfD)*1.e3, c & (Zp(jp)-zoffset-RsurfD)*1.e3 EndDo write (*,*) ' Number of particles in List=',Ncount Return End C-------------------------------------------------------------- C Find smoothing radius Rsmooth for center C C input: xc,yc - center C output: Rsmooth (Mpc units), Density estimate of density C nincr = number of cells where radius was not C sufficiently large SUBROUTINE Neib2(xi,yj,Density,Rsmooth,indx,nincr,Box) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' INCLUDE 'PMmap.h' DIMENSION iDen(Nsearch) C Initiate counters iFlag = 0 Radius = RadDens(15) Radius2= Radius**2 50 Do i=1,Nsearch iDen(i) =0 EndDo c write (*,*) '==== Neib:',Xc,Yc,Zc,dStep c limits for Label i1=INT((xi-xoffset-Radius)/Cell) j1=INT((yj-yoffset-Radius)/Cell) k1=1 i1=MIN(MAX(Nmlink,i1),Nblink) j1=MIN(MAX(Nmlink,j1),Nblink) k1=MIN(MAX(Nmlink,k1),Nblink) i2=INT((xi-xoffset+Radius)/Cell) j2=INT((yj-yoffset+Radius)/Cell) k2=1 i2=MIN(MAX(Nmlink,i2),Nblink) j2=MIN(MAX(Nmlink,j2),Nblink) k2=MIN(MAX(Nmlink,k2),Nblink) c write (*,'(" limits=",6i5)')i1,i2,j1,j2,k1,k2 C Look for neibhours Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) 10 If(jp.ne.0)Then If(iZFlag.eq.0)Then ! x-y projection dd =(xi-Xp(jp))**2+(yj-Yp(jp))**2 Else dd =(xi-Xp(jp))**2+(yj-Zp(jp))**2 EndIf If(dd.lt.Radius2.and.abs(Yp(jp)-Box/2.).lt.dZ_slice)Then c If(dd.lt.Radius2)Then disd = sqrt(dd)/dStep ind = min(INT(disd) +1,Nsearch) iDen(ind) = iDen(ind) +1 EndIf jp =Lst(jp) GoTo10 EndIf EndDo EndDo EndDo Do i=2,Nsearch ! make iDen the integral number iDen(i) =iDen(i) +iDen(i-1) EndDo ! Increase radius if there were not enough particles If(iDen(Nsearch).lt.nfilter.and.iFlag.eq.0)Then nincr = nincr +1 iFlag = 1 c write (*,'(" Increase search rad: N=",i5," R=",6g11.3)') c & iDen(Nsearch),Radius*1.e3,xi*1.e3,yj*1.e3,Cell*1.e3 Radius = RadDens(Nsearch) Radius2= Radius**2 GoTo 50 EndIf c If(iDen(Nsearch).lt.nfilter) c & write (*,'(" too small Rad=",3g11.3,2i4)') c & Radius*1.e3,xi*1.e3,yj*1.e3,iDen(Nsearch),iFlag Do i=1,Nsearch nn =iDen(i) If(nn.ge.nfilter)Then if(i.eq.1)Then Rsmooth = RadDens(1)*sqrt(nfilter/float(nn)) else Rsmooth =RadDens(i-1)+(RadDens(i)-RadDens(i-1))* & sqrt((nfilter-iDen(i-1))/ & float(max(iDen(i)-iDen(i-1),1))) endif c Rsmooth = min(max(Rsmooth,Sstep),RadDens(Nsearch)) Density = nfilter/Rsmooth**2 goto 20 EndIf EndDo Rsmooth =RadDens(Nsearch) Density = nn/RadDens(Nsearch)**2 20 Density =Density/4.188e9 ! number of particles per kpc**3 c If(indx.eq.1)Then c write(*,'(/5x,"Dens=",2g11.3,6i8)')Density,Rsmooth*1.e3,i,nn c write (*,'(25i7)')(i,i=1,Nsearch) c write (*,'(25i7)')iDen c do i=2,Nsearch c iDen(i) = iDen(i) +iDen(i-1) c EndDo c write (*,'(25i7)')iDen c endIf Return End C-------------------------------------------------------------- C Find density for a center with given smoothing C radius Rsmooth. Use Gaussian kernel with width C sigma =alpha*Rsmooth, alpha =1/3-1/2 C input: xc,yc - center C output: Density estimate of density SUBROUTINE Neib3(xi,yj,Density,Rsmooth,indx,Box) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' INCLUDE 'PMmap.h' REAL*8 W C Initiate counters c Radius = Rsmooth*2.5 Radius = Rsmooth*1.5 Radius2= Radius**2 sigma = Rsmooth aNorm = 1./(6.283185*sigma**2*(1.-exp(-(Radius/sigma)**2/2.))) W = 0. nn = 0 s2 = 2.*sigma**2 c write (*,*) '==== Neib:',Xc,Yc,Zc,dStep c limits for Label i1=INT((xi-xoffset-Radius)/Cell) j1=INT((yj-yoffset-Radius)/Cell) k1=1 i1=MIN(MAX(Nmlink,i1),Nblink) j1=MIN(MAX(Nmlink,j1),Nblink) k1=MIN(MAX(Nmlink,k1),Nblink) i2=INT((xi-xoffset+Radius)/Cell) j2=INT((yj-yoffset+Radius)/Cell) k2=1 i2=MIN(MAX(Nmlink,i2),Nblink) j2=MIN(MAX(Nmlink,j2),Nblink) k2=MIN(MAX(Nmlink,k2),Nblink) c write (*,'(" limits=",6i5)')i1,i2,j1,j2,k1,k2 C Look for neibhours Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) 10 If(jp.ne.0)Then If(iZFlag.eq.0)Then ! x-y projection dd =(xi-Xp(jp))**2+(yj-Yp(jp))**2 Else dd =(xi-Xp(jp))**2+(yj-Zp(jp))**2 EndIf If(dd.lt.Radius2.and.abs(Yp(jp)-Box/2.).lt.dZ_slice)Then c If(dd.lt.Radius2)Then W =W +exp(-dd/s2) nn = nn +1 EndIf jp =Lst(jp) GoTo10 EndIf EndDo EndDo EndDo 20 Density =aNorm*W/1.e9 If(indx.eq.1)Then write(*,'(5x,"Dens,Rs,W,xy,nn=",5g11.3,6i8)') & Density,Rsmooth*1.e3,W,xi*1.e3,yj*1.e3,nn endIf Return End