! gamsky2kap_gauss_quick_hwl16a_yu.f   Time-stamp: <2023-11-13 10:08:24 hamana>
!
! Computing the convergence map on regular grids from shear data (with weight)
! via Kaiser & Squires 1993 inversion formula with Gaussian smoothing kernel 
!
!     History
! 2023-11-13 --- updates on fits keywords "CUNIT1", "CUNIT2", "EQUINOX" and "RADESYS". Special thanks to Y. Utsumi. 
! 2020-10-05 --- a fits keyword "RADESYS" added 
! 2020-03-06 --- quick_hwl16a based on quick_lg_avng_s16a_v2.0 and quick_loc_s16a_v2.0, for making products for data release 
! 2020-02-18 --- quick_lg_avng_s16a_v2.0 based on quick_avng_s16a_v2.0 based on quick_s16a_v2.0_pzpdf, local/global shape noise output added
! 2020-02-06 --- quick_avng_s16a_v2.0 based on quick_s16a_v2.0_pzpdf, mean galaxy number density is used for the normalization
! 2019-09-18 --- quick_s16a_v2.0_pzpdf -- based on quick_v1.4psf and s16a_v2.0_pzfpdfv1.4psf, randomized noise-level computation was replaced with theoretical model 
! 2017-06-07 --- s16a_v2.0_pzfpdf -- real*8 gra,gdec,calc1,calc2,calm,erms 
! 2017-04-20 --- s16a_v2.0_calib -- corrections for calibration factors (m, c1 & c2) implemented, based on v1.4psf
! 2017-04-20 --- v1.4psf -- non-kernel-weighted weighted-galaxy number density map added
! 2016-04-13 --- v1.3psf -- radian/degree unit input option added
! 2015-07-03 --- v1.2psf -- PSF map added
! 2015-07-02 --- v1.2 -- weighted galaxy number density map added
! 2015-06-23 --- v1.13 -- a bag in computation of weigted variance is fixed
! 2015-04-24 --- v1.12 -- weight is taken into account in computation of the kappa and mean galaxy density
! 2014-03-31 --- v1.11 -- bugs in masking scheme fixed
! 2014-03-26 --- v1.1 -- masking scheme improved
! 2014-03-24 --- v1.02
! 2014-03-22 --- v1.01
! 2014-03-18 --- v1.0 -- based on gamsky2kap_gauss_v1.1.f
! 
! Required libraries: cfitsio
!
! Compile: use fortran90 compiler with an openmp option with link to cfitsio, for example,
!   gfortran -O2 -fopenmp gamsky2kap_gauss_quick_hwl16a_yu.f -o gamsky2kap_gauss_quick_hwl16a.exe -L/usr/local/lib -lcfitsio
!
! Execute: 
! gamsky2kap_gauss_quick.exe -i input.dat -tg 1.0 -to 15.0 -f output.fits [options,,,] 
! where input.dat is 10-column data having [ra,dec,e1,e2,c1,c2,em,gweight,esig,psfsize]
! Notice, shear MUST be defined in the celestial (RA-DEC) frame.
! If you do not have an weight for galaxies, give gwright=1 for all galaxies. 
! (c1,c2): additive bias, set (0.0,0.0) if not given
! em: multiplicative bias, set 0.0 if not given
! esig: intrinsic ellipticity, set 1.0 if not given
! psfsize: PSF-FWHM, , set 1.0 if not given
!
! options:   -u : input coordinate unit 'radian' or 'degree' (default = radian)
!
!           -gs : grid size of output maps [arcmin] (default = 0.15) 
!           -tg : Gaussian smoothing scale 'theta_G' [arcmin] (default = 1.0) 
!           -to : The filter truncation scale 'theta_out' [arcmin] (default = 15.0) 
!            -i : input data file name
!            -f : output fits fine name (default = gamsky2kap_gauss_quick.fits)
!           -rc : mask definition radius 'rc' [arcmin] (default = 0.4)
!                 if there is no galaxy within rc, that region is flagged as mask=2
!           -ep : length for the near edge flag (mask=1 region) [arcmin] (default = 1.0)
!           -nt : number of threads of openmp parallel computation (default = 8)
!           -ng : Maximum number of arrays for input galaxies (must be > N_input_galaxies)
!           -sw : if '-sw 1' switch the sign of e2 -> -e2, (default = 0)
!                 *** IMPORTRAN NOTE: lensfit's "e" is defined in the (-RA, Dec) frame ***
!                 ***                 thus, you need switch e2 -> -e2 to transform it  ***
!                 ***                 in the standard celestial (RA, Dec) frame.       ***
!
! ENJOY !
!
! Takashi Hamana
!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      MODULE constants
      real*8 PI,TWOPI
      real*8 deg2rad,rad2am
      END MODULE constants
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      
      program gamsky2kap_gauss_quick
      use constants
      
      implicit real (a-h,o-z)
      implicit integer (i-n)

!!! cfitsio
      integer istatus,iunit,iblocksize,ibitpix,naxis,naxes(3)
      integer i,j,igroup,ifpixel,nelements
      integer pcount,gcount
      logical simple,extend,anyf
      character*88 comment

