! gam2kap_gauss_v1.1.f   Time-stamp: <2013-10-29 14:12:03 hamana>
!
! Takashi Hamana
!
! Computing the convergence map on regular grids from shear data (with weight)
! via Kaiser & Squires 1993 inversion formula with Gaussian smoothing kernel 
!
! History 
! gam2kap_gauss_v1.1.f Time-stamp: <2013-10-29 14:11:37 hamana> --- bags fixed
! gam2kap_gauss_v1.0.f Time-stamp: <2013-10-05 14:10:46 hamana> 
! 
! Required libraries: cfitsio
!
! Compile: use fortran90 compiler, for example,
!   gfortran -O2 -fopenmp gam2kap_gauss_v1.0.f -o gam2kap_gauss.exe -L/usr/local/lib -lcfitsio
!   ifort -O2 -openmp gam2kap_gauss_v1.0.f -o gam2kap_gauss.exe -L/usr/local/lib -lcfitsio
!
! Execute: 
! gam2kap_gauss.exe -tg 1.0 -to 15.0 -r -123456789 -f output.fits [options,,,] < input.dat
! where input.dat is 5-column data having (x, y, shear0, shear1, weight) 
! If you do not have weights for galaxies, give wright=1 for all galaxies. 
!
! options: -pix : unit of input coordinates (usually CCD pixel size) [arcsec] (default = 0.202) 
!           -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) 
!           -x0 : CRPIX1 of the reference image (north is up)
!           -y0 : CRPIX2 of the reference image (north is up)
!          -ra0 : CRVAL1 of the reference image (north is up)
!         -dec0 : CRVAL2 of reference image (north is up)
!            -f : output fits fine name (default = wlwmrec_kapgauss.fits)
!           -nr : number of random realizations *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)
!
! ENJOY !!!
!
      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, allocatable,dimension(:,:)::rec ! reconstructed WL mass
      real, allocatable,dimension(:,:)::bec   ! reconstructed WL mass from 45 deg rotated shear
      real, allocatable,dimension(:,:)::gnd   ! galaxy density
      real, allocatable,dimension(:,:,:)::ofits ! output fits
      integer, allocatable,dimension(:,:)::mask
      integer, allocatable,dimension(:,:)::mask2
      real, allocatable,dimension(:,:,:)::rand

      real, allocatable,dimension(:)::gam1,gam2
      real, allocatable,dimension(:):: stat
      real*8, allocatable,dimension(:):: varan

      parameter (ndatmax=1000000) ! Maximum number of input galaxies, change this if needed

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

      character*111 fits
      character*88 opt,arg

      include 'omp_lib.h'

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

      pi=4.0*atan(1.0)
      piinv=1.0/pi
      twopi=2.0*pi

!
! set up parameters
!
      pixsize=0.202
      gsize=0.15
      thetag=1.0
      thetaout=15.0
      crpix1=0.0
      crpix2=0.0
      crval1=0.0
      crval2=0.0
      fits='wlwmrec_kapgauss.fits'
      nran=100
      idum=-123456789
      rc0=0.4
      edgepad0=1.0 ! length for the near edge flag (arcmin)
! read options
      narg=iargc()
      do i=1,narg
         call getarg(i,opt)
         call getarg(i+1,arg)
         select case (opt)
         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 ('-x0')
            read(arg,*)crpix1 ! x-pixel origin for WCS
         case ('-y0')
            read(arg,*)crpix2 ! y-pixel origin for WCS
         case ('-ra0')
            read(arg,*)crval1 ! World coordinate on this axis
         case ('-dec0')
            read(arg,*)crval2 ! World coordinate on this axis
         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)
         end select
      enddo
      
      write(*,'(a31)')'#### input parameters are #####'
      write(*,*)'input pixel zise =',pixsize,'[arcsec]'
      write(*,*)'grid size of output maps =',gsize,'[arcmin]'
      write(*,*)'theta_G (Gaussian kernel size) =',thetag,'[arcmin]'
      write(*,*)'theta_out (truncation radius) =',thetaout,'[arcmin]'
      write(*,*)'CRPIX1 of input data =',crpix1
      write(*,*)'CRPIX2 of input data =',crpix2
      write(*,*)'CRVAL1 of input data =',crval1
      write(*,*)'CRVAL2 of input data =',crval2
      write(*,*)'output fits name = ',trim(adjustl(fits))
      write(*,*)'random number seed =',idum
      
      rintmax=thetaout/thetag ! integration range in sigma (n-sigma)
      deg2grid=60.0/gsize
      pix2grid=pixsize/60.0/gsize

      cdelt1=-gsize/60.0
      cdelt2=gsize/60.0
      xwcs0=crpix1*pixsize/(60.0*gsize)
      ywcs0=crpix2*pixsize/(60.0*gsize)

      tgingrid=thetag/gsize !theta_G in grid size
      intmax=nint(rintmax*tgingrid+1.0) ! in grid unit
      rmax2=(rintmax*tgingrid)**2 ! tout^2 [grid size]^2
      s2=1.0/tgingrid**2 ! 1/theta_G^2 
                         ! (here theta_G is in the grid unit)
      
      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
      xmin=1.0e8
      xmax=-1.0e8
      ymin=1.0e8
      ymax=-1.0e8
 111  read(*,*,end=999)x,y,ex,ey,gweight
      k=k+1
      if (k.gt.ndatmax) then
         write(*,*)'!!! Number of input galaxies exceeds ndatmax'
         write(*,*)'!!! change ndatmax in the source file'
         stop
      endif
      x1(k)=x*pix2grid
      x2(k)=y*pix2grid
      eg1(k)=ex
      eg2(k)=ey
      vare=vare+ex*ex+ey*ey
      wn(k)=gweight
      wtot=wtot+gweight
      if (x1(k).lt.xmin) xmin=x1(k)
      if (x1(k).gt.xmax) xmax=x1(k)
      if (x2(k).lt.ymin) ymin=x2(k)
      if (x2(k).gt.ymax) ymax=x2(k)
      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
      
