ccc   gadget2psis_v1.1.f   Time-stamp: <2016-06-07 16:30:08 hamana>
!
! This is based on gadget2psis_v1.0.f with some modifications made.
! This assumes that input Nbody data is in the form of snapshot files of Gadget2.
! 
!* required library:
!  Healpix (gadget2psis_v1.1.f was developed and tested with Healpix_3.11)
!  cfitsio
!
!* compile with gfortran: 
! gfortran gadget2psis_v1.1.f -o gadget2psis.exe -O3 -I/usr/local/Healpix_3.11/include -DGFORTRAN -fno-second-underscore -fopenmp -fPIC -L/usr/local/Healpix_3.11/lib -L/usr/local/lib -lhealpix -lhpxgif -lcfitsio -Wl,-R/usr/local/lib
!
!!! usage:
! gadget2psis.exe [options]
! [options] are:
! -path : PATH to input Gadget2 snapshot files
! -snap : Gadget snapshot number (integer >= 0)
! -np   : [number of particles of Nbody data]^{1/3}
! -nres : Healpix's parameter "nres"
! -d    : path to output dir
! -f    : output file name
! -th   : thickness of shell [Mpc/h]
! -kp   : lens-plane number (integer >= 1 ***closest is 1***) 
! -wde  : w_DE
! -nt   : number of openmp threads
!
!!! example:
! ./gadget2psis.exe -path DATA/450Mpc/snapdir_002 -snap 2 -np 2048 -nres 12 -d DATA -f b450r004p001 -th 150.0 -kp 1 -wde -1.0 -nt 8
! ./gadget2psis.exe -path DATA/450Mpc/snapdir_001 -snap 1 -np 2048 -nres 12 -d DATA -f b450r004p002 -th 150.0 -kp 2 -wde -1.0 -nt 8
! ./gadget2psis.exe -path DATA/450Mpc/snapdir_000 -snap 0 -np 2048 -nres 12 -d DATA -f b450r004p003 -th 150.0 -kp 3 -wde -1.0 -nt 8
!
!!! copyright --- Takashi Hamana 
! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      MODULE cosmoparam
      real(8) omegam
      real(8) wde
      END MODULE cosmoparam
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      program gadget2psis
      use alm_tools
      use healpix_types
      use pix_tools
      use udgrade_nr
      use fitstools, only : output_map
      use head_fits 
      use cosmoparam

      implicit real (a-h,o-z)
      implicit integer (i-n)     
      
!!! valiables for read_snapshot
      integer(4) :: SNAPSHOT ! number of dump

      character(len=256) :: path
      character(len=256) :: filename, dirname, snumber, fnumber
      
      integer(4) :: fn,i,nstart
      integer(4) :: j,k,nlow

      integer(4) :: nsnap
      real(4)    :: pos(1:3)

!!! header structure in c
!struct io_header_1
!{
!  int      npart[6];
!  double   mass[6];
!  double   time;
!  double   redshift;
!  int      flag_sfr;
!  int      flag_feedback;
!  int      npartTotal[6];
!  int      flag_cooling;
!  int      num_files;
!  double   BoxSize;
!  double   Omega0;
!  double   OmegaLambda;
!  double   HubbleParam; 
!  char     fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8];  /* fills to 256 Bytes */
!} header1;
!*** transforming to Fortran
      integer(4) :: npart(0:5)
      real(8)    :: massarr(0:5)
      real(8)    :: a
      real(8)    :: redshift
      integer(4) :: flag_sfr
      integer(4) :: flag_feedback
      integer(4) :: nall(0:5)
      integer(4) :: flag_cooling
      integer(4) :: FILES
      real(8)    :: boxsize ! shoube Mpc/h unit ***check this*** 
!      real(8)    :: omegam ! defined in MODULE cosmoparam
      real(8)    :: OmegaLambda
      real(8)    :: HubbleParam 
      integer(4) :: unused(24)
      
!!! variables for Healpix
      integer(i4b) :: nres,nside,lmax,hporder,hpspin 
      integer(i8b) :: npix,ipring,ipnest,ipix,np
      real(dp) :: hpvector(3)
      real(sp),allocatable :: hpmap_phi(:),hpmap_alp(:,:)
      real(sp),allocatable :: hpmap_mag(:,:)
      real(sp),allocatable :: hpmap(:,:)
      
      character(len=80), dimension(1:64) :: mapheader
      real(dp), allocatable, dimension(:,:) :: dw8 
      real(dp), dimension(2) :: hpz 
      complex(spc), allocatable, dimension(:,:,:) :: hpalm
      real(sp),allocatable :: hpcl(:,:)