!!! mass maps 
      integer, allocatable,dimension(:,:)::ns ! number of galaxies in grids
      real*8, allocatable,dimension(:,:)::wns ! sum weight  in grids
      real*8, allocatable,dimension(:,:)::wns2 ! sum weight^2  in grids
      real, allocatable,dimension(:,:)::rec ! reconstructed WL mass
      real, allocatable,dimension(:,:)::bec   ! reconstructed WL mass from 45 deg rotated shear
      real, allocatable,dimension(:,:)::wgq2  ! sum w^2 gamma^2 Q^2
      real, allocatable,dimension(:,:)::wgq2loc  ! sum w^2 gamma^2 Q^2 (locally normalized)
      real, allocatable,dimension(:,:)::snglob    ! rec/sqrt(wgq2)
      real, allocatable,dimension(:,:)::snloc    ! rec(local)/sqrt(wgq2(local))
      real, allocatable,dimension(:,:)::ecal ! additive calibration for E-mode
      real, allocatable,dimension(:,:)::bcal ! additive calibration for B-mode
      real, allocatable,dimension(:,:)::mcal ! multiplicative calibration for B-mode
      real, allocatable,dimension(:,:)::resp ! R (responsibility)
      real, allocatable,dimension(:,:)::gnd   ! source galaxy density
      real, allocatable,dimension(:,:)::ggd   ! gaussian weighted galaxy density = sum_i Gaussian_at_gal_i
      real, allocatable,dimension(:,:)::gwn   ! sum of galaxy weight in grid = sum_i weight_i
      real, allocatable,dimension(:,:)::ggwn  ! sum of Gaussian weighted galaxy weight in grid = sum_i G_i weight_i
      real, allocatable,dimension(:,:)::ggwn2 ! sum of Gaussian weighted galaxy weight^2 in grid = sum_i G_i weight_i^2
      real, allocatable,dimension(:,:)::gps ! gaussian weighted average PSF size 
      real, allocatable,dimension(:,:,:)::ofits ! output fits
      integer, allocatable,dimension(:,:)::mask
      integer, allocatable,dimension(:,:)::mask2
      real, allocatable,dimension(:,:,:)::rand
      real, allocatable,dimension(:,:)::gridsx
      real, allocatable,dimension(:,:)::gridsy

      real, allocatable,dimension(:)::x1,x2,eg1,eg2,psf
      real*8, allocatable,dimension(:)::gra,gdec
      real*8, allocatable,dimension(:)::calc1,calc2,calm,erms,wn
      real*8 ra,dec,c1,c2,em,gweight,esig
      real*8 vare,vare2,wtot,wmean,wnorm,totvar,varkap,totvarloc
      real*8 sumwns,sumwns2,averesp

      character*256 fits,indat
      character*256 opt,arg
      character*6 cunit

      include 'omp_lib.h'

!!! double constants 
      PI=3.141592653589793238462643383279502884197d0
      TWOPI=2.0d0*PI
      deg2rad=PI/180.0d0
      rad2am=180.0d0*60.0d0/PI

!!! single constants
      piinv=1.0/PI
      
!
! set up parameters
!
      gsize=0.15 ! grid size in arcmin
      thetag=1.0 ! in arcmin
      thetaout=15.0 ! in arcmin
      fits='gamsky2kap_gauss_quick.fits'
      rc0=0.4
      edgepad0=1.0 ! length for the near edge flag (arcmin)
      nopenmp=56 ! number of threads
      ncat=8000000 ! Maximum number of arrays for input galaxies (must be > N_input_galaxies)
      iflage2=0
      cunit='radian'
! read options
      narg=iargc()
      do i=1,narg
         call getarg(i,opt)
         call getarg(i+1,arg)
         select case (opt)
         case ('-h')
            write(*,*)'Execute:'
            write(*,*)'gamsky2kap_gauss_quick.exe [options] < input.dat'
            write(*,*)'where input.dat is 5-column data having'
            write(*,*)'(RA, Dec, shear_1, shear_2, weight)'
            write(*,*)'Notice, shear MUST be defined in the celestial'
            write(*,*)'(RA-DEC) frame.'
            write(*,*)'If no weight for galaxies is given,' 
            write(*,*)'assign wright=1 for all galaxies.' 
            write(*,*)''
            write(*,*)'Options:'
            write(*,*)' -u : Input coordinate unit raian or degree'
            write(*,*)'      (default = radian)'
            write(*,*)'-gs : Grid size of output maps [arcmin]'
            write(*,*)'      (default = 0.15)'
            write(*,*)'-tg : Gaussian smoothing scale [arcmin]' 
            write(*,*)'      (default = 1.0)' 
            write(*,*)'-to : The filter truncation scale [arcmin]'
            write(*,*)'      (default = 15.0)' 
            write(*,*)' -f : Output fits file name'
            write(*,*)'      (default = gamsky2kap_gauss_quick.fits)'
            write(*,*)'-rc : Mask definition radius [arcmin]'
            write(*,*)'      (default = 0.4)'
            write(*,*)'      if there is no galaxy within rc,'
            write(*,*)'      that region is flagged as mask=2'
            write(*,*)'-ep : Length for the near edge flag [arcmin]'
            write(*,*)'      (default = 1.0)'
            write(*,*)'       mask=1'
            write(*,*)'-nt : Number of threads for openmp parallel job'
            write(*,*)'      (default = 8)'
            write(*,*)'-ng : number of input galaxies'
            write(*,*)'-sw : If [-sw 1] switch the sign of e2 -> -e2,'
            write(*,*)'      (default = 0)'
            stop
         case ('-u')
            cunit=trim(arg) ! Input coordinate unit raian or degree
         case ('-gs')
            read(arg,*)gsize ! grid size of output maps in arcmin
         case ('-tg')
            read(arg,*)thetag ! theta_G in arcmin
         case ('-to')
            read(arg,*)thetaout ! theta_out in arcmin
         case ('-i')
            indat=arg ! input catalog name
         case ('-f')
            fits=arg ! output fits name
         case ('-rc')
            read(arg,*)rc0 ! mask definition radius
         case ('-ep')
            read(arg,*)edgepad0 ! length for the near edge flag (arcmin)
         case ('-nt')
            read(arg,*)nopenmp ! number of threads
         case ('-ng')
            read(arg,*)ncat ! Maximum number of arrays for input galaxies (must be > N_input_galaxies)
         case ('-sw')
            read(arg,*)iflage2 ! switch the sign of e2 if -sw 1
         end select
      enddo
      
      write(*,'(a31)')'#### input coodinate unit #####'
      if (cunit.eq.'radian') then
         write(*,*)'radian unit input'
      elseif (cunit.eq.'degree') then
         write(*,*)'degree unit input'
      else
         write(*,*)'input coordinate unit undefined'
         write(*,*)'select -u radian (default)'
         write(*,*)'    or -u degree'
         stop
      endif

      allocate(gra(ncat))
      allocate(gdec(ncat))
      allocate(x1(ncat))
      allocate(x2(ncat))
      allocate(eg1(ncat))
      allocate(eg2(ncat))
      allocate(calc1(ncat))
      allocate(calc2(ncat))
      allocate(calm(ncat))
      allocate(erms(ncat))
      allocate(wn(ncat))
      allocate(psf(ncat))

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

      write(*,'(a27)')'#### input parameters #####'
      write(*,*)'      grid size of output maps =',gsize,'[arcmin]'
      write(*,*)'theta_G (Gaussian kernel size) =',thetag,'[arcmin]'
      write(*,*)' theta_out (truncation radius) =',thetaout,'[arcmin]'
      write(*,*)'   rc (mask definition radius) =',rc0,'[arcmin]'
      if (iflage2.eq.1) write(*,*)'switch sign of e2 -> -e2'
      
      deg2grid=60.0/gsize
      gsizedeg=gsize/60.0 ! grid size in deg
      grid2rad=gsizedeg*PI/180.0

      cdelt1=-gsizedeg
      cdelt2=gsizedeg

      intmax=nint(thetaout/gsize+1.0) ! in grid unit

      rmax2=thetaout**2 ! theta_out^2 [arcmin^2]
      s2=1.0/thetag**2 ! 1/theta_G^2 [arcmin^-2]      
      garea=gsize*gsize  ! arcmin^2