!$omp parallel
!$omp do
      do i=1,ncat
         x1(i)=x1(i)-xmin
         x2(i)=x2(i)-ymin
         wn(i)=wnorm*wn(i)
         eg1(i)=eg1(i)*wn(i)
         eg2(i)=eg2(i)*wn(i)
      enddo
!$omp end do
!$omp end parallel

      vare2=0.0  
!$omp parallel
!$omp do reduction(+:vare2)
      do i=1,ncat
         vare2=vare2+(eg1(i)*eg1(i)+eg2(i)*eg2(i))
      enddo 
!$omp end do
!$omp end parallel
      sigmae2=sqrt(vare2/(1.0*ncat))
      write(*,*)'weighted sigma_e =',sigmae2

ccc set grids
      write(*,'(a21)')'#### setup grid #####'
      xwcs0=xwcs0-xmin
      ywcs0=ywcs0-ymin
      nx=int(xmax-xmin+1.0)
      ny=int(ymax-ymin+1.0)
      write(*,*)'Nx =',nx
      write(*,*)'Ny =',ny
      write(*,*)'CRPIX1 of output maps =',xwcs0
      write(*,*)'CRPIX2 of output maps =',ywcs0
      write(*,*)'CRVAL1 of output maps =',crval1
      write(*,*)'CRVAL2 of output maps =',crval2
      write(*,*)'CDELT1 of output maps =',cdelt1
      write(*,*)'CDELT2 of output maps =',cdelt2

      allocate(ns(nx,ny))
      allocate(rec(nx,ny))
      allocate(bec(nx,ny))
      allocate(gnd(nx,ny))
      allocate(mask(nx,ny))
      allocate(mask2(nx,ny))
      allocate(ofits(nx,ny,6))
      allocate(rand(nran,nx,ny))

!!!! 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
      gc=0.5*sqrt(2.0)
      rc=rc0 ! arcmin
      rc=rc/gsize
      if (rc.lt.gc) rc=gc
      irc=int(rc)+1
      rc2=rc*rc
      padding=rc
      edgepad=edgepad0/gsize

c masking
      forall (i=1:nx,j=1:ny)
         mask(i,j)=2
      end forall
!$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=j-0.5
            yy=(y-x2(k))**2
            do i=iomin,iomax
               if (mask(i,j).eq.0) cycle
               x=i-0.5
               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 mask
!$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 i=1,nx
            if (mask(i,j).eq.0) then
               mask2(i,j)=0
            else
               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 (mask(ii,jj).eq.2) goto 333
                  enddo
               enddo
               mask2(i,j)=0
               goto 345
 333           mask2(i,j)=2
 345        endif
         enddo
      enddo
!$omp end do
!$omp end parallel

c padding = 2 
      rpad=padding
      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 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=sqrt(1.0*((ii-i)**2+(jj-j)**2))
                        if (sep.le.rpad) then
                           mask(i,j)=2
                           goto 444
                        endif
                     endif
                  enddo
               enddo
 444        endif
         enddo
      enddo
!$omp end do
!$omp end parallel

c mark edge flag = 1 
      rpad=edgepad
      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 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=sqrt(1.0*((ii-i)**2+(jj-j)**2))
                        if (sep.le.rpad) then
                           mask(i,j)=1
                           goto 555
                        endif
                     endif
                  enddo
               enddo
 555        endif
         enddo
      enddo
