ccc graytrix_v2.2.f   Time-stamp: <2017-05-18 11:17:03 hamana>
!
! This is the fortran program of the full-sky gravitational lensing ray-tracing
! in which the light ray path and magnification matrix are evaluated by the 
! multiple lens-plane algorithm following the CALCLENS (Becker 2012, see also Teyssier et al 2009).
! I thank Becker very much for making the source codes of CALCLENS publicity available.
! The lensing magnification matrix is evaluated by using the recurrence relation (e.g., Hilbert et al 2009).
!
!* required library:
!  Healpix (graytrix_v2.2.f was developed and tested with Healpix_3.11)
!  cfitsio
!
!* compile with gfortran: 
! gfortran graytrix_v2.2.f -o graytrix.exe -O2 -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:
! ./graytrix.exe [options] < input_parameter_file
! [options] are:
! -nres : Healpix's parameter "nres"
! -d    : path to output dir (detauld = 'DATA')
! -f    : output file name (default = 'test')
! -om   : Omega_M (default = 0.279)
! -wde  : dark energy EoS parameter "w_DE" (default = -1.0)
! -th   : thickness of shell [Mpc/h] (default = 150.0)
! -ip   : inverse distance weight power for interpolation (default = 1.0)
! -nt   : number of openmp threads (default = 16)
!
! contents of input_parameter_file should be 
! 4 - number of source planes [integer > 0]
! 0.6 - source plane redshift for 1st source plane 
! 1.0 - source plane redshift for 2nd source plane
! 1.4 - source plane redshift for 3rd source plane
! 1.8 - source plane redshift for 4th source plane (more lines follow for more source plane outputs)
! /work2/cfca2048/PSIS/nres12/b450r000p001 - file name for input potential data (b450r000p001.psis) for nearest plane
! /work2/cfca2048/PSIS/nres12/b450r000p002 - file name for input potential data (b450r000p002.psis) for 2nd plane
! /work2/cfca2048/PSIS/nres12/b450r000p003 
! /work2/cfca2048/PSIS/nres12/b900r000p004
! /work2/cfca2048/PSIS/nres12/b900r000p005
! /work2/cfca2048/PSIS/nres12/b900r000p006
! /work2/cfca2048/PSIS/nres12/b1350r000p007
! /work2/cfca2048/PSIS/nres12/b1350r000p008
! /work2/cfca2048/PSIS/nres12/b1350r000p009
! /work2/cfca2048/PSIS/nres12/b1800r000p010
! /work2/cfca2048/PSIS/nres12/b1800r000p011
! /work2/cfca2048/PSIS/nres12/b1800r000p012
! /work2/cfca2048/PSIS/nres12/b2250r000p013
! /work2/cfca2048/PSIS/nres12/b2250r000p014
! /work2/cfca2048/PSIS/nres12/b2250r000p015
! /work2/cfca2048/PSIS/nres12/b2700r000p016
! /work2/cfca2048/PSIS/nres12/b2700r000p017
! /work2/cfca2048/PSIS/nres12/b2700r000p018
! /work2/cfca2048/PSIS/nres12/b3150r000p019
! /work2/cfca2048/PSIS/nres12/b3150r000p020
! /work2/cfca2048/PSIS/nres12/b3150r000p021
! /work2/cfca2048/PSIS/nres12/b3600r000p022
! /work2/cfca2048/PSIS/nres12/b3600r000p023
! /work2/cfca2048/PSIS/nres12/b3600r000p024 - file name for input potential data (b3600r000p024.psis) for farthest plane
!
!!! example:
!
! ./graytrix.exe -nres 12 -f output_file_name < graytrix.in
!
!!! copyright --- Takashi Hamana, Ryuichi Takahashi (2016,2017)
! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      MODULE cosmoparam
      real omegam
      real wde
      END MODULE cosmoparam
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      program graytrix
      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)
      
      external dldz
      
!!! variables for fitsio      
      integer(4) :: ft_status,ft_unit,ft_blocksize,ft_bitpix
      integer(4) :: ft_naxis,ft_naxes(2)
      integer(4) :: ft_group,ft_fpixel,ft_nelements
      integer(4) :: ft_pcount,ft_gcount
      character(len=128) :: ft_comment
      logical :: ft_simple,ft_extend,ft_anyf

!!! valiables for Healpix
      integer(i4b) :: nres,nside,lmax,hporder,hpspin
      integer(i4b) :: nsidein,lmaxin
      integer(i8b) :: npix,ipring,ipnest,ipix,kpix,iray
      integer(i8b) :: npixin
      real(sp),allocatable :: hpmap(:,:) !hpmap(0:npix-1,3) output
      real(sp),allocatable :: hpmap_alp(:,:)!(1:2,0:npix-1) input potentials
      real(sp),allocatable :: hpmap_mag(:,:)!(1:3,0:npix-1) input potentials
      character(len=80), dimension(1:64) :: mapheader

      integer(i8b),allocatable :: neighlist(:,:) ! (8,0:npix-1)
      integer(i8b) :: nlist(8)
      integer(i4b) :: nneigh