!!!!! read data !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      write(*,'(a29)')'#### reading input data #####'

      vare=0.0
      wtot=0.0
      ramin=1.0e8
      ramax=-1.0e8
      decmin=1.0e8
      decmax=-1.0e8
      averesp=0.0
      open(10,file=indat,status='old',form='unformatted')
      do k=1,ncat
         read(10)ra,dec,e1,e2,c1,c2,em,gweight,esig,psfsize
         if (cunit.eq.'radian') then
            gra(k)=ra
            gdec(k)=dec
         else
            gra(k)=ra*deg2rad   ! degree -> radian
            gdec(k)=dec*deg2rad ! degree -> radian
         endif
         psf(k)=psfsize
         eg1(k)=e1
         if (iflage2.eq.1) then
            eg2(k)=-e2
         else
            eg2(k)=e2
         endif
         vare=vare+e1*e1+e2*e2
         calc1(k)=c1
         calc2(k)=c2
         calm(k)=em
         erms(k)=esig*esig
         averesp=averesp+gweight*erms(k)
         wn(k)=gweight
         wtot=wtot+gweight
         if (ra.lt.ramin) ramin=ra
         if (ra.gt.ramax) ramax=ra
         if (dec.lt.decmin) decmin=dec
         if (dec.gt.decmax) decmax=dec
      enddo
      close(10)
      sigmae=sqrt(vare/(1.0*ncat))      
      wmean=wtot/(1.0*ncat)
      wnorm=1.0/wmean     
      averesp=2.0*(1.0-averesp/wtot)
            
!!! weighted variance
      vare2=0.0  
!$omp parallel
!$omp do reduction(+:vare2)
      do i=1,ncat
         wn(i)=wnorm*wn(i)
         vare2=vare2+wn(i)*(eg1(i)*eg1(i)+eg2(i)*eg2(i))
         eg1(i)=eg1(i)*wn(i)
         eg2(i)=eg2(i)*wn(i)
         calc1(i)=calc1(i)*wn(i)
         calc2(i)=calc2(i)*wn(i)
         calm(i)=calm(i)*wn(i)
         erms(i)=erms(i)*wn(i)
      enddo 
!$omp end do
!$omp end parallel
      sigmae2=sqrt(vare2/(1.0*ncat))      
      sigmaerms=sqrt(-0.5*averesp+1.0)

      write(*,*)'num of input galaxies =',ncat
      write(*,*)'weighted sigma e_rms =',sigmaerms
      write(*,*)'          average 2R =',real(averesp)
      write(*,*)'             sigma_e =',sigmae
      write(*,*)'         sigma_gamma =',real(sigmae/averesp)
      write(*,*)'    weighted sigma_e =',sigmae2
      write(*,*)'weighted sigma_gamma =',real(sigmae2/averesp)

!!! setup grids
      if (cunit.eq.'radian') then
         ra0=0.5*(ramin+ramax) ! radian unit
         dec0=0.5*(decmin+decmax) ! radian unit
         crval1=ra0/deg2rad ! degree unit
         crval2=dec0/deg2rad ! degree unit
         ramin=ramin/deg2rad ! radian -> degree unit
         ramax=ramax/deg2rad ! radian -> degree unit
         decmin=decmin/deg2rad ! radian -> degree unit
         decmax=decmax/deg2rad ! radian -> degree unit
      else
         crval1=0.5*(ramin+ramax) ! degree unit
         crval2=0.5*(decmin+decmax) ! degree unit
         ra0=crval1*deg2rad ! radian unit
         dec0=crval2*deg2rad ! radian unit
      endif
      cosdelta0=cos(dec0)
      sindelta0=sin(dec0)
      rad2grid=180.0/(PI*gsizedeg)
!$omp parallel private(sindelta,cosdelta,capa,capf)
!$omp do
      do i=1,ncat
         sindelta=sin(gdec(i))
         cosdelta=cos(gdec(i))
         capa=cosdelta*cos(gra(i)-ra0)
         capf=rad2grid/(sindelta0*sindelta+capa*cosdelta0)
         x1(i)=-capf*cosdelta*sin(gra(i)-ra0)
         x2(i)=capf*(cosdelta0*sindelta-capa*sindelta0)
      enddo
!$omp end do
!$omp end parallel

!
!!! searching for min & max
!
      x1max=-1.0e8
      x2max=-1.0e8
!$omp parallel reduction(max:x1max,x2max)
!$omp do
      do i=1,ncat
         if (x1(i).gt.x1max) x1max=x1(i)
         if (x2(i).gt.x2max) x2max=x2(i)
      enddo
!$omp end do
!$omp end parallel

      x1min=1.0e8
      x2min=1.0e8
!$omp parallel reduction(min:x1min,x2min)
!$omp do
      do i=1,ncat
         if (x1(i).lt.x1min) x1min=x1(i)
         if (x2(i).lt.x2min) x2min=x2(i)
      enddo