!$omp end do
!$omp end parallel
      
!!! field edge => mask = 1
      forall (i=1:nx,j=1:npad,mask(i,j).eq.0) mask(i,j)=1
      forall (i=1:nx,j=ny-npad+1:ny,mask(i,j).eq.0) mask(i,j)=1
      forall (i=1:npad,j=1:ny,mask(i,j).eq.0) mask(i,j)=1
      forall (i=nx-npad+1:nx,j=1:ny,mask(i,j).eq.0) mask(i,j)=1

ccc effective area
      kmask=0
      lmask=0
      mmask=0
      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
      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) =',effarea,'[sq deg] '

!!!!! compute kappa !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      write(*,'(a42)')'#### reconstructing kappa from shear #####'
      write(*,*)' number of random realizations =',nran
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
            gnd(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)
c generate random data
ccc v1.1 ccc!$omp parallel private(rot,cs,sn)
ccc v1.1 ccc!$omp do
         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
ccc v1.1 ccc!$omp end do
ccc v1.1 ccc!$omp end parallel

         x=x1(k)+0.5
         y=x2(k)+0.5
      
         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
         endif
      
!$omp parallel private(dy,dy2,dx,dx2,r2,r4,r2s2,ers,wg,d1,d2,de)
!$omp do
         do j=iymin,iymax
            dy=1.0*j-y
            dy2=dy*dy
            do i=ixmin,ixmax
               dx=1.0*i-x
               dx2=dx*dx
               r2=dx2+dy2
               if (r2.gt.0.0.and.r2.lt.rmax2) then
                  r4=1.0/(r2*r2)
                  r2s2=r2*s2
                  ers=exp(-r2s2)
                  wg=1.0-(1.0+r2s2)*ers
                  d1=(dx2-dy2)*r4
                  d2=2.0*dx*dy*r4
                  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                  
                  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
!$omp parallel
!$omp do reduction(+:ngaleff)
      do j=1,ny
         do i=1,nx
            if (mask(i,j).eq.0) ngaleff=ngaleff+ns(i,j)
         enddo
      enddo
!$omp end do
!$omp end parallel
      gdens=1.0*ngaleff/(effarea2*60.0*60.0)
      write(*,*)'mean galaxy num density =',gdens,'[/arcmin^2]'
      
ccc source galaxy number density
!$omp parallel private(parea,dj,u2,di,z2)
!$omp do
      do j=1,ny
         do i=1,nx
            parea=0.0
            do jj=1,ny
               dj=1.0*(jj-j)
               u2=dj*dj
               if (u2.gt.rmax2) cycle
               do ii=1,nx
                  if (mask(ii,jj).ne.2) then
                     di=1.0*(ii-i)
                     z2=u2+di*di
                     if (z2.le.rmax2) parea=parea+1.0
                  endif
               enddo
            enddo
            parea=garea*parea
            if (parea.gt.0.0) then
               gnd(i,j)=gnd(i,j)/parea
            else 
               gnd(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
            if (gnd(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
               delwn=fac/(gnd(i,j)*garea)
               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))
      write(*,*)'sigma_kappa =',sigmakap
     

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(*,*)'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 (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)
               ekap=rec(i,j)/rmsnoise
               bkap=bec(i,j)/rmsnoise
               sn=rec(i,j)/sdev
            endif
            if (mask(i,j).eq.2) then
               sdev=0.0
               ekap=-32768.
               bkap=-32768.
               sn=-32768.
               gnd(i,j)=-32768.
            endif
            ofits(i,j,1)=sn
            ofits(i,j,2)=ekap
            ofits(i,j,3)=bkap
            ofits(i,j,4)=sdev
            ofits(i,j,5)=gnd(i,j)
            ofits(i,j,6)=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 about the input FITS image
c
      simple=.true.
      ibitpix=-32
      naxis=3
      extend=.true.
      istatus=0
      iblocksize=1
      naxes(1)=nx
      naxes(2)=ny
      naxes(3)=6
      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] = GAL_DENS [arcmin^-2]',istatus)
      call FTPCOM(iunit,'  data[6] = 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',xwcs0,9,'',istatus)
      call ftpkye(iunit,'CRPIX2',ywcs0,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',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 ftpkyj(iunit,'NGALAXY',ncat,'',istatus)
      call ftpkyf(iunit,'XMIN',xmin,8,'',istatus)
      call ftpkyf(iunit,'XMAX',xmax,8,'',istatus)
      call ftpkyf(iunit,'YMIN',ymin,8,'',istatus)
      call ftpkyf(iunit,'YMAX',ymax,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

      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.