!--------------------------------------------------
!
!        Read PM data and make 3D overdensity grid
!
!--------------------------------------------------
Module Scale
  Integer*4 :: Nnew,Nold, Ngal
  Real*4    :: dlog = 0.033, threshold = 20., &
               sVrand = 0.333, fraction = 0.02
  Integer*4, parameter :: Ntab = 1000
  Real*8,  dimension(:) :: PDF(-Ntab:Ntab),Delta(-Ntab:Ntab)
  Real*4, ALLOCATABLE, dimension(:,:,:) :: WW
  Real*4, ALLOCATABLE, dimension(:) :: Xgal, Ygal, Zgal, &
                                       VXgal, VYgal, VZgal, Wgal
end Module Scale
!--------------------------------------------------
Program SelectGal
  use Tools
  use Density
  use Scale
  
  Path = ''
  open(17,file='output')
  
  Call Initialize(Path)
   !ALLOCATE(WW(NGRID,NGRID,NGRID))
  print *,' ------------ Init done'
   Call DENSIT
   !Call DENSITscale
  print *,' ------------ DENSIT done'
   !Call PDF_DENSIT(20)
   !Call Transf
   !Call PDF_DENSIT
  !Call SelectGalaxies
  Call SelectParticles
   Call WriteData(2)
 Close(20)
 Close(21)
 
end Program SelectGal

!--------------------------------------------
subroutine Initialize(Path)
!--------------------------------------------
use Tools
use Scale
   logical :: exst
   character(len=*) :: Path
   character*120 :: Fname     

      WRITE (*,'(A,$)') ' Enter snapshot number to analyze => '
      READ (*,*) moment	 ! read snapshot number
      !WRITE (*,'(A,$)') ' Enter new Ngrid                  => '
      !READ (*,*) Nnew	 
      
      CALL ReadDataPM(moment,Path)
      !Nold = NGRID
      !NGRID = Nnew
      myMemory =Memory(1_8*NGRID*NGRID*NGRID)
      Allocate (FI(NGRID,NGRID,NGRID))
      
      write(Fname,'(a,2(i4.4,a))')'DensityStat.',ISTEP,'.',Nrealization,'.dat'
      open(20,file=TRIM(Fname))

      write(20,'(a,f9.5,2(a,i5),a,f8.2,a,i5)') &
           'A=',AEXPN,' Ngrid= ',NGRID, &
           ' Box= ',Box,' Realization= ',Nrealization

end subroutine Initialize

!-------------------------------------------------------------
!
!         Make density contrast FI(NGRID,NGRID,NGRID)
!         for Nparticles XPAR,YPAR,ZPAR(Nparticles) in 1-Ngrid
!
SUBROUTINE DENSITscale
    
!-------------------------------------------------------------
  use Tools
  use Scale
  real*8 :: XN,YN, &
     X,Y,Z,D1,D2,D3,T1,T2,T3,T2W,D2W
integer*8 :: IN

    XN   =Real(NGRID)+1.-1.E-7
    YN   =Real(NGRID)
    Wpar = YN**3/Real(Nparticles,8)
    Xscale = Real(NGRID)/Nold
				!       Subtract mean density
write(*,'(a,i5,a,i12)') ' Init density: Ngrid= ',Ngrid,' Nparticles= ',Nparticles
write(17,'(a,i5,a,i12)') ' Init density: Ngrid= ',Ngrid,' Nparticles= ',Nparticles
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (M1,M2,M3)
    DO M3=1,NGRID
       DO M2=1,NGRID
          DO M1=1,NGRID
	        FI(M1,M2,M3) = 0.
	  END DO
       END DO
    END DO
    print *,'---- Fi=0 done'  
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (IN,X,Y,Z,D1,D2,D3,T1,T2,T3,T2W,D2W) &
!$OMP PRIVATE (I,J,K,I1,J1,K1)    
     DO   IN=1,Nparticles         ! loop over particles 
	  X= (XPAR(IN)-1.)*Xscale +1.
	  Y= (YPAR(IN)-1.)*Xscale +1.
	  Z= (ZPAR(IN)-1.)*Xscale +1.
	  I=INT(X)
	  J=INT(Y)
	  K=INT(Z)

	  D1=X-REAL(I)
	  D2=Y-REAL(J)
	  D3=Z-REAL(K)
	  T1=1.-D1
	  T2=1.-D2
	  T3=1.-D3
	  T2W =T2*WPAR
	  D2W =D2*WPAR
	  I1=I+1
	     IF(I1.GT.NGRID)I1=1
	  J1=J+1
	     IF(J1.GT.NGRID)J1=1
	  K1=K+1
          IF(K1.GT.NGRID)K1=1