!$omp end do
!$omp end parallel

      x1min=x1min-0.5
      x1max=x1max+0.5
      x2min=x2min-0.5
      x2max=x2max+0.5
      nx=int(x1max-x1min+1.0)
      ny=int(x2max-x2min+1.0)
!$omp parallel
!$omp do
      do i=1,ncat
         x1(i)=x1(i)-x1min
         x2(i)=x2(i)-x2min         
      enddo
!$omp end do
!$omp end parallel

      crpix1=-x1min
      crpix2=-x2min
      write(*,'(a21)')'#### setup grid #####'      
      write(*,*)'             R.A. min =',ramin
      write(*,*)'             R.A. max =',ramax
      write(*,*)'             Dec. min =',decmin
      write(*,*)'             Dec. max =',decmax
      write(*,*)'CRVAL1 of output data =',crval1
      write(*,*)'CRVAL2 of output data =',crval2
      write(*,*)'CRPIX1 of output data =',crpix1
      write(*,*)'CRPIX2 of output data =',crpix2      
      write(*,*)'    N1 of output data =',nx
      write(*,*)'    N2 of output data =',ny

      allocate(ns(nx,ny))
      allocate(wns(nx,ny))
      allocate(wns2(nx,ny))
      allocate(rec(nx,ny))
      allocate(bec(nx,ny))
      allocate(wgq2(nx,ny))
      allocate(wgq2loc(nx,ny))
      allocate(snglob(nx,ny))
      allocate(snloc(nx,ny))
      allocate(ecal(nx,ny))
      allocate(bcal(nx,ny))
      allocate(mcal(nx,ny))
      allocate(resp(nx,ny))
      allocate(gnd(nx,ny))
      allocate(ggd(nx,ny))
      allocate(gwn(nx,ny))
      allocate(ggwn(nx,ny))
      allocate(ggwn2(nx,ny))
      allocate(mask(nx,ny))
      allocate(mask2(nx,ny))
      allocate(gps(nx,ny))
      allocate(ofits(nx,ny,10))
      allocate(gridsx(nx,ny))
      allocate(gridsy(nx,ny))

!!! RA, Dec of grid points
!$omp parallel private(dy,dx,capd,capb,
!$omp+coscapb,sincapb,coscapd,sincapd,xx,yy)
!$omp do
      do j=1,ny
         dy=-(1.0*j-crpix2)*grid2rad
         do i=1,nx
           dx=(1.0*i-crpix1)*grid2rad
           capd=atan(sqrt(dx*dx+dy*dy))
           capb=atan2(-dx,dy)
           coscapb=cos(capb)
           sincapb=sin(capb)
           coscapd=cos(capd)
           sincapd=sin(capd)           
           xx=sindelta0*sincapd*coscapb+cosdelta0*coscapd
           yy=sincapd*sincapb
           gridsx(i,j)=ra0+atan2(yy,xx)
           gridsy(i,j)=asin(sindelta0*coscapd-cosdelta0*sincapd*coscapb)
        enddo
      enddo
!$omp end do
!$omp end parallel


!!!! set mask !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      write(*,'(a31)')'#### checking mask regions #####'
! mask=0 for data region
! mask=1 for near edge region
! mask=2 for mask region
      rc=rc0 ! arcmin
      rc=rc/gsize
      irc=int(rc)+1
      rc2=rc*rc
      padding=rc
      edgepad=edgepad0/gsize

!!! masking
!$omp parallel
!$omp do
      do j=1,ny
         do i=1,nx
            mask(i,j)=2
         enddo
      enddo
!$omp end do
!$omp end parallel

!!! Fits's coordinate definition
!!! left-bottom corner of a pixel [i,j] is (i-0.5, j-0.5).
!!! Thus the left-bottom corner of an image is (0.5, 0.5).

!$omp parallel private(iom,iomin,iomax,jom,jomin,jomax,y,yy,x,r2)
!$omp do
      do k=1,ncat
         iom=nint(x1(k))
         iomin=iom-irc
         if (iomin.le.0) iomin=1
         iomax=iom+irc
         if (iomax.gt.nx) iomax=nx
         jom=nint(x2(k))
         jomin=jom-irc
         if (jomin.le.0) jomin=1
         jomax=jom+irc
         if (jomax.gt.ny) jomax=ny
         do j=jomin,jomax
            y=1.0*j ! pixel center is not (i-0.5,j-0.5) but (i,j) 
            yy=(y-x2(k))**2
            do i=iomin,iomax
               if (mask(i,j).eq.0) cycle
               x=1.0*i ! pixel center is not (i-0.5,j-0.5) but (i,j) 
               r2=(x-x1(k))**2+yy
               if (r2.lt.rc2) mask(i,j)=0
            enddo
         enddo
      enddo
!$omp end do
!$omp end parallel

!!! remove isolated data pixel
! 222    222
! 202 => 222
! 222    222
!!! and remove isolated mask pixel
! 000    000
! 020 => 000
! 000    000
!$omp parallel private(j,j1,j2,i,i1,i2,jj,ii)
!$omp do
      do j=1,ny
         j1=j-1
         if (j1.lt.1) j1=1
         j2=j+1
         if (j2.gt.ny) j2=ny
         do 345 i=1,nx
            i1=i-1
            if (i1.lt.1) i1=1
            i2=i+1
            if (i2.gt.nx) i2=nx
            if (mask(i,j).eq.0) then ! remove isolated data pixel
               mask2(i,j)=0
               do jj=j1,j2
                  do ii=i1,i2
                     if (i.ne.ii.and.j.ne.jj.and.mask(ii,jj).eq.0) 
     $                    goto 345
                  enddo
               enddo
               mask2(i,j)=2 ! surrounding 8 pixels are all 2 
            else ! remove isolated mask pixel
               mask2(i,j)=2
            endif
 345     continue
      enddo
!$omp end do
!$omp end parallel

! padding = 2
      rpad=padding
      rpad2=rpad*rpad
      npad=int(rpad+0.9)
