! gamsky2kap_gauss_v1.3.f   Time-stamp: <2016-02-10 09:32:13 hamana>
!
! Computing the convergence map on regular grids from shear data (with weight)
! via Kaiser & Squires 1993 inversion formula with Gaussian smoothing kernel 
!
! History 
! 2016-02-10 --- v1.3 -- radian/degree unit input option 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, for example,
!   gfortran -O2 -fopenmp gamsky2kap_gauss_v1.3.f -o gamsky2kap_gauss.exe -L/usr/local/lib -lcfitsio
!
! Execute: 
! gamsky2kap_gauss.exe -tg 1.0 -to 15.0 -r -123456789 -f output.fits [options,,,] < input.dat
! where input.dat is 5-column data having (RA, Dec, shear_1, shear_2, weight)
! Notice, shear MUST be defined in the celestial (RA-DEC) frame.
! If you do not have an weight for galaxies, give wright=1 for all galaxies. 
!
! 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) 
!            -f : output fits fine name (default = gamsky2kap_gauss.fits)
!           -nr : number of random realizations which *MUST BE > 2* (default = 100)
!            -r : random number seed [a muinus 9-digit number] (default = -123456789)
!           -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
      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(:,:)::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(:,:,:)::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(:)::gam1,gam2
      real, allocatable,dimension(:):: stat
      real*8, allocatable,dimension(:):: varan

      real, allocatable,dimension(:)::gra,gdec,x1,x2,eg1,eg2,wn
      real*8 vare,vare2,wtot,wmean,wnorm,totvarran,varkap
      real*8 sumwns,sumwns2

      character*256 fits
      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.fits'
      nran=100
      idum=-123456789
      rc0=0.4
      edgepad0=1.0 ! length for the near edge flag (arcmin)
      nopenmp=16 ! number of threads
      ndatmax=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.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 weight=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.fits)'
            write(*,*)'-nr : Number of random realizations'
            write(*,*)'      which *MUST BE > 2* (default = 100)'
            write(*,*)' -r : Random number seed [a muinus 9-digit]' 
            write(*,*)'      (default = -123456789)'
            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 : Max number of arrays for input galaxies'
            write(*,*)'      (must be > N_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 ('-pix')
            read(arg,*)pixsize ! pizel size (unit of input coordinate) in arcsec
         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 ('-f')
            fits=arg ! output fits name
         case ('-nr')
            read(arg,*)nran  ! number of random realizations > 2         
         case ('-r')
            read(arg,*)idum  ! random number seed in a muinus 9-digit number 
         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,*)ndatmax ! 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

      if (nran.le.1) then
         write(*,*)'*** ERROR ***'
         write(*,*)'number of random realizations should be > 2'
         write(*,*)'reset the option -nr'
         stop
      endif

      allocate(gra(ndatmax))
      allocate(gdec(ndatmax))
      allocate(x1(ndatmax))
      allocate(x2(ndatmax))
      allocate(eg1(ndatmax))
      allocate(eg2(ndatmax))
      allocate(wn(ndatmax))

!!! 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]'
      write(*,*)' number of random realizations =',nran
      write(*,*)'            random number seed =',idum
      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

      allocate(gam1(nran))
      allocate(gam2(nran))
      allocate(stat(nran))
      allocate(varan(nran))

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

      k=0
      vare=0.0
      wtot=0.0
      ramin=1.0e8
      ramax=-1.0e8
      decmin=1.0e8
      decmax=-1.0e8
 111  read(*,*,end=999)ra,dec,e1,e2,gweight
!      gweight=1.0
      ee=e1*e1+e2*e2
      if (ee.ge.1.0) goto 111
      k=k+1
      if (k.gt.ndatmax) then
         write(*,*)'!!! Number of input galaxies exceeds Nmax'
         write(*,*)'!!! set option -ng Nmax so that'
         write(*,*)'!!! Nmax > N_input_galaxy'
         stop
      endif
      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
      eg1(k)=e1
      if (iflage2.eq.1) then
         eg2(k)=-e2
      else
         eg2(k)=e2
      endif
      vare=vare+e1*e1+e2*e2
      wn(k)=gweight
      wtot=wtot+gweight      
      if (ra.lt.ramin) ramin=ra ! input unit 
      if (ra.gt.ramax) ramax=ra ! input unit 
      if (dec.lt.decmin) decmin=dec ! input unit 
      if (dec.gt.decmax) decmax=dec ! input unit
      
      goto 111
 999  ncat=k
      sigmae=sqrt(vare/(1.0*ncat))
      write(*,*)'num of input galaxies =',ncat
      write(*,*)'         sigma_e =',sigmae
      wmean=wtot/(1.0*ncat)
      wnorm=1.0/wmean

