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

      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
c      Ndisk = Ndisk +Nbulge
 
      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 = 40.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)
         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)
         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)
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)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 integal 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)
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)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