!$omp parallel private(j,j1,j2,i,i1,i2,jj,ii,sep)
!$omp do
      do j=1,ny
         j1=j-npad
         if (j1.lt.1) j1=1
         j2=j+npad
         if (j2.gt.ny) j2=ny
         do 444 i=1,nx
            if (mask2(i,j).eq.2) then
               mask(i,j)=2
            else
               mask(i,j)=0
               i1=i-npad
               if (i1.lt.1) i1=1
               i2=i+npad
               if (i2.gt.nx) i2=nx               
               do jj=j1,j2
                  do ii=i1,i2
                     if (mask2(ii,jj).eq.2) then
                        sep=1.0*((ii-i)**2+(jj-j)**2)
                        if (sep.le.rpad2) then
                           mask(i,j)=2
                           goto 444
                        endif
                     endif
                  enddo
               enddo
            endif
 444     continue
      enddo
!$omp end do
!$omp end parallel

!!! remove isolated data pixel
! 222    222
! 202 => 222
! 222    222
!$omp parallel private(j,j1,j2,i,i1,i2,jj,ii)
!$omp do
      do j=1,ny
         j1=j-1
         if (j1.lt.1) j1=1
         j2=j+1
         if (j2.gt.ny) j2=ny
         do 543 i=1,nx            
            if (mask(i,j).eq.0) then ! remove isolated data pixel
               i1=i-1
               if (i1.lt.1) i1=1
               i2=i+1
               if (i2.gt.nx) i2=nx
               do jj=j1,j2
                  do ii=i1,i2
                     if (i.ne.ii.and.j.ne.jj.and.mask(ii,jj).eq.0) 
     $                    goto 543
                  enddo
               enddo
               mask2(i,j)=-1 ! surrounding 8 pixels are all 2 
            endif
 543        continue
      enddo
!$omp end do
!$omp end parallel

!$omp parallel
!$omp workshare  
      forall (i=1:nx,j=1:ny,mask2(i,j).eq.-1) mask(i,j)=2
!$omp end workshare
!$omp end parallel

!!! mark edge flag = 1 
      rpad=edgepad
      rpad2=rpad*rpad
      npad=int(rpad+0.9)
!$omp parallel private(j,j1,j2,i,i1,i2,jj,ii,sep)
!$omp do
      do j=1,ny
         j1=j-npad
         if (j1.lt.1) j1=1
         j2=j+npad
         if (j2.gt.ny) j2=ny
         do 555 i=1,nx
            if (mask(i,j).eq.0) then
               i1=i-npad
               if (i1.lt.1) i1=1
               i2=i+npad
               if (i2.gt.nx) i2=nx
               do jj=j1,j2
                  do ii=i1,i2
                     if (mask(ii,jj).eq.2) then
                        sep=1.0*((ii-i)**2+(jj-j)**2)
                        if (sep.le.rpad2) then
                           mask(i,j)=1
                           goto 555
                        endif
                     endif
                  enddo
               enddo
            endif
 555     continue
      enddo
!$omp end do
!$omp end parallel
      
!!! field edge => mask = 1
!$omp parallel
!$omp workshare      
      forall (i=1:nx,j=1:npad,mask(i,j).eq.0) mask(i,j)=1
!$omp end workshare
!$omp end parallel

!$omp parallel
!$omp workshare  
      forall (i=1:nx,j=ny-npad+1:ny,mask(i,j).eq.0) mask(i,j)=1
!$omp end workshare
!$omp end parallel

!$omp parallel
!$omp workshare  
      forall (i=1:npad,j=1:ny,mask(i,j).eq.0) mask(i,j)=1
!$omp end workshare
!$omp end parallel

!$omp parallel
!$omp workshare  
      forall (i=nx-npad+1:nx,j=1:ny,mask(i,j).eq.0) mask(i,j)=1
!$omp end workshare
!$omp end parallel

!!! effective area
      kmask=0
      lmask=0
      mmask=0
!$omp parallel reduction(+:kmask,lmask,mmask)
!$omp do
      do j=1,ny
         do i=1,nx
            if (mask(i,j).eq.0) then
               kmask=kmask+1
            elseif (mask(i,j).eq.1) then
               lmask=lmask+1
            else
               mmask=mmask+1
            endif
         enddo
      enddo
!$omp end do
!$omp end parallel
      effarea2=garea*kmask/(60.0*60.0) ! clean area in deg^2
      effarea=garea*(kmask+lmask)/(60.0*60.0) ! effective area in deg^2
      
      write(*,*)'mask=0 for data region -------',kmask,'[grids]'
      write(*,*)'mask=1 for near mask region --',lmask,'[grids]'
      write(*,*)'mask=2 for mask region -------',mmask,'[grids]'
      write(*,*)'effective area (mask=0)     =',effarea2,'[sq deg]'
      write(*,*)'effective area (mask=0 & 1) =',effarea,'[sq deg]'

!!!!! compute kappa !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      write(*,'(a42)')'#### reconstructing kappa from shear #####'      
!!! clear arrays
!$omp parallel
!$omp do
      do j=1,ny
         do i=1,nx
            rec(i,j)=0.0 ! clear array
            bec(i,j)=0.0 ! clear array
            wgq2(i,j)=0.0 ! clear array
            wgq2loc(i,j)=0.0 ! clear array
            snglob(i,j)=0.0 ! clear array
            snloc(i,j)=0.0 ! clear array
            ecal(i,j)=0.0 ! clear array
            bcal(i,j)=0.0 ! clear array
            mcal(i,j)=0.0 ! clear array
            resp(i,j)=0.0 ! clear array
            ns(i,j)=0    ! clear array
            wns(i,j)=0.0d0    ! clear array
            wns2(i,j)=0.0d0    ! clear array
            gnd(i,j)=0.0 ! clear array
            ggd(i,j)=0.0 ! clear array
            gwn(i,j)=0.0 ! clear array
            ggwn(i,j)=0.0 ! clear array
            ggwn2(i,j)=0.0 ! clear array
            gps(i,j)=0.0 ! clear array
         enddo
      enddo