!!! 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)         
      enddo 
!$omp end do
!$omp end parallel
      sigmae2=sqrt(vare2/(1.0*ncat))
      write(*,*)'weighted sigma_e =',sigmae2

!!! 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(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(ofits(nx,ny,8))
      allocate(rand(nran,nx,ny))
      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)')'#### cheking mask regions #####'
c mask=0 for data region
c mask=1 for near edge region
c mask=2 for mask region
      rc=rc0 ! arcmin
      rc=rc/gsize
      irc=int(rc)+1
      rc2=rc*rc
      padding=rc
      edgepad=edgepad0/gsize

c 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

ccc Fits's coordinate definition
ccc left-bottom corner of a pixel [i,j] is (i-0.5, j-0.5).
ccc 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

c remove isolated data pixel
c 222    222
c 202 => 222
c 222    222
cc and remove isolated mask pixel
cc 000    000
cc 020 => 000
cc 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
c               do jj=j1,j2
c                  do ii=i1,i2
c                     if (i.ne.ii.and.j.ne.jj.and.mask(ii,jj).eq.2) 
c     $                    goto 345
c                  enddo
c               enddo
c               mask2(i,j)=0 ! surrounding 8 pixels are all 0
            endif
 345     continue
      enddo
!$omp end do
!$omp end parallel

c 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

c remove isolated data pixel
c 222    222
c 202 => 222
c 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

c 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

ccc 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 #####'      
c 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
            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
            do k=1,nran
               rand(k,i,j)=0.0  ! clear array
            enddo
         enddo
      enddo
!$omp end do
!$omp end parallel

      do k=1,ncat  ! loop for galaxies
         e1=eg1(k)
         e2=eg2(k)
         gweight=wn(k)
c generate random data
         do iran=1,nran
            rot=TWOPI*ran2(idum)
            cs=cos(rot)
            sn=sin(rot)
            gam1(iran)=e1*cs-e2*sn
            gam2(iran)=e2*cs+e1*sn
         enddo
         
         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
                  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
                  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
                  do iran=1,nran
                     de=gam1(iran)*d1+gam2(iran)*d2
                     rand(iran,i,j)=rand(iran,i,j)+wg*de
                  enddo
               endif
            enddo
         enddo
!$omp end do
!$omp end parallel
      enddo

ccc 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)      

ccc 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
               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 
               gwn(i,j)=-32768.
               gnd(i,j)=-32768.
               ggd(i,j)=-32768.
               ggwn(i,j)=-32768.
            endif
         enddo
      enddo
!$omp end do
!$omp end parallel

ccc normalization
      fac=-piinv
!$omp parallel private(delwn)
!$omp do
      do j=1,ny
         do i=1,nx
c            if (gnd(i,j).le.0.0) then
            if (gwn(i,j).le.0.0) then
               rec(i,j)=-32768.
               bec(i,j)=-32768.
               do k=1,nran
                  rand(k,i,j)=0.0
               enddo
            else
c               delwn=fac/gnd(i,j)
               delwn=fac/gwn(i,j)
               rec(i,j)=delwn*rec(i,j)
               bec(i,j)=delwn*bec(i,j)
               do k=1,nran
                  rand(k,i,j)=delwn*rand(k,i,j)
               enddo
            endif
         enddo
      enddo
!$omp end do
!$omp end parallel

ccc 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))
      
ccc   global noise RMS
      totvarran=0.0
