#
# QDPHOT: Quick & Dirty Photometry
#
# 1999 Copyleft (L) Kenneth J. Mighell (Kitt Peak National Observatory)
#
# 19990913: QDPHOT_V1.0
# 19990914: QDPHOT_V1.1
# 19990920: QDPHOT_V1.2
#
# task qdphot = t_qdphot
#
#

include	<imhdr.h>
include <imio.h>
include <mach.h>

procedure t_qdphot()

pointer	sp

# 1st image
pointer im1                     # descriptor
pointer im1fn 			# filename string pointer
pointer im1sect			# section
bool    im1all                  # entire image?
int     im1pixtype		# pixel type
int     im1nx			# number of columns
int     im1ny			# number of rows
int     im1dx                   # column offset
int     im1dy                   # row offset
pointer im1ptr			# pixel pointer
real    im1gain	                # gain in electrons per adu

# 2nd image
pointer im2                     # descriptor
pointer im2fn 			# filename string pointer
pointer im2sect			# section
bool    im2all                  # entire image?
int     im2pixtype		# pixel type
int     im2nx			# number of columns
int     im2ny			# number of rows
int     im2dx                   # column offset
int     im2dy                   # row offset
pointer im2ptr			# pixel pointer
real    im2gain	                # gain in electrons per adu

real    mnpk1
real    mxpk1

bool    debug
int     status
bool    ok

pointer outptr			# results pointer
int     mxnstr
int     mxncol
int     szoutv
int     nstr

bool    both

pointer out1
int	out1fd

pointer out2
int	out2fd

int     ofd
int     j

# functions
real    clgetr()
int     open()

define	bye_ 91

begin

   debug = true

   call smark (sp)
   call salloc (im1fn, SZ_LINE, TY_CHAR)
   call salloc (im1sect, SZ_FNAME, TY_CHAR)
   call salloc (im2fn, SZ_LINE, TY_CHAR)
   call salloc (im2sect, SZ_FNAME, TY_CHAR)
   call salloc (out1, SZ_FNAME, TY_CHAR)
   call salloc (out2, SZ_FNAME, TY_CHAR)

   # get 1st image
   call qdgtim ("image1", im1, im1fn, im1sect, im1all, im1pixtype, 
      im1nx, im1ny, im1dx, im1dy, im1ptr, status)
   if (status != 0) goto bye_
   # get gain of 1st image
   im1gain = clgetr("gain1")
   # convert 1st image from ADU to electrons
   call amulkr(Memr[im1ptr], im1gain, Memr[im1ptr], im1nx*im1ny)

   # Get the minimum peak
   mnpk1 = clgetr("mnpeak1")

   # Get the minimum peak
   mxpk1 = clgetr("mxpeak1")

   # get 2nd image
   call qdgtim ("image2", im2, im2fn, im2sect, im2all, im2pixtype, 
      im2nx, im2ny, im2dx, im2dy, im2ptr, status)
   if (status != 0) goto bye_
   # get gain of 2nd image
   im2gain = clgetr("gain2")
   # convert 2nd image from ADU to electrons
   call amulkr(Memr[im2ptr], im2gain, Memr[im2ptr], im2nx*im2ny)

   ok = (im1nx==im2nx) && (im1ny==im2ny) 
   if (!ok) {
      call printf ("\nQDPHOT: *** ERROR *** 8=X\n")
      call printf (  "QDPHOT: Image section sizes must be identical!\n\n")
      goto bye_
   }

   # Get 1st output filename
   call clgstr("out1", Memc[out1], SZ_FNAME)
   # Open it
   if (Memc[out1] != EOS) {
      out1fd = open (Memc[out1], NEW_FILE, TEXT_FILE)
   } else {
      out1fd = STDOUT
   }

   # Get 2nd output filename
   call clgstr("out2", Memc[out2], SZ_FNAME)
   # Open it
   both = false
   if (Memc[out2] != EOS) {
      out2fd = open (Memc[out2], NEW_FILE, TEXT_FILE)
      if (out2fd == out1fd) both = true
   } else {
      out2fd = out1fd
      both = true
   }

   # Allocate the results vector
   mxnstr = 50000
   mxncol = 11
   szoutv = mxnstr * mxncol
   call malloc (outptr, szoutv, TY_DOUBLE)

   # quick and dirty peaker square aperture photometry
   call qdpsap(Memr[im1ptr], Memr[im2ptr], 
      im1nx, im1ny, im1dx, im1dy, mnpk1, mxpk1,
      Memd[outptr], szoutv, nstr)

   # print the photometry results
   call qdt1(Memd[outptr], szoutv, nstr, out1fd, out2fd)

   # print information about the program and parameters
   do j = 1, 2 {
      ok = true
      if (j == 1) { ofd = out1fd }
      if (j == 2) {
         if (both) {
            ok = false
         } else {
            ofd = out2fd
         }
      }
      if (ok) {
         call fprintf (ofd, "#\n")
         call fprintf (ofd, "# QDPHOT_V1.3 (21SEP99)\n")
         call fprintf (ofd, "#\n")
         call fprintf (ofd, "# Copyleft (L) 1999 Kenneth J. Mighell\n")

         call fprintf (ofd, "#\n")
         call fprintf (ofd, "# qdphot\n")
         call fprintf (ofd, "#    image1 = \"%s\"\n")
           call pargstr (Memc[im1fn])
         call fprintf (ofd, "#     gain1 = %g\n")
           call pargr (im1gain)
         call fprintf (ofd, "#   mnpeak1 = %g\n")
           call pargr (mnpk1)
         call fprintf (ofd, "#   mxpeak1 = %g\n")
           call pargr (mxpk1)
         call fprintf (ofd, "#    image2 = \"%s\"\n")
           call pargstr (Memc[im2fn])
         call fprintf (ofd, "#     gain2 = %g\n")
           call pargr (im1gain)
         call fprintf (ofd, "#      out1 = \"%s\"\n")
           call pargstr (Memc[out1])
         call fprintf (ofd, "#      out2 = \"%s\"\n")
           call pargstr (Memc[out2])

         call fprintf (ofd, "#\n")
         call fprintf (ofd, "# qdphot \"%s\" %g %g %g \"%s\" %g \"%s\" \"%s\"\n")
           call pargstr (Memc[im1fn])
           call pargr (im1gain)        	
           call pargr (mnpk1)        	
           call pargr (mxpk1)        	
           call pargstr (Memc[im2fn])
           call pargr (im2gain)        	
           call pargstr (Memc[out1])        	
           call pargstr (Memc[out2])        	

         call fprintf (ofd, "#\n")

      }
   }

bye_
   # free up allocated memory
   call mfree (im1ptr, TY_REAL)
   call imunmap (im1)
   call mfree (im2ptr, TY_REAL)
   call imunmap (im2)
   call mfree (outptr, TY_DOUBLE)
   call close (out1fd)
   call close (out2fd)
   call sfree (sp)
   call flush (STDOUT)

end
#