!$omp end do
!$omp end parallel

      do k=1,ncat  ! loop for galaxies
         e1=eg1(k)
         e2=eg2(k)
         c1=calc1(k)
         c2=calc2(k)
         em=calm(k)
         esig=erms(k)
         psfsize=psf(k)
         gweight=wn(k)
         re1=e1-c1*averesp
         re2=e2-c2*averesp
         eehalf=0.5*(re1*re1+re2*re2)

         sx=gra(k)
         sy=gdec(k)
         x=x1(k)
         y=x2(k)

         ix=nint(x)
         ixmin=ix-intmax
         if (ixmin.lt.1) ixmin=1
         ixmax=ix+intmax
         if (ixmax.gt.nx) ixmax=nx
      
         iy=nint(y)
         iymin=iy-intmax
         if (iymin.lt.1) iymin=1
         iymax=iy+intmax
         if (iymax.gt.ny) iymax=ny
      
         if (ix.ge.1.and.ix.le.nx.and.
     $        iy.ge.1.and.iy.le.ny) then
            ns(ix,iy)=ns(ix,iy)+1            
            wns(ix,iy)=wns(ix,iy)+gweight
            wns2(ix,iy)=wns2(ix,iy)+gweight*gweight
         endif

!$omp parallel private(sep,gcos,gsin,r2,r2s2,ers,wg,r2inv,d1,d2,de)
!$omp do
         do j=iymin,iymax
            do i=ixmin,ixmax
               call courseangle(sx,sy,gridsx(i,j),gridsy(i,j),
     $              sep,gcos,gsin)
               sep=sep*rad2am
               r2=sep*sep
               if (r2.gt.0.0.and.r2.lt.rmax2) then
                  r2s2=r2*s2
                  ers=exp(-r2s2)
                  wg=1.0-(1.0+r2s2)*ers
                  r2inv=1.0/(r2*PI)
                  d1=(gcos*gcos-gsin*gsin)*r2inv
                  d2=2.0*gcos*gsin*r2inv
                  de=e1*d1+e2*d2
                  rec(i,j)=rec(i,j)+wg*de
                  de=e2*d1-e1*d2
                  bec(i,j)=bec(i,j)+wg*de
                  wgq2(i,j)=wgq2(i,j)+eehalf*(wg*r2inv)**2
c
                  de=c1*d1+c2*d2
                  ecal(i,j)=ecal(i,j)+wg*de
                  de=c2*d1-c1*d2
                  bcal(i,j)=bcal(i,j)+wg*de
                  mcal(i,j)=mcal(i,j)+em
                  resp(i,j)=resp(i,j)+esig
c
                  gnd(i,j)=gnd(i,j)+1.0
                  ggd(i,j)=ggd(i,j)+ers
                  gwn(i,j)=gwn(i,j)+gweight
                  ggwn(i,j)=ggwn(i,j)+ers*gweight
                  ggwn2(i,j)=ggwn2(i,j)+ers*gweight*gweight
                  gps(i,j)=gps(i,j)+ers*psfsize
               endif
            enddo
         enddo
!$omp end do
!$omp end parallel
      enddo

!!! galaxy number
      ngaleff=0
      sumwns=0.0d0
      sumwns2=0.0d0
!$omp parallel
!$omp do reduction(+:ngaleff,sumwns,sumwns2)
      do j=1,ny
         do i=1,nx            
            if (mask(i,j).eq.0) then
               ngaleff=ngaleff+ns(i,j)
               sumwns=sumwns+wns(i,j)
               sumwns2=sumwns2+wns2(i,j)
            endif
         enddo
      enddo
!$omp end do
!$omp end parallel
      gdens=1.0*ngaleff/(effarea2*60.0*60.0)
      wgdens=sumwns*sumwns/(sumwns2*effarea2*60.0*60.0)      
      avgwn=sumwns/(effarea2*60.0*60.0) 

!!! source galaxy number density
!$omp parallel private(parea,pgarea,dj,u2,di,z2)
!$omp do
      do j=1,ny
         do i=1,nx
            parea=0.0
            pgarea=0.0
            do jj=1,ny
               dj=gsize*(jj-j)
               u2=dj*dj
               if (u2.gt.rmax2) cycle
               do ii=1,nx
                  if (mask(ii,jj).ne.2) then
                     di=gsize*(ii-i)
                     z2=u2+di*di
                     if (z2.le.rmax2) then
                        parea=parea+1.0
                        pgarea=pgarea+exp(-z2*s2)
                     endif
                  endif
               enddo
            enddo
            parea=garea*parea
            pgarea=garea*pgarea
            if (parea.gt.0.0) then
               resp(i,j)=1.0-resp(i,j)/gwn(i,j)
               resp(i,j)=0.5/resp(i,j)
               mcal(i,j)=1.0+mcal(i,j)/gwn(i,j)
               mcal(i,j)=1.0/mcal(i,j)
!
               gps(i,j)=gps(i,j)/ggd(i,j)
               gwn(i,j)=gwn(i,j)/parea
               gnd(i,j)=gnd(i,j)/parea
               ggd(i,j)=ggd(i,j)/pgarea    
               ggwn(i,j)=ggwn(i,j)*ggwn(i,j)/(ggwn2(i,j)*pgarea)
            else 
               resp(i,j)=0.5
               mcal(i,j)=1.0
!
               gwn(i,j)=-32768.
               gnd(i,j)=-32768.
               ggd(i,j)=-32768.
               gps(i,j)=-32768.
               ggwn(i,j)=-32768.
            endif
         enddo
      enddo
!$omp end do
!$omp end parallel

!$omp parallel private(delwn)
!$omp do
      do j=1,ny
         do i=1,nx
            if (gwn(i,j).le.0.0) then
               rec(i,j)=-32768.
               bec(i,j)=-32768.
               wgq2(i,j)=-32768.
               wgq2loc(i,j)=-32768.
               snglob(i,j)=-32768.
               snloc(i,j)=-32768.
            else
!!! local estimators
               wgq2loc(i,j)=
     $              wgq2(i,j)*(resp(i,j)*mcal(i,j))**2/gwn(i,j)**2 ! !!! TH20200218 locally defined 
               delwn=-1.0/gwn(i,j) ! the local normalization
               snloc(i,j)=delwn*(rec(i,j)*resp(i,j)-ecal(i,j))*mcal(i,j)
     $              /sqrt(wgq2loc(i,j))