!$omp parallel
!$omp do reduction(+:totvarran)
      do k=1,nran
         varan(k)=0.0 ! clear array
         do j=1,ny
            do i=1,nx
               if (mask(i,j).eq.0) 
     $              varan(k)=varan(k)+rand(k,i,j)**2
            enddo
         enddo
         varan(k)=varan(k)/(1.0*kmask)
         totvarran=totvarran+varan(k)
      enddo
!$omp end do
!$omp end parallel
      totvarran=totvarran/(1.0*nran)
      rmsnoise=sqrt(totvarran)

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

ccc compute local noise RMS
!$omp parallel private(sdev,ekap,bkap,sn,stat)
!$omp do
      do j=1,ny
         do i=1,nx
            if (mask(i,j).eq.2) then
               sdev=0.0
               ekap=-32768.
               bkap=-32768.
               sn=-32768.
               gnd(i,j)=-32768.
               ggd(i,j)=-32768.
               gwn(i,j)=-32768.
               ggwn(i,j)=-32768.
               goto 376
            endif
            if (gnd(i,j).le.0.0) then
               sdev=0.0
               ekap=-32768.           
               bkap=-32768.
               sn=-32768.
            else
               do k=1,nran
                  stat(k)=rand(k,i,j)
               enddo
               call moment(stat,nran,sdev)
               if (sdev.le.0.0) then
                  sdev=0.0
                  ekap=-32768.           
                  bkap=-32768.
                  sn=-32768.
               else
                  ekap=rec(i,j)/rmsnoise
                  bkap=bec(i,j)/rmsnoise
                  sn=rec(i,j)/sdev
               endif
            endif
c            
 376        ofits(i,j,1)=sn
            ofits(i,j,2)=ekap
            ofits(i,j,3)=bkap
            ofits(i,j,4)=sdev
            ofits(i,j,5)=gwn(i,j)!gnd(i,j)
            ofits(i,j,6)=ggd(i,j)
            ofits(i,j,7)=ggwn(i,j)
            ofits(i,j,8)=1.0*mask(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
c
C initialize parameters for the output FITS image
c
      simple=.true.
      ibitpix=-32
      naxis=3
      extend=.true.
      istatus=0
      iblocksize=1
      naxes(1)=nx
      naxes(2)=ny
      naxes(3)=8
      nelements=naxes(1)*naxes(2)*naxes(3)
c
      call ftgiou(iunit,istatus)
      call ftopen(iunit,fits,1,iblocksize,istatus)
c      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)
C     write the additional header keywords
      call FTPCOM(iunit,'  data[1] = kappa/LOCAL_SIGMA',istatus)
      call FTPCOM(iunit,'  data[2] = kappa/GLOBAL_SIGMA',istatus)
      call FTPCOM(iunit,'  data[3] = B-map/GLOBAL_SIGMA',istatus)
      call FTPCOM(iunit,'  data[4] = LOCAL_SIGMA',istatus)
      call FTPCOM(iunit,'  data[5] = W-GAL_DENS [arcmin^-2]',istatus)
      call FTPCOM(iunit,'  data[6] = G-GAL_DENS [arcmin^-2]',istatus)
      call FTPCOM(iunit,'  data[7] = G-W-GAL_DENS [arcmin^-2]',istatus)
      call FTPCOM(iunit,'  data[8] = 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 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','degree','',istatus)
      call ftpkys(iunit,'CUNIT2','degree','',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)
      call ftpkyf(iunit,'GLOBAL_S',rmsnoise,8,'GLOBAL_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 ftpkyj(iunit,'NRAN',nran,'',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)
ccc fits end
      write(*,'(a15)')'#### done #####'

      end program gamsky2kap_gauss

ccc
      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
c         write(*,*)'input points are the identical'
c         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


      SUBROUTINE moment(data,n,sdev)
      INTEGER n
      REAL ave,sdev,var,data(n)
      INTEGER j
      REAL p,s,ep
      s=0.
      do 11 j=1,n
         s=s+data(j)
 11   continue
      ave=s/n
      var=0.
      ep=0.
      do 12 j=1,n
         s=data(j)-ave
         ep=ep+s
         p=s*s
         var=var+p
 12   continue
      var=(var-ep**2/n)/(n-1)
      sdev=sqrt(var)
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.

      FUNCTION ran2(idum)
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL ran2,AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 41m.