/* MX NATIVE IRAF TASK: mxqdphot
 *
 * PURPOSE:  MX version of QDPHOT 
 *
 * AUTHOR:  K. J. Mighell (mighell@noao.edu)
 *
 * LANGUAGE:  ANSI C
 *
 * DOCUMENTATION:  http://www.noao.edu/staff/mighell/mx
 *
 * DATE:  20000316
 *
 * MOD:  3
 * 
 * Copyleft (L) 2000 Kenneth John Mighell
 */

#include <mx.h>
#include "mxqdphot.h"

#define MXQDPHOT_DATE "16MAR2000"

IRAFTASK (mxqdphot) {
 int i;
 extern int iraf_mxqdphot0 (void);
 i = iraf_mxqdphot0 ();
}

/* 
 * ############################################################################
 */
int iraf_mxqdphot0(void) {

  char mxfunc[] = "iraf_mxqdphot0\7"; /* \7 rings the bell! */
  int status;

  struct mxiraf_image_s *im1=NULL; /* 1st image */
  char im1fn[MX_IRAF_SZ_FNAME+1];  /* .. filename */
  double im1gain;                  /* .. gain */
  double mnpk1;                    /* .. minimum peak */
  double mxpk1;                    /* .. maximum peak */

  struct mxiraf_image_s *im2=NULL; /* 2nd image */
  char im2fn[MX_IRAF_SZ_FNAME+1];  /* .. filename */
  double im2gain;                  /* .. gain */

  FILE *ofp=NULL;                  /* output pointer */
  FILE *out1fp=NULL;               /* 1st output pointer */
  char out1fn[MX_IRAF_SZ_FNAME+1]; /* .. filename */
  FILE *out2fp=NULL;               /* 2nd output pointer */
  char out2fn[MX_IRAF_SZ_FNAME+1]; /* .. filename */

  int ok;
  int both;
  int j;

  double *outd; /* output vector */
  int szoutd;   /* .. dimension */
  int mxnstr;   /* .. maximum number of stars */
  int mxncol;   /* .. maximum number of columns */

  int nstr;

/* CVOS forgot to define c_clgetr! */
#define c_clgetr c_clgetf

  status = mxp_init_f0();
  if (status) {status=1; goto mx_error;}

  /*
   * Get 1st image parameters
   */
  c_clgstr ("image1", im1fn, sizeof(im1fn));
  if (mxiraf_cvos_error_check_f0()) {status=2; goto mx_error;}
  im1gain = (double) c_clgetr ("gain1"); 
  if (mxiraf_cvos_error_check_f0()) {status=3; goto mx_error;}
  mnpk1 = (double) c_clgetr ("mnpeak1"); 
  if (mxiraf_cvos_error_check_f0()) {status=4; goto mx_error;}
  mxpk1 = (double) c_clgetr ("mxpeak1"); 
  if (mxiraf_cvos_error_check_f0()) {status=5; goto mx_error;}

  /*
   * Get 2nd image parameters
   */
  c_clgstr ("image2", im2fn, sizeof(im2fn));
  if (mxiraf_cvos_error_check_f0()) {status=6; goto mx_error;}
  im2gain = (double) c_clgetr ("gain2"); 
  if (mxiraf_cvos_error_check_f0()) {status=7; goto mx_error;}

  /*
   * Open 1st image and multiply it by gain1
   */
  status = mxiraf_image_s_open_f3 (&im1, im1fn, IRAF_READ_ONLY);
  if (status) {status = 8; goto mx_error;}
  status = mxip_image_s_const_mul_f2 (im1->data, im1gain);
  if (status) {status = 9; goto mx_error;}

  /*
   * Open 2nd image and multiply it by gain2
   */
  status = mxiraf_image_s_open_f3 (&im2, im2fn, IRAF_READ_ONLY);
  if (status) {status = 10; goto mx_error;}
  status = mxip_image_s_const_mul_f2 (im2->data, im2gain);
  if (status) {status = 11; goto mx_error;}

  ok = (im1->data->nxi==im2->data->nxi) 
    && (im1->data->nyi==im2->data->nyi);
  if (!ok) {
    status = 12;
    sprintf (MX.tmpmsg, "#   image section sizes must be identical!\n");
    goto mx_error2;
   }

  c_clgstr ("out1", out1fn, sizeof(out1fn));
  if (mxiraf_cvos_error_check_f0()) {status=12; goto mx_error;}
  if (out1fn[0]!='\0') {
    out1fp = fopen (out1fn, "w");
    if ((FILE *)NULL==out1fp) {
      status = 13; 
      sprintf (
        MX.tmpmsg, 
        "#   Could not open 1st output file:\n"
        "#   \"%s\"\n"
        "#   Does this file already exist?\n",
        out1fn
      );
      goto mx_error2;
    }
  } else {
    out1fp = stdout;
  }

  c_clgstr ("out2", out2fn, sizeof(out2fn));
  if (mxiraf_cvos_error_check_f0()) {status=14; goto mx_error;}
  both = MX_FALSE;
  if (out2fn[0]!='\0') {
    out2fp = fopen (out2fn, "w");
    if ((FILE *)NULL==out2fp) {
      status = 14; 
      sprintf (
        MX.tmpmsg, 
        "#   Could not open 1st output file:\n"
        "#   \"%s\"\n"
        "#   Does this file already exist?\n",
        out2fn
      );
      goto mx_error2;
    }
  } else {
    out2fp = out1fp;
    both = MX_TRUE;
  }

  /* Allocate and initialize the results vector */
  mxnstr = 50000; 
  mxncol = 11;
  szoutd = mxnstr * mxncol;
  outd = calloc ((mxnstr*mxncol), sizeof(double));
  if ((double *)NULL==outd) {
    status = 15; 
    sprintf (
      MX.tmpmsg, 
      "#   could not allocate output vector of %d bytes\n",
      szoutd
    );
    goto mx_error2;
  }
  for (j=0; j<szoutd; ++j) outd[j] = 0.0;

  /* quick and dirty peaker square aperture photometry */
  mxqdphot0_qdpsap_vf11 (
    im1->data->matrixd,
    im2->data->matrixd,
    im1->data->nxi,
    im1->data->nyi,
    im1->data->llx0i,
    im1->data->lly0i,
    mnpk1,
    mxpk1,
    outd,
    szoutd,
    &nstr
  );

  /* print the photometry results */
  mxqdphot0_qdt1_vf5 (outd, szoutd, nstr, out1fp, out2fp);

  for (j=1; j<=2; ++j) {
    ok = MX_TRUE;
    if (j == 1) ofp = out1fp;
    if (j == 2) {
      if (both) {
        ok = MX_FALSE;
      } else {
        ofp = out2fp;
      }
    }
    if (ok) {
      fprintf (ofp, "#\n");
      fprintf (ofp, 
	"# MXQDPHOT (" 
        MXQDPHOT_DATE 
        ")  "
        MX_COPYLEFT
	"\n"
      );
      fprintf (ofp, "#\n");
      fprintf (ofp, "# mxqdphot\n");
      fprintf (ofp, "#    image1 = \"%s\"\n", im1fn);
      fprintf (ofp, "#     gain1 = %g\n", im1gain);
      fprintf (ofp, "#   mnpeak1 = %g\n", mnpk1);
      fprintf (ofp, "#   mxpeak1 = %g\n", mxpk1);
      fprintf (ofp, "#    image2 = \"%s\"\n", im2fn);
      fprintf (ofp, "#     gain2 = %g\n", im1gain);
      fprintf (ofp, "#      out1 = \"%s\"\n", out1fn);
      fprintf (ofp, "#      out2 = \"%s\"\n", out2fn);
    }
  }

ok:
  status = 0;
  goto bye;

mx_error:
  mxp_errmsg_append_f3 (mxfunc, status, "");
  mxp_errmsg_dump_f0();

  goto bye;

mx_error2:
  mxp_errmsg_append_f3 (mxfunc, status, MX.tmpmsg);
  mxp_errmsg_dump_f0();
  goto bye;

bye:
  fflush (stdout);
  mxiraf_image_s_close_f1 (&im1);
  mxiraf_image_s_close_f1 (&im2);
  if ( ((FILE *)NULL != out1fp) && (stdout != out1fp) ) fclose(out1fp);
  if ( ((FILE *)NULL != out2fp) && (stdout != out2fp) ) fclose(out2fp);
  if ((double *)NULL!=outd) free(outd);
  return(status);

}