!!! variables for Ray-tracing
      real(sp),allocatable :: sdens(:)
      
      real(4) :: rmin2,rmax2,pixpart
      real(8) :: sdensnorm
      real(4) :: alens

      character(len=256) :: pnumber
      character(len=256) :: pfname
      character(len=256) :: fname,indir,outfile,outdir,fin,filein,outlog

      character*256 opt,arg

      include 'omp_lib.h'
      
!
!!!!!! definitions  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
      chubble0=2.9979e3  ! c/H0 in Mpc/h (H0=100h km/sex/Mpc)

!!! Healpix related
      hporder=1 ! integer =1 for ring, =2 for nested
      hpz(1)=-1.0_dp
      hpz(2)=1.0_dp

!
!!!!!!! set up parameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
      narg=iargc()
c      write(*,*)narg
!!! cosmological parameters *** FLAT COSMOLOGY ONLY ***
      wde=-1.0

!!! Gadget related
      path='/data1/DATA_HSC/SNAP/L450/004/snapdir_000'
      SNAPSHOT=0
      np=1024

!!! Healpix related
      nres=12

!!! output file names
      outdir='/home/hamana/T2K_LENS/PSIS'
      fname='test'

!!! RayTrace related 
      thick=150.0 !Mpc/h
      kplane=2

!!! number of threds
      nopenmp=8

      do i=1,narg
         call getarg(i,opt)
         call getarg(i+1,arg)
         select case (opt)
         case ('-path')
            path=arg
         case ('-snap')
            read(arg,*)SNAPSHOT  ! gadget snapshot number
         case ('-np')
            read(arg,*)np ! [number of particles]^{1/3} = ll^3
         case ('-nres')
            read(arg,*)nres  ! Healpix parameter
         case ('-d')
            outdir=arg ! output dir
         case ('-f')
            fname=arg ! output file name
         case ('-th')
            read(arg,*)thick  ! RayTrace parameter
         case ('-kp')
            read(arg,*)kplane  ! RayTrace parameter *** >0 ***
         case ('-wde')
            read(arg,*)wde
         case ('-nt')
            read(arg,*)nopenmp
         end select
      enddo
!!! set parameters for openmp
!$    call omp_set_num_threads(nopenmp)

      if (kplane.le.0) then
         write(*,*)'*** ERROR ***'
         write(*,*)'*** kp MUST be > 0 ***'
         stop
      endif

!!! Healpix related
      nside=2**nres
      npix=12_i8b*nside*nside
      write(*,*)'#'
      write(*,*)'### HEALPIX parameters'
      write(*,*)'#'
      write(*,*)'nres  =',nres
      write(*,*)'nside =',nside
      write(*,*)'npix  =',npix
      
      allocate(sdens(0:npix-1))
!$omp parallel private(l)
!$omp do
      do ipix=0,npix-1
         sdens(ipix)=0.0
      enddo
!$omp end do
!$omp end parallel

!!! Gadget related file path
      path=path(1:len_trim(path))
 
!!! RayTrace related
      rmin=(kplane-1)*thick
      rmax=kplane*thick
      rmin2=rmin*rmin
      rmax2=rmax*rmax
            
!
!!!!!! read Gadget snapshot !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! first we just see how many particles there are of each type
      write(snumber,'(I3)')SNAPSHOT
      write(fnumber,'(I3)')0
      do i=1,3
         if (snumber(i:i) .eq. ' ') then
            snumber(i:i)='0'
         end if
      end do

      dirname=path(1:len_trim(path))
      if (FILES.eq.1) then
         filename= dirname(1:len_trim(dirname)) // '/snapshot_' // 
     $        snumber(1:len_trim(snumber))
      else
         filename= dirname(1:len_trim(dirname)) // '/snapshot_' // 
     $        snumber(1:len_trim(snumber)) // '.' // 
     $        fnumber(verify(fnumber,' '):3)
      endif

      write(*,*)'#'
      write(*,*)'### reading GADGET header from snapshot :',
     $     filename(1:len_trim(filename))
      write(*,*)'#'
! now, read in the header
      open (1, file=filename, form='unformatted')
      read(1)npart,massarr,a,redshift,flag_sfr,flag_feedback,nall
     $     ,flag_cooling,FILES,boxsize,omegam,OmegaLambda,HubbleParam
     $     ,unused
      close(1)
      write(*,*)'npart =',npart
      write(*,*)'massarr =',massarr
      write(*,*)'a =',a
      write(*,*)'redshift =',redshift
      write(*,*)'flag_sfr =',flag_sfr
      write(*,*)'flag_feedback =',flag_feedback
      write(*,*)'nall =',nall      
      write(*,*)'flag_cooling  =',flag_cooling
      write(*,*)'FILES  =',FILES
      write(*,*)'boxsize  =',boxsize
      write(*,*)'omegam  =',omegam
      write(*,*)'OmegaLambda  =',OmegaLambda
      write(*,*)'HubbleParam  =',HubbleParam
      write(*,*)'unused =',unused
      