!$OMP ATOMIC
          FI(I ,J ,K ) =FI(I ,J ,K ) +T3*T1*T2W
!$OMP ATOMIC
	  FI(I1,J ,K ) =FI(I1,J ,K ) +T3*D1*T2W
!$OMP ATOMIC
	  FI(I ,J1,K ) =FI(I ,J1,K ) +T3*T1*D2W
!$OMP ATOMIC
	  FI(I1,J1,K ) =FI(I1,J1,K ) +T3*D1*D2W
   
!$OMP ATOMIC
	  FI(I ,J ,K1) =FI(I ,J ,K1) +D3*T1*T2W
!$OMP ATOMIC
	  FI(I1,J ,K1) =FI(I1,J ,K1) +D3*D1*T2W
!$OMP ATOMIC
	  FI(I ,J1,K1) =FI(I ,J1,K1) +D3*T1*D2W
!$OMP ATOMIC
	  FI(I1,J1,K1) =FI(I1,J1,K1) +D3*D1*D2W
	   
ENDDO
print *,' Do WW <- FI'
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (M1,M2,M3)
    DO M3=1,NGRID
       DO M2=1,NGRID
          DO M1=1,NGRID
	        WW(M1,M2,M3) = FI(M1,M2,M3)
	  END DO
       END DO
    END DO


  END SUBROUTINE DENSITscale
!-------------------------------------------------------------
!
!         PDF of density field
!         
!
SUBROUTINE PDF_DENSIT(ifile)
    
!-------------------------------------------------------------
  use Tools
  use Scale
  real*8 ::     D1,D2
	
  D1 = 0. ; D2 =0.
  PDF(:)   = 0.
  Delta(:) = 0.
  Av = Nparticles/(float(NGRID))**3
  print *,' Inside PDF '
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (M1,M2,M3,d,ind) REDUCTION(+:D1,D2)
  DO M3=1,NGRID
    DO M2=1,NGRID
    DO M1=1,NGRID
	D1 = D1 + FI(M1,M2,M3)
        D2 = D2 + FI(M1,M2,M3)**2
        d  = FI(M1,M2,M3)
        If(d>1.e-5)Then
         ind = FLOOR(log10(d)/dlog)
         ind = min(max(ind,-Ntab),Ntab)
!$OMP ATOMIC
         Delta(ind) = Delta(ind) + d
!$OMP ATOMIC
         PDF(ind)   = PDF(ind)  + 1
         !if(ind==2.and.Delta(ind)/PDF(ind)<10.**((ind)*dlog)/1.01)Then
         !   print '(3i6)',M1,M2,M3
         !   print '(3es12.4,2i6)',d,10.**((ind)*dlog),Delta(ind)/PDF(ind),ind
         !   stop
         !end if
      end If
     END DO
     END DO
     END DO
    D1 = D1/(Real(NGRID))**3
    D2 = sqrt(D2/(Real(NGRID))**3 -D1**2)
    write(*,'(a,2es13.4)') ' density: average/rms= ',D1,D2
    write(ifile,'(a,2es12.4)') ' density: average/rms= ',D1,D2
    write(ifile,'(a)') ' 1+delta       dN/ddelta/N    dN'
    Do i=-10,Ntab
       If(PDF(i)>100.)Then
          dup  = 10.**((i+1)*dlog)
          down = 10.**((i)*dlog)
          write(ifile,'(3es12.4,i5)') Delta(i)/PDF(i),PDF(i)/(dup-down)/(Real(NGRID))**3,PDF(i)
       end If
    end Do
    
  END SUBROUTINE PDF_DENSIT