!!! global estimators
               delwn=-1.0/avgwn !!! TH20200206
               wgq2(i,j)=wgq2(i,j)*(resp(i,j)*mcal(i,j))**2/avgwn**2
               rec(i,j)=delwn*(rec(i,j)*resp(i,j)-ecal(i,j))*mcal(i,j)
               bec(i,j)=delwn*(bec(i,j)*resp(i,j)-bcal(i,j))*mcal(i,j)
               snglob(i,j)=rec(i,j)/sqrt(wgq2(i,j))               
            endif
         enddo
      enddo
!$omp end do
!$omp end parallel

!!! sigma_kappa
      varkap=0.0
!$omp parallel
!$omp do reduction(+:varkap)
      do j=1,ny
         do i=1,nx
            if (mask(i,j).eq.0) varkap=varkap+rec(i,j)**2
         enddo
      enddo
!$omp end do
!$omp end parallel
      sigmakap=sqrt(varkap/(1.0*kmask))
      
!!!   global noise RMS
      totvar=0.0
      totvarloc=0.0 ! TH20200218 locally defined 
!$omp parallel
!$omp do reduction(+:totvar,totvarloc)
      do j=1,ny
         do i=1,nx
            if (mask(i,j).eq.0) then
               totvar=totvar+wgq2(i,j)
               totvarloc=totvarloc+wgq2loc(i,j)
            endif
         enddo
      enddo
!$omp end do
!$omp end parallel
      rmsnoise=sqrt(totvar/(1.0*kmask))
      rmsnoiseloc=sqrt(totvarloc/(1.0*kmask))

      write(*,*)'mean galaxy num density =',gdens,'[/arcmin^2]'    
      write(*,*)'weighted galaxy density =',wgdens,'[/arcmin^2]'    
      write(*,*)'            sigma_kappa =',sigmakap
      write(*,*)'       global noise rms =',rmsnoise
      write(*,*)'global noise rms:locally=',rmsnoiseloc

!!! compute local noise RMS
!$omp parallel private(sdev,sdevloc,ekap,bkap,snpixg,snpixl)
!$omp do
      do j=1,ny
         do i=1,nx
            if (mask(i,j).eq.2) then
               sdev=-32768.
               sdevloc=-32768.
               ekap=-32768.
               bkap=-32768.
               snpixg=-32768.
               snpixl=-32768.
               gnd(i,j)=-32768.
               ggd(i,j)=-32768.
               gwn(i,j)=-32768.
               ggwn(i,j)=-32768.
               gps(i,j)=-32768.
               goto 376
            endif
            if (gnd(i,j).le.0.0) then
               sdev=-32768.
               sdevloc=-32768.
               ekap=-32768.           
               bkap=-32768.
               snpixg=-32768.
               snpixl=-32768.
            else
               sdev=sqrt(wgq2(i,j))
               sdevloc=sqrt(wgq2loc(i,j))
               ekap=rec(i,j)/rmsnoise
               bkap=bec(i,j)/rmsnoise
               snpixg=snglob(i,j)
               snpixl=snloc(i,j)
            endif
!            
 376        ofits(i,j,1)=ekap
            ofits(i,j,2)=snpixl            
            ofits(i,j,3)=bkap
            ofits(i,j,4)=sdev
            ofits(i,j,5)=sdevloc
            ofits(i,j,6)=real(wns(i,j))!gwn(i,j)!gnd(i,j)
            ofits(i,j,7)=gwn(i,j)
            ofits(i,j,8)=ggwn(i,j)
            ofits(i,j,9)=1.0*mask(i,j)
            ofits(i,j,10)=gps(i,j)
         enddo
      enddo
!$omp end do
!$omp end parallel     

      write(*,'(a27)')'#### output fits data #####'
      write(*,*)'output fits name = ',trim(adjustl(fits))
ccc fits begin
!
!!! initialize parameters for the output FITS image
!
      simple=.true.
      ibitpix=-32
      naxis=3
      extend=.true.
      istatus=0
      iblocksize=1
      naxes(1)=nx
      naxes(2)=ny
      naxes(3)=10
      nelements=naxes(1)*naxes(2)*naxes(3)
!
      call ftgiou(iunit,istatus)
      call ftopen(iunit,fits,1,iblocksize,istatus)
!      write(*,*)istatus
      if (istatus.eq.0) then 
         call ftdelt(iunit,istatus)
      elseif (istatus.eq.104)then 
         istatus=0
         call ftcmsg
      else !there was some other error opening the file; delete the file anyway
         istatus=0
         call ftcmsg
         call ftdelt(iunit,istatus)
      endif
      call ftinit(iunit,fits,iblocksize,istatus)
      call ftphpr(iunit,simple,ibitpix,naxis,naxes
     $     ,0,1,extend,istatus)
