  ! This code randomly distributes mock galaxies for a given model (its redshift distribution and galaxy bias) based on Poisson distribution.
  ! Redshift distribution of galaxies is dndz(z) in Line 73; Galaxy bias is bias(z) in Line 66; Galaxy number density is gal_dens in Line 19.
  ! User may edit these input mode parameters for his/her purpose. 
  ! This code was used in Namikawa & Takahashi (1810.03346, section 3).
  
  implicit none
  integer(4) :: nside, kplane, ngal
  integer(8) :: i, npix, ipix, sum_gal
  real(4) :: z(42)
  real(8) :: mu, mu_mean, gal_dens, tot_ngal, norm, pi
  real(4), allocatable :: delta_shell(:)
  character(len=200) :: fname
 
  pi=4.d0*atan(1.d0)

  nside=4096
  npix=12*nside*nside

  gal_dens=20. ! galaxy number density [1/arcmin^2]
  tot_ngal=4.*pi*(180.*60/pi)**2*gal_dens  ! total galaxy number

  allocate(delta_shell(0:npix-1))

  fname="allskymap_nres12r000.delta_shell.dat"   ! input file name, this file contains surface mass densities of the 42 lens shells from z=0 to 7.1
  open(12,file=fname,status="old",form='unformatted')


  norm=0.d0;
  do i=1,42
    read(12) kplane,z(i)      ! kplane(=1-42) is the index of the lens shell, z(i) is the corresponding redshift, list of the 42 redshifts is given in Table 1 (column of Dumped Redshift; z=0.025, 0.076 and so on) in our paper  
    !write(*,*) kplane,z(i)
    norm=norm+5.004e-2*hubble(z(i))*dndz(z(i))  ! 5.004e-2=(the shell thickness 150Mpc/h)/(the Hubble distance 2997.9Mpc/h)
  enddo
  norm=tot_ngal/norm


  sum_gal=0
  do i=1,42
    read(12) kplane  ! here kplane=i
    read(12) delta_shell(0:npix-1)  ! projected density contrast (=\Sigma/\bar{\Sigma}-1) on the i-th lens shell given in the Healpix ring ordering 
    
    mu_mean=norm*5.004e-2*hubble(z(i))*dndz(z(i))/dble(npix)

    do ipix=0,npix-1
      mu=mu_mean*(1.+bias(z(i))*delta_shell(ipix))

      if(mu.gt.0.) then
        ngal=random_poisson(mu)        
      else
        ngal=0
      end if
      sum_gal=sum_gal+ngal  

      if(ngal.gt.0) write(*,*) kplane,z(i),ipix,ngal   ! ngal is the number of mock galaxies in the pixel of ipix
    end do

  end do
  !write(*,*) "# Ngal=",sum_gal

  deallocate(delta_shell)

  contains


  !! galaxy bias
  real function bias(z)
  real(4) :: z
    bias=sqrt(1.0+z)
  return
  end function

  !! galaxy number distribution as a function of z, dn/dz
  real function dndz(z)
  real(4) :: z,z0,zm=1.
    z0=0.69*zm
    dndz=z*z*exp(-(z/z0)**(1.5))
  return
  end function

  !! Hubble parameter
  real function hubble(z)
  real(4) :: z, omegam=0.279
    hubble=sqrt(omegam*(1.+z)**3+(1.-omegam))
  return
  end function

  !! generate ramdom numbers from Poisson distribution
  integer function random_poisson(mu)
  real(8) :: mu,p,rnd
    p=mu; random_poisson=-1
    do while (p.ge.0.d0)
      call random_number(rnd)
      if(rnd.le.0.d0) call random_number(rnd)
      p=p+log(rnd); random_poisson=random_poisson+1
    end do
  return 
  end function

  end program