!!! setup Gadget related parameters
      boxsize=boxsize/1000.0 !kpc/h => Mpc/h
      boxinv=1.0/boxsize
      boxvol=boxsize**3 ! in (Mpc/h)^3

!!! setup raytracing related parameters
      shellvol=4.0*PI*(rmax**3-rmin**3)/(8.0*3.0) !!! 1/8 because of shell-volume in Octant !!! 
      shellpart=np**3*shellvol/boxvol
      pixpart=shellpart/(npix/8.0)
      dcent=(kplane-0.5)*thick ! Mpc/h
      dlplane=0.75*(rmax**4-rmin**4)/(rmax**3-rmin**3) ! particle number weighted d_lensplane in Mpc/h ! change made in v1.1     

      write(*,*)'#'
      write(*,*)'### RayTracing parameters'
      write(*,*)'#'
      write(*,*)'kplane =',kplane
      write(*,*)'shell center =',dcent
      write(*,*)'weighted dlens =',dlplane
      write(*,*)'pixpart =',pixpart

      ain=a
      dlens=dlplane/chubble0
      call newlap(dlens,aout,ain)
      alens=aout ! lens redshift
      zlens=1.0/alens-1.0
      write(*,*)'z_plane =',zlens
      sdensnorm=1.5*omegam*boxvol*npix/
     $     (alens*dlplane*chubble0**2*np**3*4.0*PI)
      write(*,*)'sdensnorm =',sdensnorm

! Now we just read in the coordinates, and only keep those of type=1
      write(*,*)'#'
      do fn=0,FILES-1
         write(snumber,'(I3)') SNAPSHOT
         write(fnumber,'(I4)') fn
         do i=1,3
            if (snumber(i:i) .eq. ' ') then
               snumber(i:i)='0'
            end if
         end do
         if (FILES.eq.1) then
            filename= dirname(1:len_trim(dirname)) // '/snapshot_' // 
     $           snumber(1:len_trim(snumber))
         else
            filename= dirname(1:len_trim(dirname)) // '/snapshot_' // 
     $           snumber(1:len_trim(snumber)) // '.' // 
     $           fnumber(verify(fnumber,' '):4)
         endif

! now, read in the header
         open(1,file=filename,form='unformatted')
         read(1)npart,massarr,a,redshift,flag_sfr,flag_feedback,nall
     $        ,flag_cooling,FILES,boxsize,omegam,OmegaLambda,HubbleParam
     $        ,unused
         close(1)
         boxsize=boxsize/1000.0 !kpc/h => Mpc/h
                  
! now, read in the snapshot data
         open(1,file=filename,form='unformatted',
     $        access='direct',recl=4)
         nsnap=sum(npart)         
         write(*,*)'### reading Gadget snapshot :',         
     $        filename(1:len_trim(filename))
         write(*,*)'Num of particles in this file =',nsnap
         
         do i=1,nsnap
            read(1,rec=3*i+65)pos(1)
            read(1,rec=3*i+66)pos(2)
            read(1,rec=3*i+67)pos(3)
c            write(*,*)pos
            forall (ii=1:3) pos(ii)=1.e-3*pos(ii) ! kpc => Mpc
!!!!!! place particles onto 8 octant-spheres 
            do kz=1,2
               if (kz.eq.1) then
                  p3=pos(3)
               else
                  p3=pos(3)-boxsize
               endif
               p33=p3*p3
               if (p33.ge.rmax2) cycle
               do ky=1,2
                  if (ky.eq.1) then
                     p2=pos(2)
                  else
                     p2=pos(2)-boxsize
                  endif
                  p23=p2*p2+p33
                  if (p23.ge.rmax2) cycle
                  do kx=1,2
                     if (kx.eq.1) then
                        p1=pos(1)
                     else
                        p1=pos(1)-boxsize
                     endif
                     rr=p1*p1+p23
                     
                     if (rr.ge.rmin2.and.rr.lt.rmax2) then
                        hpvector(1)=p1
                        hpvector(2)=p2
                        hpvector(3)=p3
                        call vec2pix_ring(nside,hpvector,ipring)
                        sdens(ipring)=sdens(ipring)+1.0
                     endif
                  enddo
               enddo
            enddo
         enddo
