ccc dat2fits.f   Time-stamp: <2016-04-07 15:45:52 hamana>
!
!* required library:
!  Healpix_3.11
!  cfitsio
!
!* compile with gfortran: 
! gfortran dat2fits.f -o dat2fits.exe -O2 -I/usr/local/Healpix_3.11/include -DGFORTRAN -fno-second-underscore -fopenmp -fPIC -L/usr/local/Healpix_3.11/lib -L/usr/local/lib -lhealpix -lhpxgif -lcfitsio -Wl,-R/usr/local/lib
!
!!! usage:
! ./dat2fits.exe [options]
! [options] are:
! -d    : file.dat
! -f    : file.fits
!
!!! example:
!
! ./dat2fits.exe -d allskymap_nres12r014.zs16.mag.dat -f test.fits
!
!!! copyright --- Takashi Hamana 
! 
      program dat2fits
      use alm_tools
      use healpix_types
      use pix_tools
      use udgrade_nr
      use fitstools, only : output_map
      use head_fits 

      implicit real (a-h,o-z)
      implicit integer (i-n)
      
      external dldz
      
!!! variables for fitsio      
      integer(4) :: ft_status,ft_unit,ft_blocksize,ft_bitpix
      integer(4) :: ft_naxis,ft_naxes(2)
      integer(4) :: ft_group,ft_fpixel,ft_nelements
      integer(4) :: ft_pcount,ft_gcount
      character(len=88) :: ft_comment
      logical ft_simple,ft_extend,ft_anyf

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

!!! valiables for Ray-tracing
      real(sp),allocatable :: conv(:),shear1(:),shear2(:)

      character(len=256) :: indat,outfits,fin,filein,outlog

      character*256 opt,arg

!
!!!!!! definitions  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!

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

!
!!!!!!! set up parameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
      narg=iargc()
      indat='allskymap_nres12r014.zs16.mag.dat'
      outfits='test.fits'

      do i=1,narg
         call getarg(i,opt)
         call getarg(i+1,arg)
         select case (opt)
         case ('-d')
            indat=arg
         case ('-f')
            outfits=arg
         end select
      enddo


       
!!!! output kappa & gamma map
      open(10,file=indat,status='old',form='unformatted')
      read(10)nside,npix
      allocate(conv(0:npix-1))
      allocate(shear1(0:npix-1))
      allocate(shear2(0:npix-1))
      read(10)conv
      read(10)shear1
      read(10)shear2
      close(10)
      
!!!! output fits 
      allocate(hpmap(0:npix-1,3))
!$omp parallel
!$omp do
      do ipix=0,npix-1 
         hpmap(ipix,1)=conv(ipix)
         hpmap(ipix,2)=shear1(ipix)
         hpmap(ipix,3)=shear2(ipix)
      enddo
!$omp end do
!$omp end parallel
      deallocate(conv)
      deallocate(shear1)
      deallocate(shear2)

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

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

      call output_map(hpmap,mapheader,outfits)
      deallocate(hpmap)
      
      call ftclos(ft_unit,ft_status)

      end program dat2fits