! 
! Program to read a GREAT08 fits file
!
! Uses cfitsio which can be downloaded from
! http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html
!
! Compile with a fortran 90 compiler 
! 
! Example:  insert the following line in .cshrc
! alias f90io 'ifort -o \!:1.a \!:1.f90  -L/usr/lib -lcfitsio'
!
! To compile:
! f90io read_GREAT08_fits
!
! To run:
! read_GREAT08_fits.a -image $GREAT08DIR/set0001.fit
!
! Author: Catherine Heymans on behalf of GREAT08
! F90 fits code questions?  e-mail heymans@roe.ac.uk
! GREAT08 questions?   e-mail sarah@sarahbridle.net
!
! http://www.great08challenge.info/
! 
!-------------------------------------------------------------

Module Main

!fitsio parameters
Integer, Dimension(1:2)               :: naxes   ! fits image dimension

! fits image
real, dimension(:,:), allocatable     :: image  ! array to contain fits image

! GREAT08 postage stamp properties
integer                               :: x_size,y_size
parameter(x_size = 40, y_size = 40)   ! dimension
integer                             :: Nobj_x, Nobj_y
parameter(Nobj_x = 100, Nobj_y = 100) ! Number

!postage stamp array
real, dimension(1:x_size,1:y_size):: object

end Module Main
!-------------------------------------------------------------

Program read_GREAT08_fits
use main
implicit none

integer:: i,j,ii,jj

! command line input
character*150 arg,opt
integer iargc,narg

! input fits file
character*250 filefits

narg=iargc()

do i=1,narg
   call getarg(i,opt)
   call getarg(i+1,arg)
   select case (opt)
   case ('-h')
      print *,'read_GREAT08_fits.a'
      print *,'Example: '
      print *,'read_GREAT08_fits.a -image set0001.fits'
      print *,''
      STOP

   case ('-image')
      filefits=arg
   end select

enddo

! open the image files
call openfits(filefits)

! loop over image array to extract GREAT08 postage stamps

do i = 1, Nobj_x
   do j = 1, Nobj_y
      
      do ii = 1, x_size
         do jj = 1, y_size

            object(ii,jj) = image((i-1)*x_size + ii, (j-1)*y_size + jj)

         end do
      end do

      ! pixelated centroid/peak finder

      print *, i,j,'Centroid', MAXLOC(object)

      !* Add shape measurement method here, or *
      !* write postage stamps out to file of your choice format *

   end do
end do

deallocate(image)

end Program read_GREAT08_fits

!***************************************************************
! FITSIO ROUTINES 
! Modified from cookbook.f
! http://heasarc.gsfc.nasa.gov/FTP/software/fitsio/c/cookbook.f
!***************************************************************

Subroutine openfits(filename)
use main
implicit none

integer                    :: status,unit
integer                    :: group,naxis
real                       :: nullval
logical anyf
integer, Dimension(1:2)    :: incs

character(len=200) filename
integer             :: readwrite,blocksize
integer             :: nfound
integer, dimension(1:2) :: fpixels,lpixels

! The status and file unit parameters initialized.

status = 0      

! Open the read only FITS file 
! Note blocksize is an obsolete ignorable parameter

readwrite=0

unit = 99
call ftopen(unit,filename,readwrite,blocksize,status)

! Determine the size of the image.

call ftgknj(unit,'NAXIS',1,2,naxes,nfound,status)

! Check that it found both NAXIS1 and NAXIS2 keywords.

if (nfound .ne. 2)then
   write(*,*) filename
   print *,'READIMAGE failed to read the fits file',filename
   return
end if

! Initialize variables

group=1
nullval=-999.0
naxis = 2                 ! no of dimensions
incs(1) = 1
incs(2) = 1               ! sampling interval of fits file

allocate (image(1:naxes(1),1:naxes(2)))

! FITSIO routine

fpixels = 1
lpixels = naxes

call ftgsve(unit,group,naxis,naxes,fpixels,lpixels,incs, &
     nullval,image,anyf,status)

! Check for any error, and if so print out error messages.

if (status .gt. 0)then
   write(*,*) trim(filename),&
        'Openfits error problem, fitsio status = ', status
   call printerror(status)
end if

! Closing FITS file 
      
call ftflus(unit,status)
call ftclos(unit,status)

end Subroutine openfits

!-----------------------------------------------------------------------

subroutine printerror(status)

!  Print out the FITSIO error messages to the user

  integer status
  character errtext*30,errmessage*80

!     check if status is OK (no error); if so, simply return
  if (status .le. 0)return

!     get the text string which describes the error
  call ftgerr(status,errtext)
  print *,'FITSIO Error Status =',status,': ',errtext

!     read and print out all the error messages on the FITSIO stack
  call ftgmsg(errmessage)
  do while (errmessage .ne. ' ')
     print *,errmessage
     call ftgmsg(errmessage)
  end do

end subroutine printerror
!-----------------------------------------------------