!-------------------------------------------------------------
!
!        modify density field by applying a biasing scheme
!         
!
SUBROUTINE Transf
    
!-------------------------------------------------------------
  use Tools
  use Scale
  real*8 ::     D1,D2,D3,d,w,ss,lambda
  Integer(kind=4) :: Npoisson

  Nseed  =1230917
  D1 = 0.
  D2 = 0.

!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (M1,M2,M3,d) REDUCTION(+:D1,D2)
    DO M3=1,NGRID
    DO M2=1,NGRID
    DO M1=1,NGRID
       d  = FI(M1,M2,M3)
       D1 = D1 + d
       D2 = D2 + d**2
    end DO
    end DO
 end DO
 
    D1 = D1/(Real(NGRID))**3
    D2 = sqrt(D2/(Real(NGRID))**3 -D1**2) 
    write(*,'(a,2es13.4)') ' density: average/rms= ',D1,D2
    write(17,'(a,2es12.4)') ' density: average/rms= ',D1,D2 

    ss = 0.
    !If(slope>1.)Then
       D3 = densAv/(pivot*D2)**(slope-1.)
    !else
    !   D3 = densAv
    !end If

       ss  = 0.
       w = 0.001
    DO M3=1,NGRID          !-----------  Bias density to get galaxies
    DO M2=1,NGRID
    DO M1=1,NGRID
       d  = FI(M1,M2,M3)
       !w = d**slope*D3                        !--- galaxy field
       !w = d**slope*0.05                        !--- galaxy field
       !w = d                       !--- galaxy field
       !if(w<0.250)w =0.                       !--- no galaxy if density too small
       !FI(M1,M2,M3) = INT(w)       !--- number of galaxies
       FI(M1,M2,M3) = Npoisson(Nseed,w)       !--- number of galaxies
       !If(FI(M1,M2,M3)<0)print *,' ---- error ---'
       !if(WW(M1,M2,M3)>2.e3)print '(2es12.4,3x,3i5)',WW(M1,M2,M3),FI(M1,M2,M3),M1,M2,M3
       ss = max(ss,FI(M1,M2,M3))
       !if(FI(M1,M2,M3)<0.)print '(2es12.4,3x,3i5)',WW(M1,M2,M3),FI(M1,M2,M3),M1,M2,M3
     END DO
     END DO
  END DO
  
  d  = 0.
  D1 = 0.
  D2 = 0.
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE (M1,M2,M3,w)  REDUCTION(+:D1,D2)
    DO M3=1,NGRID
    DO M2=1,NGRID
       DO M1=1,NGRID
       w  = FI(M1,M2,M3)
       D1 = D1 + w
       D2 = D2 + w**2
    end DO
    end DO
 end DO
    Ngal = D1
    D1 = D1/(Real(NGRID))**3
    D2 = sqrt(D2/(Real(NGRID))**3 -D1**2)
    Ngalaxies = INT(D1)
    print '(a,es12.4)',' maximum density of galaxies = ',ss  
    write(*,'(a,2es13.4)') ' density galaxies: average/rms= ',D1,D2
    write(17,'(a,2es12.4)') ' density galaxies: average/rms= ',D1,D2 

     print *,' Done LOGtrandf'
   END SUBROUTINE Transf

 !-----------------------------------------------
!
!
!
Integer(kind=4) Function Npoisson(Nseed,lambda)
  Real(kind=8) :: U,Prod,e, lambda

  If(lambda<1.e-5)Then
     Npoisson = 0
     return
  end If
  N = -1
  Prod = 1.d0
  e    = exp(-lambda)
  Do while(Prod> e)
     U = RANDd(Nseed)
     Prod = Prod*U
     N = N +1
  end Do
  !print *,N,lambda
  Npoisson = N