c         read(1)idamy
         close(1)
      enddo
      write(*,*)'#'
      
!!! setup output filename
      write(pnumber,'(I3)')kplane
      do i=1,3
         if (pnumber(i:i).eq.' ') then
            pnumber(i:i)='0'
         end if
      end do
      pfname=fname(1:len_trim(fname))//
     $     'p'//pnumber(1:len_trim(pnumber))
     $     
      write(*,*)'output file name =',pfname(1:len_trim(pfname))
      write(*,*)'pixpart =',pixpart

!
! write log
!
      outlog=ADJUSTL(ADJUSTR(pfname)//'.log')
      outlog=ADJUSTL(ADJUSTR(outdir)//'/'//outlog)
      open(10,file=outlog,status='unknown')
      write(10,*)'pfname = ',
     $     pfname(1:len_trim(pfname))
      write(10,*)'path = ',path(1:len_trim(path))
      write(10,*)'SNAPSHOT =',SNAPSHOT
      write(10,*)'np       = ',np
      write(10,*)'boxsize  = ',boxsize
      write(10,*)'nres  =',nres
      write(10,*)'nside =',nside
      write(10,*)'npix  =',npix
      write(10,*)'kplane  =',kplane
      write(10,*)'thick   =',thick
      write(10,*)'dlens   =',dlplane
      write(10,*)'alens   =',alens
      write(10,*)'zlens   =',1.0/alens-1.0
      write(10,*)'dcent   =',dcent
      write(10,*)'pixpart =',pixpart
      write(10,*)'omegam =',omegam
      write(10,*)'wde    =',wde
      close(10)

!
!!!!!!!!!!!! spherical harmonics transform --- sdens =>  alm_sdens !!!!!!!!!!!!!!!!
!
      lmax=3*nside
      allocate(hpalm(1:1,0:lmax,0:lmax)) 
      allocate(dw8(1:2*nside,1:1))
      allocate(hpcl(0:lmax,1:1))
      write(*,*)'*** spherical harmonics transformation for ',
     $     pfname(1:len_trim(pfname))
         
!!! particle number density => delta Sigma
      forall (ipix=0:npix-1) sdens(ipix)=sdensnorm*(sdens(ipix)-pixpart)
      outfile=ADJUSTL(ADJUSTR(pfname)//'.spkap.dat')
      outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile)
      open(10,file=outfile,form='unformatted',STATUS='unknown')
      write(10)nside,npix
      write(10)sdens
      close(10)
      
      dw8=1.0_dp
      call map2alm(nside,lmax,lmax,sdens,hpalm,hpz,dw8) 
      deallocate(sdens)
      deallocate(dw8)
      
      outfile=ADJUSTL(ADJUSTR(pfname)//'.spkap.alm')
      outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile)
      open(10,file=outfile,form='unformatted',STATUS='unknown')
      write(10)nside,npix,lmax
      write(10)hpalm
      close(10)

!!! alm_kappa => cl
      call alm2cl(lmax,lmax,hpalm,hpcl)
      outfile=ADJUSTL(ADJUSTR(pfname)//'.spkap.cl')
      outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile)
      open(10,file=outfile,status='unknown')
      do i=1,lmax
         write(10,*)i,hpcl(i,1)
      enddo
      close(10)
      deallocate(hpcl)

!!!!!!!!! alm_kappa => alm_phi (alm_phi=2*alm_kappa/l(l+1)) !!!!!!!!!!!!
      write(*,*)'*** spkap_lm => psi,i, psi,ij: processing...'
!$omp parallel private(l)
!$omp do
      do m=0,lmax
         hpalm(1,0,m)=0.0 !psi_lm=0 for l=0
         do l=1,lmax
            hpalm(1,l,m)=hpalm(1,l,m)/(0.5*l*(l+1)) ! psi_lm = 2 kappa_lm/l/(1+1) for  l not= 0
         enddo
      enddo
!$omp end do
!$omp end parallel

      allocate(hpmap_phi(0:npix-1))
      allocate(hpmap_alp(0:npix-1,1:2))
      allocate(hpmap_mag(0:npix-1,1:3))
       
      call alm2map_der(nside,lmax,lmax,hpalm,hpmap_phi,
     $     hpmap_alp,hpmap_mag)
      deallocate(hpmap_phi)
      deallocate(hpalm)
      
ccc changing hpmap_alp(0:npix-1,1:2) => hpmap(1:2,0:npix-1)
      allocate(hpmap(1:2,0:npix-1))
!$omp parallel private(j)
!$omp do
      do ipix=0,npix-1
         do j=1,2
            hpmap(j,ipix)=hpmap_alp(ipix,j)
         enddo
      enddo
!$omp end do
!$omp end parallel
      deallocate(hpmap_alp)

!!!!!!!! output PSIs
      outfile=ADJUSTL(ADJUSTR(pfname)//'.psis')
      outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile)
      write(*,*)'#'
      write(*,*)'### output'
      write(*,*)'#'
      write(*,*)'output PSIs file = ',outfile(1:len_trim(outfile))
      open(10,file=outfile,status='unknown',form='unformatted')
      write(10)nside,npix
      write(10)hpmap
      deallocate(hpmap)

ccc changing hpmap_mag(0:npix-1,1:3) => hpmap(1:3,0:npix-1)
      allocate(hpmap(1:3,0:npix-1))
!$omp parallel private(j)
!$omp do
      do ipix=0,npix-1
         do j=1,3
            hpmap(j,ipix)=hpmap_mag(ipix,j)
         enddo
      enddo
!$omp end do
!$omp end parallel
      deallocate(hpmap_mag)
      write(10)hpmap
      close(10)

      write(*,*)'*** gadget2psis done'

      end program gadget2psis

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
ccc  Newton-Raphson method
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine newlap(alm,aout,ain)
      implicit real (a-h,o-z)
      implicit integer (i-n)
      parameter (eps=1.0e-5)
      external dldz
      rhs=alm
      a=ain
 356  call qromb(dldz,a,1.0,solv)
      err=abs(solv/rhs-1.0)
c      write(*,*)solv,err
      if (err.gt.eps) then
         dela=-(rhs-solv)/dldz(a)
         a2=a+dela
         if (a2.le.0.0) a2=a+0.5*dela
         a=a2
c         write(*,*)a
         goto 356
      endif
      aout=a
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
ccc  package of solving the integral equation of 
ccc  affine parameter - redshit (scale factor) relation
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      function dldz(a)
      use cosmoparam
      implicit real (a-h,o-z)
      implicit integer (i-n)
      dldz=1.0/sqrt(a*omegam+(1.0-omegam)*a**(1.0-3.0*wde))
      return
      end
ccc
      SUBROUTINE qromb(func,a,b,ss)
      INTEGER JMAX,JMAXP,K,KM
      REAL a,b,func,ss,EPS
      EXTERNAL func
      PARAMETER (EPS=1.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1)
CU    USES polint,trapzd
      INTEGER j
      REAL dss,h(JMAXP),s(JMAXP)
      h(1)=1.
      do 11 j=1,JMAX
        call trapzd(func,a,b,s(j),j)
        if (j.ge.K) then
          call polint(h(j-KM),s(j-KM),K,0.,ss,dss)
          if (abs(dss).le.EPS*abs(ss)) return
        endif
        s(j+1)=s(j)
        h(j+1)=0.25*h(j)
11    continue
      write(*,*)'too many steps in qromb'
      stop
      END
ccc
      SUBROUTINE trapzd(func,a,b,s,n)
      INTEGER n
      REAL a,b,s,func
      EXTERNAL func
      INTEGER it,j
      REAL del,sum,tnm,x
      if (n.eq.1) then
        s=0.5*(b-a)*(func(a)+func(b))
      else
        it=2**(n-2)
        tnm=it
        del=(b-a)/tnm
        x=a+0.5*del
        sum=0.
        do 11 j=1,it
          sum=sum+func(x)
          x=x+del
11      continue
        s=0.5*(s+(b-a)*sum/tnm)
      endif
      return
      END
ccc
      SUBROUTINE polint(xa,ya,n,x,y,dy)
      INTEGER n,NMAX
      REAL dy,x,y,xa(n),ya(n)
      PARAMETER (NMAX=10)
      INTEGER i,m,ns
      REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
      ns=1
      dif=abs(x-xa(1))
      do 11 i=1,n
        dift=abs(x-xa(i))
        if (dift.lt.dif) then
          ns=i
          dif=dift
        endif
        c(i)=ya(i)
        d(i)=ya(i)
11    continue
      y=ya(ns)
      ns=ns-1
      do 13 m=1,n-1
        do 12 i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          if(den.eq.0.)then
             write(*,*)'failure in polint'
             stop
          endif
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
12      continue
        if (2*ns.lt.n-m)then
          dy=c(ns+1)
        else
          dy=d(ns)
          ns=ns-1
        endif
        y=y+dy
13    continue
      return
      END