!!!     write the additional header keywords
      call FTPCOM(iunit,'  data[1] = kappa/GLOBAL_SIGMA',istatus)
      call FTPCOM(iunit,'  data[2] = kappa/local_SIGMA',istatus)
      call FTPCOM(iunit,'  data[3] = B-map/GLOBAL_SIGMA',istatus)
      call FTPCOM(iunit,'  data[4] = GLOBAL_SIGMA',istatus)
      call FTPCOM(iunit,'  data[5] = LOCAL_SIGMA',istatus)
      call FTPCOM(iunit,'  data[6] = W-NUM_G [in grid]',istatus)
      call FTPCOM(iunit,'  data[7] = W-N_G(<15) [arcmin^-2]',istatus)
      call FTPCOM(iunit,'  data[8] = Gauss-W-N_G [arcmin^-2]',istatus)
      call FTPCOM(iunit,'  data[9] = MASK',istatus)
      call FTPCOM(iunit,'            MASK = 0 for data',istatus)
      call FTPCOM(iunit,'            MASK = 1 for edge',istatus)
      call FTPCOM(iunit,'            MASK = 2 for no data',istatus)
      call FTPCOM(iunit,'  data[10] = G-PSF [arcsec]',istatus)
      call ftpkye(iunit,'CRPIX1',crpix1,9,'',istatus)
      call ftpkye(iunit,'CRPIX2',crpix2,9,'',istatus)
      call ftpkyf(iunit,'CDELT1',cdelt1,8,'',istatus)
      call ftpkyf(iunit,'CDELT2',cdelt2,8,'',istatus)
      call ftpkye(iunit,'CRVAL1',crval1,9,'',istatus)
      call ftpkye(iunit,'CRVAL2',crval2,9,'',istatus)
      call ftpkyf(iunit,'PC001001',1.0,8,'',istatus)
      call ftpkyf(iunit,'PC001002',0.0,8,'',istatus)
      call ftpkyf(iunit,'PC002001',-0.0,8,'',istatus)
      call ftpkyf(iunit,'PC002002',1.0,8,'',istatus)
      call ftpkyf(iunit,'LONGPOLE',180.0,5,'',istatus)
      call ftpkys(iunit,'CTYPE1','RA---TAN','',istatus)
      call ftpkys(iunit,'CTYPE2','DEC--TAN','',istatus)
      call ftpkys(iunit,'CUNIT1','deg','',istatus)
      call ftpkys(iunit,'CUNIT2','deg','',istatus)
      call ftpkys(iunit,'EQUINOX','2000.0','',istatus)
      call ftpkys(iunit,'RADECSYS','FK5','',istatus)
      call ftpkyf(iunit,'EFF_AREA',effarea2,4,'[deg^2]',istatus)
      call ftpkyf(iunit,'EFFAREA2',effarea,4,'[deg^2]',istatus)
      call ftpkyf(iunit,'GRIDSIZE',gsize,4,'[arcmin]',istatus)
      call ftpkyf(iunit,'THETA_G',thetag,4,'[arcmin]',istatus)
      call ftpkyf(iunit,'THETAOUT',thetaout,4,'[arcmin]',istatus)
      call ftpkyf(iunit,'SIGMA_E',sigmae,8,'',istatus)
      call ftpkyf(iunit,'WSIGMA_E',sigmae2,8,'',istatus)
      sigmag=sigmae/averesp
      sigmag2=sigmae2/averesp
      averespreal=real(averesp)
      call ftpkyf(iunit,'SIGMA_G',sigmag,8,'',istatus)
      call ftpkyf(iunit,'WSIGMA_G',sigmag2,8,'',istatus)   
      call ftpkyf(iunit,'SIG_RMS',sigmaerms,8,'',istatus)
      call ftpkyf(iunit,'MEAN_2R',averespreal,8,'',istatus)
      call ftpkyf(iunit,'GLOBAL_S',rmsnoise,8,'GLOBAL_SIGMA',istatus)
      call ftpkyf(iunit,'LOCAL_S',rmsnoiseloc,8,'LOCAL_SIGMA',istatus)
      call ftpkyf(iunit,'SIGMAKAP',sigmakap,8,'',istatus)
      call ftpkyf(iunit,'GALDENS',gdens,4,'[arcmin^-2]',istatus)
      call ftpkyf(iunit,'WGALDENS',wgdens,4,'[arcmin^-2]',istatus)
      call ftpkyj(iunit,'NGALAXY',ncat,'',istatus)
      call ftpkyf(iunit,'RA_MIN',ramin,8,'',istatus)
      call ftpkyf(iunit,'RA_MAX',ramax,8,'',istatus)
      call ftpkyf(iunit,'Dec_MIN',decmin,8,'',istatus)
      call ftpkyf(iunit,'Dec_MAX',decmax,8,'',istatus)
      call ftpkyf(iunit,'R_MASK',rc0,4,'[arcmin]',istatus)
      call ftpkyf(iunit,'EDGEPAD',edgepad0,4,'[arcmin]',istatus)

      igroup=1
      ifpixel=1
      call ftppre(iunit,igroup,ifpixel,nelements,ofits,istatus)
      call ftclos(iunit, istatus)
      call ftfiou(iunit, istatus)
!!! fits end
      write(*,'(a15)')'#### done #####'

      end program gamsky2kap_gauss_quick

!!!
      subroutine courseangle(a1,d1,a2,d2,sep,cs2,sn2)
      use constants

      implicit real (a-h,o-z)
      implicit integer (i-n)
      
      real*8 alp1,alp2,dta1,dta2
      real*8 theta,phi1,phi2
      real*8 dsep,ctheta,stheta,calp21,salp21
      real*8 cphi1,sphi1,cphi2,sphi2
      real*8 c1s2,s1c2
      
      if (a2.ge.a1) then
         jflag=0
         alp1=a1
         dta1=d1
         alp2=a2
         dta2=d2
      else
         jflag=1
         alp1=a2
         dta1=d2
         alp2=a1
         dta2=d1
      endif

      calp21=dcos(alp2-alp1)
      salp21=dsin(alp2-alp1)
      ctheta=dsin(dta1)*dsin(dta2)+dcos(dta1)*dcos(dta2)*calp21
      stheta=dsqrt(1.0-ctheta*ctheta)
      dsep=datan2(stheta,ctheta)
      sep=abs(dsep)
      if (stheta.eq.0.0d0) then
!         write(*,*)'input points are the identical'
!         stop
         sep=0.0
         cs1=0.0
         sn1=0.0
         cs2=0.0
         sn2=0.0
         return
      endif

      cphi2=salp21/stheta
      cphi1=cphi2*dcos(dta2)
      cphi2=cphi2*dcos(dta1)

      c1s2=dcos(dta1)*dsin(dta2)
      s1c2=dsin(dta1)*dcos(dta2)
      sphi1=(c1s2-s1c2*calp21)/stheta
      sphi2=(c1s2*calp21-s1c2)/stheta

      if (jflag.eq.0) then
         cs1=cphi1
         sn1=sphi1
         cs2=cphi2
         sn2=sphi2
      else
         cs1=cphi2
         sn1=sphi2
         cs2=cphi1
         sn2=sphi1
      endif
      return
      end