end Function Npoisson


FUNCTION RANDd (M) 
!------------------------------------------------                       
      DATA LC, AM, KI, K1, K2, K3, K4, L1, L2, L3, L4	 / 453815927,     &
      2147483648., 2147483647, 536870912, 131072, 256, 	16777216, 4,    &
      16384, 8388608, 128 /
      ML = M / K1 * K1 
      M1 = (M - ML) * L1 
      ML = M / K2 * K2 
      M2 = (M - ML) * L2 
      ML = M / K3 * K3 
      M3 = (M - ML) * L3 
      ML = M / K4 * K4 
      M4 = (M - ML) * L4 
      M5 = KI - M 
      IF (M1.GE.M5) M1 = M1 - KI - 1 
      ML = M + M1 
      M5 = KI - ML 
      IF (M2.GE.M5) M2 = M2 - KI - 1 
      ML = ML + M2 
      M5 = KI - ML 
      IF (M3.GE.M5) M3 = M3 - KI - 1 
      ML = ML + M3 
      M5 = KI - ML 
      IF (M4.GE.M5) M4 = M4 - KI - 1 
      ML = ML + M4 
      M5 = KI - ML 
      IF (LC.GE.M5) ML = ML - KI - 1 
      M = ML + LC 
      RANDd = M / AM 
      RETURN 
    END FUNCTION RANDd
    !--------------------------------------                                 
!				normal random numbers                                              
      FUNCTION GAUSS (M) 
!--------------------------------------                                 
      X = 0. 
      DO I = 1, 5 
      X = X + RANDd (M) 
      EndDo 
      X2 = 1.5491933 * (X - 2.5) 
      GAUSS = X2 * (1. - 0.01 * (3. - X2 * X2) ) 
      RETURN 
    END FUNCTION GAUSS

!-------------------------------------------------------------
!
!                 
!
    SUBROUTINE RandomParticles
      
!-------------------------------------------------------------
  use Tools
  use Scale
  real*8 ::  X,Y,Z,D1,D2,D3
  Integer*8 :: IN, ip

   Ngal = 125000
  
    ALLOCATE(Xgal(Ngal),Ygal(Ngal),Zgal(Ngal))
    ALLOCATE(VXgal(Ngal),VYgal(Ngal),VZgal(Ngal))
    ALLOCATE(Wgal(Ngal))

    Xscale = Real(NGRID)/Nold
    Vscale = 100.*Xscale/AEXPN          ! Scale for velocities
    Ng = 0
    Nseed = 9845619

    print '(a,i11,a,2es12.4)',' Ngalaxies= ',Ngal,' Scales = ',Xscale,Vscale
    DO   IN=1,Nparticles         ! loop over particles
       IF(Randd(Nseed)>fraction)cycle

           Ng = Ng +1
           If(Ng>Ngal)Stop ' Error in number of galaxies in SelectGalaxies'
           Xgal(Ng) = (XPAR(IN)-1.)*Xscale
           Ygal(Ng) = (YPAR(IN)-1.)*Xscale
           Zgal(Ng) = (ZPAR(IN)-1.)*Xscale
           VXgal(Ng) = VX(IN)*Vscale ! + sVrand*GAUSS(Nseed)*WW(I,J,K)
           VYgal(Ng) = VY(IN)*Vscale ! + sVrand*GAUSS(Nseed)*WW(I,J,K)
           VZgal(Ng) = VZ(IN)*Vscale ! + sVrand*GAUSS(Nseed)*WW(I,J,K)
           Wgal(Ng)  = WW(I,J,K) 
           !print '(2i11,3x,3f8.2,f8.1)',IN,Ng,X,Y,Z,FI(I,J,K)
        end DO
      print '(a,i9)',' Number of particles selected = ',Ng
     Ngal = Ng
  end SUBROUTINE 
!-------------------------------------------------------------
!
!                 
!
    SUBROUTINE SelectParticles
      