!!! valiables for Ray-tracing
      real(sp),allocatable :: ampmat(:,:,:,:) !ampmat(2,2,3,0:npix-1) magnification matrix for three planes
      real(sp),allocatable :: umat(:,:) !(3,0:npix-1) phi,ij at ray positions
      real(sp) :: ampmat1(2,2) ! tmp data

      real(sp),allocatable :: conv(:),shear(:,:),omeg(:) !(0:npix-1,1/2) output
      real(sp),allocatable :: dflang(:,:) !(0:npix-1,1:2) output

      real(dp),allocatable :: theta_i(:),phi_i(:) !(0:npix-1) angular ray position on image-plane
      real(dp),allocatable :: theta_s(:),phi_s(:) !(0:npix-1) output - angular ray position on source-plane
      real(dp) :: theta_ray,phi_ray ! angular position of the ray 
      real(dp),allocatable :: xyzray(:,:) ! 3dim coordinate of ray (x,y,z) = xyzray(3,0:npix-1)
      real(dp) :: dvec(3)
      real(sp) :: xbase(3)
      real(dp),allocatable :: dvecpix(:,:) ! dvec of pixel dvecpix(3,0:npix-1)
      real(sp),allocatable :: betaray(:,:) ! beta of ray betaray(3,0:npix-1)
      real(sp) :: betatmp(3)

      real(sp) :: alphanorm,xxyy,tpnorm
      real(sp),allocatable :: phihat(:,:),thetahat(:,:) ! (3,0:npix-1)
      real(sp),allocatable :: alpxyz(:,:) !(3,0:npix-1))
      real(sp),allocatable :: alpxyzpix(:,:) ! alpha of pixels in xyz alpxyzpix(3,0:npix-1)
      real(sp) :: rotvec(3),rotvecnorm,rotmat(3,3)
      
      real(sp) :: dlambda,b4lambda,c4lambda,dxyzray

      real(dp) :: deltaphi,cosalp,sinalp,alpha,cosdelta,sindelta
      real(dp) :: cosoversin
      real(dp) :: cossep(8),cosray,sinray
      real(dp) :: cosmax(4),wgt(4),wgtnorm,dray
      integer(i4b) :: i4neigh(4)
      
!!! source redshifts
      real,allocatable :: dists(:)
      character(len=128) :: snum

      character(len=128) :: fname,fnamei
      character(len=128) :: outfile,outdir,fin,filein,outlog

      character*128 opt,arg

      include 'omp_lib.h'

!
!!!!!! definitions  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
      chubble0=2.9979e3 ! c/H0 in Mpc/h (H0=100h km/sec/Mpc)

!!! Healpix related
      hporder=1 ! =1 for ring, =2 for nested

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

!!! Healpix related
      nres=13

!!! output file names
      outdir='DATA'
      fname='test'

!!! inverse distance weight power 
!!! f(x) = sum_i w_i f_i(x_i) / sum_i w_i
!!! with w_i = 1/d^pidw
!!! where d is the distance between x and x_i
      pidw=4.0
      
!!! RayTrace related: thickness of shells
      thick=150.0 !Mpc/h

!!! number of threds
      nopenmp=32

      do i=1,narg
         call getarg(i,opt)
         call getarg(i+1,arg)
         select case (opt)
         case ('-nres')
            read(arg,*)nres  ! Healpix parameter
         case ('-d')
            outdir=arg
         case ('-f')
            fname=arg
         case ('-om')
            read(arg,*)omegam
         case ('-wde')
            read(arg,*)wde
         case ('-th')
            read(arg,*)thick
         case ('-ip')
            read(arg,*)pidw
         case ('-nt')
            read(arg,*)nopenmp
         end select
      enddo
      pidw=0.5*pidw ! a numerical trick: d^-pidw = (d*d)^(-0.5*pidw)

!!! set parameters for openmp
!$    call omp_set_num_threads(nopenmp)

!!! Healpix related
      nside=2**nres
      npix=12_i8b*nside*nside
      lmax=3*nside
      write(*,*)'##### Healpix parameters:'
      write(*,*)'# nres  =',nres
      write(*,*)'# nside =',nside
      write(*,*)'# npix  =',npix
      write(*,*)'# lmax  =',lmax

!!! setup neighbours list
      allocate(neighlist(8,0:npix-1))
!$omp parallel private(ipnest,nlist,nneigh,j)
!$omp do
      do ipring=0,npix-1
         do j=1,8
            neighlist(j,ipring)=-1
         enddo
         call ring2nest(nside,ipring,ipnest)
         call neighbours_nest(nside,ipnest,nlist,nneigh)
         do j=1,nneigh
            call nest2ring(nside,nlist(j),neighlist(j,ipring))
         enddo
      enddo
!$omp end do
!$omp end parallel

!!! RayTrace related ::: set initial values
      allocate(ampmat(2,2,3,0:npix-1))
      allocate(umat(3,0:npix-1))
      allocate(theta_i(0:npix-1))
      allocate(phi_i(0:npix-1))
      allocate(betaray(3,0:npix-1))
      allocate(dvecpix(3,0:npix-1))
      allocate(xyzray(3,0:npix-1))
      allocate(phihat(3,0:npix-1))
      allocate(thetahat(3,0:npix-1))
      allocate(alpxyzpix(3,0:npix-1))
      allocate(alpxyz(3,0:npix-1))