!-------------------------------------------------------------
  use Tools
  use Scale
  real*8 ::  X,Y,Z,D1,D2,D3
  Integer*8 :: IN, ip

  Ngal = Nparticles*fraction*1.05
  
    ALLOCATE(Xgal(Ngal),Ygal(Ngal),Zgal(Ngal))
    ALLOCATE(VXgal(Ngal),VYgal(Ngal),VZgal(Ngal))
    ALLOCATE(Wgal(Ngal))

    XscaleBox = Box/NGRID
    Vscale = 100.*XscaleBox/AEXPN          ! Scale for velocities
    Ng = 0
    Nseed = 9845619

    print '(a,i11,a,2es12.4)',' Ngalaxies= ',Ngal,' Scales = ',XscaleBox,Vscale
    DO   IN=1,Nparticles         ! loop over particles
 	X= XPAR(IN)  !  grid units
	Y= YPAR(IN)
	Z= ZPAR(IN)
	I=INT(X)
	J=INT(Y)
	K=INT(Z)

        D1=X-REAL(I)            !--- find nearest cell center
        If(D1>0.5)I = I+1
	D2=Y-REAL(J)
        If(D2>0.5)J = J+1
	D3=Z-REAL(K)
        If(D3>0.5)K = K+1
        IF(I.GT.NGRID)I=1
	IF(J.GT.NGRID)J=1
        IF(K.GT.NGRID)K=1

           IF(FI(I,J,K) +1.<threshold)cycle
           IF(Randd(Nseed)>fraction)cycle
           Ng = Ng +1
           If(Ng>Ngal)Stop ' Error in number of galaxies in SelectGalaxies'           
           Xgal(Ng) = (XPAR(IN)-1.)*XscaleBox
           Ygal(Ng) = (YPAR(IN)-1.)*XscaleBox
           Zgal(Ng) = (ZPAR(IN)-1.)*XscaleBox
           VXgal(Ng) = VX(IN)*Vscale ! + sVrand*GAUSS(Nseed)*WW(I,J,K)
           VYgal(Ng) = VY(IN)*Vscale ! + sVrand*GAUSS(Nseed)*WW(I,J,K)
           VZgal(Ng) = VZ(IN)*Vscale ! + sVrand*GAUSS(Nseed)*WW(I,J,K)
           Wgal(Ng)  = FI(I,J,K) +1. 
           !print '(2i11,3x,3f8.2,f8.1)',IN,Ng,X,Y,Z,FI(I,J,K)
        end DO
      print '(a,i9)',' Number of particles selected = ',Ng
     Ngal = Ng
  end SUBROUTINE SelectParticles
!-------------------------------------------------------------
!
!                 
!
SUBROUTINE SelectGalaxies
    
!-------------------------------------------------------------
  use Tools
  use Scale
  real*8 ::  X,Y,Z,D1,D2,D3
  Integer*8 :: IN

  Ngal = 1.1*Ngal
    ALLOCATE(Xgal(Ngal),Ygal(Ngal),Zgal(Ngal))
    ALLOCATE(VXgal(Ngal),VYgal(Ngal),VZgal(Ngal))
    ALLOCATE(Wgal(Ngal))

    Xscale = Real(NGRID)/Nold
    Vscale = 100.*Xscale/AEXPN          ! Scale for velocities
    Ng = 0
    Nseed = 9845619
    print '(a,i11,a,2es12.4)',' Ngalaxies= ',Ngal,' Scales = ',Xscale,Vscale
     DO   IN=1,Nparticles         ! loop over particles 
 	X= (XPAR(IN)-1.)*Xscale +1.  ! rescale to grid units
	Y= (YPAR(IN)-1.)*Xscale +1.
	Z= (ZPAR(IN)-1.)*Xscale +1.
	I=INT(X)
	J=INT(Y)
	K=INT(Z)

        D1=X-REAL(I)            !--- find nearest cell center
        If(D1>0.5)I = I+1
	D2=Y-REAL(J)
        If(D2>0.5)J = J+1
	D3=Z-REAL(K)
        If(D3>0.5)K = K+1
        IF(I.GT.NGRID)I=1
	IF(J.GT.NGRID)J=1
        IF(K.GT.NGRID)K=1

        If(FI(I,J,K)>1.e-3)Then    !--- non-empty cell
           Ng = Ng +1
           If(Ng>Ngal)Stop ' Error in number of galaxies in SelectGalaxies'
           Xgal(Ng) = (XPAR(IN)-1.)*Xscale
           Ygal(Ng) = (YPAR(IN)-1.)*Xscale
           Zgal(Ng) = (ZPAR(IN)-1.)*Xscale
           VXgal(Ng) = VX(IN)*Vscale  + sVrand*GAUSS(Nseed)*WW(I,J,K)
           VYgal(Ng) = VY(IN)*Vscale  + sVrand*GAUSS(Nseed)*WW(I,J,K)
           VZgal(Ng) = VZ(IN)*Vscale  + sVrand*GAUSS(Nseed)*WW(I,J,K)
           Wgal(Ng)  = WW(I,J,K) 
           !print '(2i11,3x,3f8.2,f8.1)',IN,Ng,X,Y,Z,FI(I,J,K)
           FI(I,J,K) = FI(I,J,K) -1
        end If
     end DO
     print '(a,i9)',' Number of galaxies requested= ',Ngal
     print '(a,i9)',' Number of galaxies selected = ',Ng
     Ngal = Ng
  END SUBROUTINE SelectGalaxies
!-------------------------------------------------------------
!
!                 
!
SUBROUTINE WriteData(iSample)
    
!-------------------------------------------------------------
  use Tools
  use Scale
  character :: Fname*120, txt1*5 

  write(Fname,'(3(a,i4.4),a)')'Galaxies.',ISTEP,'.',iSample,'.',Nrealization,'.dat'
  open(50,file=TRIM(Fname),form='unformatted')

  !write(50,'(3a)')      '  X           Y         Z         ', &
  !                      '  Vx         Vy        Vz        Density'
  !Do ip = 1, Ngal
  !   write(50,'(3f10.3,4f10.2)') Xgal(ip), Ygal(ip), Zgal(ip), &
  !                 VXgal(ip), VYgal(ip), VZgal(ip), Wgal(ip)
  !end Do
  write(50) AEXPN,ISTEP,Box,Ngal,NGRID
  write(50) (Xgal(i),i=1,Ngal)
  write(50) (Ygal(i),i=1,Ngal)
  write(50) (Zgal(i),i=1,Ngal)
  write(50) (VXgal(i),i=1,Ngal)
  write(50) (VYgal(i),i=1,Ngal)
  write(50) (VZgal(i),i=1,Ngal)
  write(50) (Wgal(i),i=1,Ngal)
  
  close(50)
  
  write(Fname,'(3(a,i4.4),a)')'GalaxiesShort.',ISTEP,'.',iSample,'.',Nrealization,'.dat'
  open(50,file=TRIM(Fname))
  write(50,'(a,f9.4)') 'Aexpn       = ',AEXPN
  write(50,'(a,i4)')   'ISTEP       = ',ISTEP
  write(50,'(a,f8.2)') 'BOX         = ',BOX
  write(50,'(a,i10)')  'Ngalx       = ',Ngal
  write(50,'(a,i4)')   'Mesh        = ',NGRID
  write(50,'(a,f8.4)') 'fraction    = ',fraction
  write(50,'(a,f8.4)') 'threshold   = ',threshold
  write(50,'(3a)')      '  X           Y         Z         ', &
                        '  Vx         Vy        Vz        Density'
  Do ip = 1, min(100,Ngal)
     write(50,'(3f10.3,4f10.2)') Xgal(ip), Ygal(ip), Zgal(ip), &
                   VXgal(ip), VYgal(ip), VZgal(ip), Wgal(ip)
  end Do
  close(50)
  DEALLOCATE(Xgal,Ygal,Zgal,VXgal,VYgal,VZgal,Wgal)
end SUBROUTINE WriteData