!!! log
      outlog=ADJUSTL(ADJUSTR(fname)//'.log')
      outlog=ADJUSTL(ADJUSTR(outdir)//'/'//outlog)
      open(20,file=outlog,status='unknown')
      write(20,*)'outdir = ',outdir(1:len_trim(outdir))
      write(20,*)'fname  = ',fname(1:len_trim(fname))
      write(20,*)'nres  =',nres
      write(20,*)'nside =',nside
      write(20,*)'npix  =',npix      
      write(20,*)'omegam =',omegam
      write(20,*)'wde    =',wde
      write(20,*)'thick =',thick
      write(20,*)'pidw =',pidw
      write(20,*)'source redshifts:'      
! source redshifts
      read(*,*)nsplanes
      write(20,*)'n_zs =',nsplanes      
      allocate(dists(nsplanes))
      do izs=1,nsplanes
         read(*,*)zs
         write(20,*)'zs [',izs,'] =',zs
         azs=1.0/(zs+1.0)
         call qromb(dldz,azs,1.0,ss)
         dists(izs)=ss*chubble0 ! distance to source plane
         write(20,*)'ds [',izs,'] =',dists(izs)
      enddo
      write(20,*)'input files:'

!!! Healpix 
      allocate(hpmap_alp(1:2,0:npix-1))
      allocate(hpmap_mag(1:3,0:npix-1))
      
!
!!! start ray-tracing 
!
      koutflag=0
      isplanes=1
      kplane=0
!
!!!! initial setting
!
c      distl=0.5*thick ! distance to the 1st lens-plane
      distlmin=0.0
      distlmax=thick
      distl=0.75*(distlmax**4-distlmin**4)/(distlmax**3-distlmin**3) ! light-cone weighted mean distance 
                                                                     ! change made in v1.1
!$omp parallel private(dvec,xxyy,tpnorm,j)
!$omp do      
      do ipix=0,npix-1
         call pix2vec_ring(nside,ipix,dvec)
         do j=1,3
            xyzray(j,ipix)=distl*dvec(j)
            betaray(j,ipix)=dvec(j)
            dvecpix(j,ipix)=dvec(j)
         enddo
ccc   theta-phi unit vector
         xxyy=dvec(1)*dvec(1)+dvec(2)*dvec(2)
         tpnorm=sqrt(xxyy)
         thetahat(3,ipix)=-tpnorm
         tpnorm=1.0/tpnorm
         phihat(1,ipix)=-dvec(2)*tpnorm
         phihat(2,ipix)=dvec(1)*tpnorm
         phihat(3,ipix)=0.0
         tpnorm=dvec(3)*tpnorm
         thetahat(1,ipix)=dvec(1)*tpnorm
         thetahat(2,ipix)=dvec(2)*tpnorm
         call vec2ang(dvec,theta_i(ipix),phi_i(ipix)) ! image position
      enddo
!$omp end do
!$omp end parallel

!$omp parallel private(j)
!$omp do 
      do ipix=0,npix-1
         do j=1,2
            ampmat(1,1,j,ipix)=1.0 ! A(1,1)
            ampmat(2,1,j,ipix)=0.0 ! A(2,1)
            ampmat(1,2,j,ipix)=0.0 ! A(1,2)
            ampmat(2,2,j,ipix)=1.0 ! A(2,2)
         enddo
      enddo
!$omp end do
!$omp end parallel

!
!!!!! loop for lens planes
!     
 111  read(*,'(a128)',end=999)fin ! loop for lens planes
      kplane=kplane+1
      write(*,*)'#**********',kplane,'lens-plane **********'      
c      distl=(kplane-0.5)*thick ! distance to this lens-plane
      distlmin=(kplane-1)*thick
      distlmax=kplane*thick
      distl=0.75*(distlmax**4-distlmin**4)/(distlmax**3-distlmin**3) ! change made in v1.1
      distlinv=1.0/distl
      write(*,*)'# distance to this lens-plane =',distl
c      distllast=distl-thick ! distance to the last lens-plane      
      distlmin=(kplane-2)*thick
      distlmax=(kplane-1)*thick
      distllast=0.75*(distlmax**4-distlmin**4)
     $     /(distlmax**3-distlmin**3) ! change made in v1.1
c      distlnext=(kplane+0.5)*thick ! distance to the next lens-plane
      distlmin=kplane*thick
      distlmax=(kplane+1)*thick
      distlnext=0.75*(distlmax**4-distlmin**4)
     $     /(distlmax**3-distlmin**3) ! change made in v1.1
      if (distlnext.ge.dists(isplanes)) then
         koutflag=1
c         distlnext=dists
      endif

!
!!!!!! reading psis file
!
      write(*,*)'##### reading psis data: ',filein(1:len_trim(filein))
      filein=ADJUSTL(ADJUSTR(fin)//'.psis')
      write(20,*)kplane,filein(1:len_trim(filein))
      open(10,file=filein,status='old',form='unformatted')
      read(10)nsidein,npixin
      write(*,*)'# nsidein = ',nsidein
      write(*,*)'# npixin = ',npixin
      read(10)hpmap_alp
      read(10)hpmap_mag
      close(10)

!
!!! covariant derivatives 
!!! see Appendic A of Becker (2013) MNRAS, 435, 115
!
!$omp parallel private(cosoversin)
!$omp do
      do ipix=0,npix-1
         cosoversin=dcos(theta_i(ipix))/dsin(theta_i(ipix))
         hpmap_mag(2,ipix)=hpmap_mag(2,ipix)
     $        -cosoversin*hpmap_alp(2,ipix)
         hpmap_mag(3,ipix)=hpmap_mag(3,ipix)
     $        +cosoversin*hpmap_alp(1,ipix)
      enddo
!$omp end do
!$omp end parallel

ccc alpha in 3-vector
!$omp parallel private(j)
!$omp do
      do ipix=0,npix-1
         do j=1,3
            alpxyzpix(j,ipix)=hpmap_alp(1,ipix)*thetahat(j,ipix)
     $           +hpmap_alp(2,ipix)*phihat(j,ipix)
         enddo
      enddo
!$omp end do
!$omp end parallel
!
!!!!!!! for each ray, searching for 4 nearest pixels from the ray 
!!!!!!! and interpolating the alpha and magnification matrix
!
      write(*,*)'##### interpolating psi,i & psi,ij at ray positions'
!
      if (kplane.eq.1) then ! for the first lens plane (image plane) no interpolation is needed 
!$omp parallel private(j)
!$omp do
         do ipix=0,npix-1 
            umat(1,ipix)=-hpmap_mag(1,ipix) ! umat(1)=phi_11=kappa+gamma1
            umat(2,ipix)=-hpmap_mag(2,ipix) ! umat(2)=phi_12=phi_12=-gamma2
            umat(3,ipix)=-hpmap_mag(3,ipix) ! umat(3)=phi_22=kappa-gamma1
            do j=1,3
               alpxyz(j,ipix)=alpxyzpix(j,ipix)
            enddo
         enddo
!$omp end do
!$omp end parallel

      else ! for i>1 plane, serching for neighbor pixels and doing interpolation 

c clear array
!$omp parallel
!$omp do
         do ipix=0,npix-1
            do j=1,3
               umat(j,ipix)=0.0
               alpxyz(j,ipix)=0.0
            enddo
         enddo
!$omp end do
!$omp end parallel
!
!$omp parallel
!$omp+private(dvec,ipring,i4neigh,theta_ray,phi_ray,cosray,sinray,dray,
!$omp+cosmax,wgt,kpix,nlist,cossep,nneigh,wgtnorm,iray,ampmat1,
!$omp+j,k,l,m,n)
!$omp do
         do ipix=0,npix-1
            do j=1,3 
               dvec(j)=xyzray(j,ipix)*distlinv
            enddo
            call vec2pix_ring(nside,dvec,ipring) !!! pixel nearest to the ray
            i4neigh(1)=ipring   ! nearest pixel
            call vec2ang(dvec,theta_ray,phi_ray) !!! (theta,phi) of ray
            cosray=dcos(theta_ray)
            sinray=dsin(theta_ray)
            cosmax(1)=dcos(theta_i(ipring))*cosray
     $           +dsin(theta_i(ipring))*sinray
     $           *dcos(phi_i(ipring)-phi_ray)
            if (cosmax(1).gt.1.0_dp) then
               dray=0.0_dp
            else
               dray=dacos(cosmax(1))
            endif
            wgt(1)=(dray*dray)**pidw

!!! searching for next nearest 3 pixels            
            cosmax(2)=-1.1_dp  ! searching for 2nd nearest pixel 
            do k=1,7 ! loop for surrounding 8 pixels (1-7)
               kpix=neighlist(k,ipring)
               nlist(k)=kpix
               cossep(k)=dcos(theta_i(kpix))*cosray
     $              +dsin(theta_i(kpix))*sinray
     $              *dcos(phi_i(kpix)-phi_ray)
               if (cossep(k).gt.cosmax(2)) then
                  cosmax(2)=cossep(k)
                  i4neigh(2)=kpix ! 2nd nearest pixel
               endif
            enddo
            nneigh=7
            if (neighlist(8,ipring).ge.0) then ! for the surrounding 8th pixel, if exits
               nneigh=8
               kpix=neighlist(8,ipring)
               nlist(8)=kpix
               cossep(8)=dcos(theta_i(kpix))*cosray
     $              +dsin(theta_i(kpix))*sinray
     $              *dcos(phi_i(kpix)-phi_ray)
               if (cossep(8).gt.cosmax(2)) then
                  cosmax(2)=cossep(8)
                  i4neigh(2)=kpix ! 2nd nearest pixel
               endif
            endif
            if (cosmax(2).gt.1.0_dp) then
               dray=0.0_dp
            else
               dray=dacos(cosmax(2))
            endif
            wgt(2)=(dray*dray)**pidw
c     
            cosmax(3)=-1.1_dp  ! searching for 3rd nearest pixel
            do k=1,nneigh
               if (nlist(k).ne.i4neigh(2).and.
     $              cossep(k).gt.cosmax(3)) then
                  cosmax(3)=cossep(k)
                  i4neigh(3)=nlist(k) ! 3rd nearest pixel
               endif
            enddo
            if (cosmax(3).gt.1.0_dp) then
               dray=0.0_dp
            else
               dray=dacos(cosmax(3))
            endif
            wgt(3)=(dray*dray)**pidw
c     
            cosmax(4)=-1.1_dp ! searching for 4th nearest pixel
            do k=1,nneigh
               if (nlist(k).ne.i4neigh(2).and.nlist(k).ne.i4neigh(3).
     $              and.cossep(k).gt.cosmax(4)) then
                  cosmax(4)=cossep(k)
                  i4neigh(4)=nlist(k) ! 4th nearest pixel
               endif
            enddo
            if (cosmax(4).gt.1.0_dp) then
               dray=0.0_dp
            else
               dray=dacos(cosmax(4))
            endif
            wgt(4)=(dray*dray)**pidw
c     
            wgtnorm=1.0_dp
            do l=2,4
               wgt(l)=wgt(1)/wgt(l)
               wgtnorm=wgtnorm+wgt(l)
            enddo
            wgt(1)=1.0_dp
            do l=1,4
               wgt(l)=wgt(l)/wgtnorm
            enddo
ccc transport the magnification matrix from ray's xbase() to the image position's xbase() 
ccc and computing alpxyz at ray-position from 4 neighbors with the weight wgt
            do l=1,4
               iray=i4neigh(l)
c for magnification matrix
               ampmat1(1,1)=-hpmap_mag(1,iray)
               ampmat1(2,1)=-hpmap_mag(2,iray)
               ampmat1(1,2)=-hpmap_mag(2,iray)
               ampmat1(2,2)=-hpmap_mag(3,iray)
               if (wgt(l).lt.1.0_dp) 
     $              call transampmat(dvecpix(1:3,iray),
     $              dvecpix(1:3,ipix),ampmat1)
               umat(1,ipix)=umat(1,ipix)+wgt(l)*ampmat1(1,1)
               umat(2,ipix)=umat(2,ipix)+wgt(l)*
     $              0.5*(ampmat1(1,2)+ampmat1(2,1))
               umat(3,ipix)=umat(3,ipix)+wgt(l)*ampmat1(2,2)
c for deflection angle in xyz coordinate (thus no base transformation is needed) 
               do j=1,3
                  alpxyz(j,ipix)=alpxyz(j,ipix)+wgt(l)*alpxyzpix(j,iray)
               enddo
            enddo
         enddo
!$omp end do
!$omp end parallel
      endif

      if (koutflag.eq.1) goto 511
!     
!!!!!!!!! computing kappa & gamma on the next place
!!!!!!!!! with multiple-lens plane recurrence relation !!!!
!
 115  distl2distlnext=distlnext/distl
      write(*,*)'##### computing kappa and gamma for the plane =',kplane
      dratinrec2=distl*(distlnext-distllast)
     $     /(distlnext*(distl-distllast))
      dratinrec1=1.0-dratinrec2
      dratinrec3=-(distlnext-distl)/distlnext
!$omp parallel private(m,n)
!$omp do
      do ipix=0,npix-1
c computing A(i,j) on this source plane
         ampmat(1,1,3,ipix)=dratinrec1*ampmat(1,1,1,ipix) ! A11
     $        +dratinrec2*ampmat(1,1,2,ipix)
     $        +dratinrec3*(umat(1,ipix)*ampmat(1,1,2,ipix)
     $        +umat(2,ipix)*ampmat(2,1,2,ipix))
         ampmat(2,1,3,ipix)=dratinrec1*ampmat(2,1,1,ipix) ! A21
     $        +dratinrec2*ampmat(2,1,2,ipix)
     $        +dratinrec3*(umat(2,ipix)*ampmat(1,1,2,ipix)
     $        +umat(3,ipix)*ampmat(2,1,2,ipix))
         ampmat(1,2,3,ipix)=dratinrec1*ampmat(1,2,1,ipix) ! A12
     $        +dratinrec2*ampmat(1,2,2,ipix)
     $        +dratinrec3*(umat(1,ipix)*ampmat(1,2,2,ipix)+
     $        umat(2,ipix)*ampmat(2,2,2,ipix))
         ampmat(2,2,3,ipix)=dratinrec1*ampmat(2,2,1,ipix) ! A22
     $        +dratinrec2*ampmat(2,2,2,ipix)
     $        +dratinrec3*(umat(2,ipix)*ampmat(1,2,2,ipix)+
     $        umat(3,ipix)*ampmat(2,2,2,ipix))
c updating
         do n=1,2
            do m=1,2
               ampmat(m,n,1,ipix)=ampmat(m,n,2,ipix)
               ampmat(m,n,2,ipix)=ampmat(m,n,3,ipix)
            enddo
         enddo
      enddo
!$omp end do
!$omp end parallel
!     
!!!!!!!!! computing ray positions on the next lens plane 
!!!!!!!!! using the ray position and deflection angle at this plane  
!
      write(*,*)'##### ray-position at the next plane'
      c4lambda=distl*distl-distlnext*distlnext
!$omp parallel 
!$omp+private(alphanorm,xbase,rotvec,rotvecnorm,rotmat,betatmp,
!$omp+b4lambda,dlambda,dxyzray,dvec,j,k)
!$omp do
      do ipix=0,npix-1
         alphanorm=sqrt(alpxyz(1,ipix)*alpxyz(1,ipix)
     $        +alpxyz(2,ipix)*alpxyz(2,ipix)
     $        +alpxyz(3,ipix)*alpxyz(3,ipix))
         if (alphanorm.gt.0.0) then
            do j=1,3
               xbase(j)=xyzray(j,ipix)/distl
            enddo
ccc   compute xbase x alpha
            rotvec(1)=xbase(2)*alpxyz(3,ipix)-xbase(3)*alpxyz(2,ipix)
            rotvec(2)=xbase(3)*alpxyz(1,ipix)-xbase(1)*alpxyz(3,ipix)
            rotvec(3)=xbase(1)*alpxyz(2,ipix)-xbase(2)*alpxyz(1,ipix)
            rotvecnorm=1.0/sqrt(rotvec(1)*rotvec(1)
     $           +rotvec(2)*rotvec(2)+rotvec(3)*rotvec(3))
            do j=1,3
               rotvec(j)=rotvec(j)*rotvecnorm
            enddo
ccc   compute rotation matrix
            call gen_rotmat(rotvec,alphanorm,rotmat)
ccc   compute beta for the next plane
            do j=1,3
               betatmp(j)=0.0
               do k=1,3
                  betatmp(j)=betatmp(j)+rotmat(j,k)*betaray(k,ipix)
               enddo
            enddo
         else
            do j=1,3
               betatmp(j)=betaray(j,ipix)  ! RT
            enddo
         endif
ccc   compute lambda
         b4lambda=0.0
         do j=1,3 
            b4lambda=b4lambda+xyzray(j,ipix)*betatmp(j)
         enddo
         dlambda=sqrt(b4lambda*b4lambda-c4lambda)-b4lambda
         if (dlambda.lt.0.0) dlambda=c4lambda/dlambda
         dxyzray=0.0
         do j=1,3
            xyzray(j,ipix)=xyzray(j,ipix)+dlambda*betatmp(j)
            dxyzray=dxyzray+xyzray(j,ipix)*xyzray(j,ipix)
            betaray(j,ipix)=betatmp(j)
         enddo
         dxyzray=1.0/sqrt(dxyzray)
         do j=1,3
            dvec(j)=xyzray(j,ipix)*dxyzray ! normalized position
            xyzray(j,ipix)=dvec(j)*distlnext
         enddo
      enddo
!$omp end do
!$omp end parallel
      goto 111 ! loop for lens planes
      close(20)

!     
!!!!!!!!! computing kappa & gamma on the next place
!!!!!!!!! with multiple-lens plane recurrence relation !!!!
!
 511  write(*,*)'##### output for the source plane =',isplanes
      write(*,*)'# redshift of this source plane = ',zs
      write(*,*)'# distance to this source plane = ',dists(isplanes)
      write(snum,*)isplanes
      fnamei=ADJUSTL(ADJUSTR(fname)//'.zs')
      fnamei=ADJUSTL(ADJUSTR(fnamei)//ADJUSTL(snum))
      write(*,*)'# file name = ',fnamei(1:len_trim(fnamei))
      write(*,*)'# computing kappa & gamma'
      distl2distlnext=dists(isplanes)/distl      
      dratinrec2=distl*(dists(isplanes)-distllast)
     $     /(dists(isplanes)*(distl-distllast))
      dratinrec1=1.0-dratinrec2
      dratinrec3=-(dists(isplanes)-distl)/dists(isplanes)
!$omp parallel private(m,n)
!$omp do
      do ipix=0,npix-1
c computing A(i,j) on this source plane
         ampmat(1,1,3,ipix)=dratinrec1*ampmat(1,1,1,ipix) ! A11
     $        +dratinrec2*ampmat(1,1,2,ipix)
     $        +dratinrec3*(umat(1,ipix)*ampmat(1,1,2,ipix)
     $        +umat(2,ipix)*ampmat(2,1,2,ipix))
         ampmat(2,1,3,ipix)=dratinrec1*ampmat(2,1,1,ipix) ! A21
     $        +dratinrec2*ampmat(2,1,2,ipix)
     $        +dratinrec3*(umat(2,ipix)*ampmat(1,1,2,ipix)
     $        +umat(3,ipix)*ampmat(2,1,2,ipix))
         ampmat(1,2,3,ipix)=dratinrec1*ampmat(1,2,1,ipix) ! A12
     $        +dratinrec2*ampmat(1,2,2,ipix)
     $        +dratinrec3*(umat(1,ipix)*ampmat(1,2,2,ipix)+
     $        umat(2,ipix)*ampmat(2,2,2,ipix))
         ampmat(2,2,3,ipix)=dratinrec1*ampmat(2,2,1,ipix) ! A22
     $        +dratinrec2*ampmat(2,2,2,ipix)
     $        +dratinrec3*(umat(2,ipix)*ampmat(1,2,2,ipix)+
     $        umat(3,ipix)*ampmat(2,2,2,ipix))
      enddo
!$omp end do
!$omp end parallel

cccc output
      allocate(conv(0:npix-1))
      allocate(shear(0:npix-1,1:2))
      allocate(omeg(0:npix-1))
!$omp parallel
!$omp do
      do ipix=0,npix-1
         conv(ipix)=1.0-0.5*(ampmat(1,1,3,ipix)+ampmat(2,2,3,ipix))
         shear(ipix,1)=0.5*(ampmat(2,2,3,ipix)-ampmat(1,1,3,ipix))
         shear(ipix,2)=-0.5*(ampmat(1,2,3,ipix)+ampmat(2,1,3,ipix))
         omeg(ipix)=0.5*(ampmat(2,1,3,ipix)-ampmat(1,2,3,ipix))
      enddo
!$omp end do
!$omp end parallel
       
!!!! output kappa & gamma map
      outfile=ADJUSTL(ADJUSTR(fnamei)//'.mag.dat')
      outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile)
      open(10,file=outfile,status='unknown',form='unformatted')
      write(10)nside,npix
      write(10)conv
      write(10)(shear(ipix,1),ipix=0,npix-1)
      write(10)(shear(ipix,2),ipix=0,npix-1)
      write(10)omeg
      close(10)
      deallocate(omeg)
      
!!!! output fits 
      allocate(hpmap(0:npix-1,3))
!$omp parallel
!$omp do
      do ipix=0,npix-1 
         hpmap(ipix,1)=conv(ipix)
         hpmap(ipix,2)=shear(ipix,1)
         hpmap(ipix,3)=shear(ipix,2)
      enddo
!$omp end do
!$omp end parallel
      deallocate(conv)
      deallocate(shear)      

!!! check if fitsfile exists, and if it does, delete it 
      outfile=ADJUSTL(ADJUSTR(fnamei)//'.mag.fits')
      outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile)
      ft_status=0
      ft_blocksize=1
      call ftgiou(ft_unit,ft_status)
      call ftopen(ft_unit,outfile,1,ft_blocksize,ft_status)
      if (ft_status.eq.0) then 
         call ftdelt(ft_unit,ft_status)
      elseif (ft_status.eq.104)then 
         ft_status=0
         call ftcmsg
      else ! there was some other error opening the file; delete the file anyway
         ft_status=0
         call ftcmsg
         call ftdelt(ft_unit,ft_status)
      endif
      call ftclos(ft_unit,ft_status)
      call ftfiou(ft_unit,ft_status)

!!! output fits
      call write_minimal_header(mapheader, 'MAP', nside=nside,
     $     order=hporder)
      call add_card(mapheader,'TTYPE1','kappa','label for field 1')
      call add_card(mapheader,'TTYPE2','shear1','label for field 2')
      call add_card(mapheader,'TTYPE3','shear2','label for field 3')

      call output_map(hpmap,mapheader,outfile)
      deallocate(hpmap)

!     
!!!!!!!!! computing ray positions on the next lens plane 
!!!!!!!!! using the ray position and deflection angle at this plane  
!
      write(*,*)'##### ray-position at the source plane'
      allocate(theta_s(0:npix-1))
      allocate(phi_s(0:npix-1))
      c4lambda=distl*distl-dists(isplanes)*dists(isplanes)
!$omp parallel 
!$omp+private(alphanorm,xbase,rotvec,rotvecnorm,rotmat,betatmp,
!$omp+b4lambda,dlambda,dxyzray,dvec,j,k)
!$omp do
      do ipix=0,npix-1
         alphanorm=sqrt(alpxyz(1,ipix)*alpxyz(1,ipix)
     $        +alpxyz(2,ipix)*alpxyz(2,ipix)
     $        +alpxyz(3,ipix)*alpxyz(3,ipix))
         if (alphanorm.gt.0.0) then
            do j=1,3
               xbase(j)=xyzray(j,ipix)/distl
            enddo
ccc   compute xbase x alpha
            rotvec(1)=xbase(2)*alpxyz(3,ipix)-xbase(3)*alpxyz(2,ipix)
            rotvec(2)=xbase(3)*alpxyz(1,ipix)-xbase(1)*alpxyz(3,ipix)
            rotvec(3)=xbase(1)*alpxyz(2,ipix)-xbase(2)*alpxyz(1,ipix)
            rotvecnorm=1.0/sqrt(rotvec(1)*rotvec(1)
     $           +rotvec(2)*rotvec(2)+rotvec(3)*rotvec(3))
            do j=1,3
               rotvec(j)=rotvec(j)*rotvecnorm
            enddo
ccc   compute rotation matrix
            call gen_rotmat(rotvec,alphanorm,rotmat)
ccc   compute beta for the next plane
            do j=1,3
               betatmp(j)=0.0
               do k=1,3
                  betatmp(j)=betatmp(j)+rotmat(j,k)*betaray(k,ipix)
               enddo
            enddo
         else
            do j=1,3
               betatmp(j)=betaray(j,ipix)  ! RT
            enddo
         endif
ccc   compute lambda
         b4lambda=0.0
         do j=1,3 
            b4lambda=b4lambda+xyzray(j,ipix)*betatmp(j)
         enddo
         dlambda=sqrt(b4lambda*b4lambda-c4lambda)-b4lambda
         if (dlambda.lt.0.0) dlambda=c4lambda/dlambda
         dxyzray=0.0
         do j=1,3
            dvec(j)=xyzray(j,ipix)+dlambda*betatmp(j)
            dxyzray=dxyzray+dvec(j)*dvec(j)
         enddo
         dxyzray=1.0/sqrt(dxyzray)
         do j=1,3
            dvec(j)=dvec(j)*dxyzray ! normalized position
         enddo
         call vec2ang(dvec,theta_s(ipix),phi_s(ipix))
      enddo
!$omp end do
!$omp end parallel

!!!! output source position
      outfile=ADJUSTL(ADJUSTR(fnamei)//'.src.dat')
      outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile)
      open(10,file=outfile,status='unknown',form='unformatted')
      write(10)nside,npix
      write(10)(theta_s(ipix),ipix=0,npix-1)   ! RT
      write(10)(phi_s(ipix),ipix=0,npix-1)
      close(10)
      deallocate(theta_s)
      deallocate(phi_s)
      
      koutflag=0
      isplanes=isplanes+1
      if (isplanes.gt.nsplanes) goto 999
      goto 115 ! return to the loop for the lens planes

c!!!! compute deflection angle
c      allocate(dflang(0:npix-1,1:2))
c!$omp parallel private(deltaphi,cosalp,alpha,sinalp,cosdelta,sindelta)
c!$omp do
c      do ipix=0,npix-1
c         deltaphi=phi_s(ipix)-phi_i(ipix)
c         cosalp=dcos(theta_i(ipix))*dcos(theta_s(ipix))
c     $        +dsin(theta_i(ipix))*dsin(theta_s(ipix))*dcos(deltaphi)
c         if (cosalp.gt.1.0_dp) then
c            alpha=0.0_dp
c         else
c            alpha=dacos(cosalp)
c         endif
c         sinalp=dsin(alpha)
c         if (sinalp.eq.0.0_dp) then
c            dflang(ipix,1)=0.0
c            dflang(ipix,2)=0.0
c         else
c            cosdelta=(dcos(theta_s(ipix))-dcos(theta_i(ipix))*cosalp)
c     $           /(dsin(theta_i(ipix))*sinalp)
c            sindelta=dsin(deltaphi)*dsin(theta_s(ipix))/sinalp
c            dflang(ipix,1)=-alpha*cosdelta
c            dflang(ipix,2)=alpha*sindelta
c         endif
c      enddo
c!$omp end do
c!$omp end parallel
c!!!! output deflection angle
c      outfile=ADJUSTL(ADJUSTR(fname)//'.alp.dat')
c      outfile=ADJUSTL(ADJUSTR(outdir)//'/'//outfile)
c      open(10,file=outfile,status='unknown',form='unformatted')
c      write(10)nside,npix
c      write(10)(dflang(ipix,1),ipix=0,npix-1)
c      write(10)(dflang(ipix,2),ipix=0,npix-1)
c      close(10)


!!! ray-tracing done, closing arrays
 999  deallocate(neighlist)
      deallocate(betaray)
      deallocate(dvecpix)
      deallocate(phihat)
      deallocate(thetahat)
      deallocate(alpxyzpix)
      deallocate(alpxyz)
      deallocate(umat)
      deallocate(hpmap_alp)
      deallocate(hpmap_mag)
      deallocate(xyzray)
      deallocate(theta_i)
      deallocate(phi_i)

      end program graytrix

ccc
      subroutine transampmat(xbase,dvec,ampmat)
ccc This is the FORTRAN version of the paratrans_tangtensor of CALCLENS (Becker 2012).
ccc I thank Becker very much for making the source codes of CALCLENS publicity available.
      use healpix_types
      implicit real (a-h,o-z)
      implicit integer (i-n)
      real(dp) :: xbase(3),dvec(3) !!!! should be normalized !!!
      real(sp) :: ampmat(2,2),t1(2,2)
      real(dp) :: axis(3),cosangle,sinangle,p(3)
      real(dp) :: ephi_dvec(3),etheta_dvec(3),rephi_xbase(3)
      real(dp) :: norm,sinpsi,cospsi,r(2,2),rt(2,2)
      
      axis(1)=xbase(2)*dvec(3)-xbase(3)*dvec(2)
      axis(2)=xbase(3)*dvec(1)-xbase(1)*dvec(3)
      axis(3)=xbase(1)*dvec(2)-xbase(2)*dvec(1)
      cosangle=xbase(1)*dvec(1)+xbase(2)*dvec(2)+xbase(3)*dvec(3)
      sinangle=dsqrt(axis(1)*axis(1)+axis(2)*axis(2)+axis(3)*axis(3))
      if (sinangle.ne.0.0_dp) then
         forall (j=1:3) axis(j)=axis(j)/sinangle
      else
         axis(1)=1.0_dp
         axis(2)=0.0_dp
         axis(3)=0.0_dp
      endif

      p(1)=-xbase(2)
      p(2)=xbase(1)
      p(3)=0.0_dp
      call gen_crotmax(p,rephi_xbase,axis,cosangle,sinangle)
      
      ephi_dvec(1)=-dvec(2)
      ephi_dvec(2)=dvec(1)
      ephi_dvec(3)=0.0_dp
      
      etheta_dvec(1)=dvec(3)*dvec(1)
      etheta_dvec(2)=dvec(3)*dvec(2)
      etheta_dvec(3)=-1.0_dp*(dvec(1)*dvec(1)+dvec(2)*dvec(2))

      if (dvec(3).gt.1.0_dp) dvec(3)=1.0_dp
      if (xbase(3).gt.1.0_dp) xbase(3)=1.0_dp
      norm=dsqrt((1.0_dp-dvec(3)*dvec(3))*(1.0_dp-xbase(3)*xbase(3)))

      sinpsi=(rephi_xbase(1)*etheta_dvec(1)
     $     +rephi_xbase(2)*etheta_dvec(2)
     $     +rephi_xbase(3)*etheta_dvec(3))/norm
      cospsi=(rephi_xbase(1)*ephi_dvec(1)
     $     +rephi_xbase(2)*ephi_dvec(2)
     $     +rephi_xbase(3)*ephi_dvec(3))/norm
      
      r(1,1)=cospsi
      r(1,2)=-1.0_dp*sinpsi
      r(2,1)=sinpsi
      r(2,2)=cospsi
      rt(1,1)=r(1,1)
      rt(1,2)=r(2,1)
      rt(2,1)=r(1,2)
      rt(2,2)=r(2,2)
      
      do i=1,2
         do j=1,2
            t1(i,j)=ampmat(i,1)*r(1,j)+ampmat(i,2)*r(2,j)
         enddo
      enddo
      do i=1,2
         do j=1,2
            ampmat(i,j)=rt(i,1)*t1(1,j)+rt(i,2)*t1(2,j)
         enddo
      enddo

      return
      end

      subroutine gen_crotmax(vec,rvec,axis,cosangle,sinangle)
      use healpix_types
      implicit real (a-h,o-z)
      implicit integer (i-n)
      real(dp) :: vec(3),rvec(3),axis(3),cosangle,sinangle
      real(dp) :: axisdotvec,axiscrossvec(3)
      axisdotvec=axis(1)*vec(1)+axis(2)*vec(2)+axis(3)*vec(3)
      axiscrossvec(1)=axis(2)*vec(3)-axis(3)*vec(2)
      axiscrossvec(2)=axis(3)*vec(1)-axis(1)*vec(3)
      axiscrossvec(3)=axis(1)*vec(2)-axis(2)*vec(1)
      rvec(1)=vec(1)*cosangle+axis(1)*axisdotvec*(1.0_dp-cosangle)
     $     +axiscrossvec(1)*sinangle
      rvec(2)=vec(2)*cosangle+axis(2)*axisdotvec*(1.0_dp-cosangle)
     $     +axiscrossvec(2)*sinangle
      rvec(3)=vec(3)*cosangle+axis(3)*axisdotvec*(1.0_dp-cosangle)
     $     +axiscrossvec(3)*sinangle
      return
      end

ccc   
      SUBROUTINE gen_rotmat(v,t,r)
      use healpix_types
      implicit real (a-h,o-z)
      implicit integer (i-n)
      real(sp) :: v(3),t,r(3,3)
      real(sp) :: c,s,c1
      real(sp) :: v1,v2,v3,s1,s2,s3
      c=cos(t)
      s=sin(t)
      c1=1.0-c
      v1=v(1)*c1
      v2=v(2)*c1
      v3=v(3)*c1
      s1=v(1)*s
      s2=v(2)*s
      s3=v(3)*s
      r(1,1)=v(1)*v1+c
      r(2,2)=v(2)*v2+c
      r(3,3)=v(3)*v3+c
      r(1,2)=v(1)*v2-s3
      r(2,1)=v(2)*v1+s3
      r(1,3)=v(1)*v3+s2
      r(3,1)=v(3)*v1-s2
      r(2,3)=v(2)*v3-s1
      r(3,2)=v(3)*v2+s1
      return
      END


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)
      if (err.gt.eps) then
         dela=-(rhs-solv)/dldz(a)
         a=a+dela
         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