/*      FILE: mxtools/src/psfstar/psfstar.c
 *   PURPOSE: Shift a MATPHOT PSF file
 *    AUTHOR: Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE: ANSI C
 *      DATE: 2001NOV06
 * COPYRIGHT: (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 */

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

IRAFTASK (psfstar) {
 int i;
 extern int iraf_psfstar (void);
 i = iraf_psfstar ();
}

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

  char mxfunc[] = "iraf_psfstar";
  int mxstatus=0;
  int status;
  int debug = MX_TRUE;
#define DEBUG (debug==MX_TRUE)

  double pixeld;

  double sxp;
  double syp;
  int sxj;
  int syj;
  double sxf;
  double syf;

  double pxp;
  double pyp;
  int pxj;
  int pyj;
  double pxf;
  double pyf;

  double dxf;
  double dyf;

  double dmxf;
  double dmyf;

  int ix;
  int iy;
  int jx;
  int jy;
  int kx;
  int ky;

  inc_ParS ParS;
  inc_ImageS ImageS;
  inc_PsfS PsfS;

  mxstatus++;
  status = mxp_init_f0();
  if (status) goto mx_error;

  mxstatus++;
  status = inc_ParS_Get_i1 (&ParS);
  if (status) goto mx_error;

  mxstatus++;
  status = inc_ImageS_Open_i2 (&ParS, &ImageS);
  if (status) goto mx_error;

  mxstatus++;
  status = inc_PsfS_Open_i2 (&ParS, &PsfS);
  if (status) goto mx_error;

  /*
  ** =========================================================================
  */

  sxp = ParS.sx - ImageS.llx0i - 0.5;
  syp = ParS.sy - ImageS.lly0i - 0.5;
  if (DEBUG) printf("# sxp is %g\n",sxp);
  if (DEBUG) printf("# syp is %g\n",syp);
  sxj = sxp;
  syj = syp;
  if (DEBUG) printf("# sxj is %d\n",sxj);
  if (DEBUG) printf("# syj is %d\n",syj);
  sxf = sxp - sxj;
  syf = syp - syj;
  if (DEBUG) printf("# sxf is %g\n",sxf);
  if (DEBUG) printf("# syf is %g\n",syf);

  pxp = PsfS.PSF_CPXP;
  pyp = PsfS.PSF_CPYP;
  if (DEBUG) printf("# pxp is %g\n",pxp);
  if (DEBUG) printf("# pyp is %g\n",pyp);
  pxj = pxp;
  pyj = pyp;
  if (DEBUG) printf("# pxj is %d\n",pxj);
  if (DEBUG) printf("# pyj is %d\n",pyj);
  pxf = pxp - pxj;
  pyf = pyp - pyj;
  if (DEBUG) printf("# pxf is %g\n",pxf);
  if (DEBUG) printf("# pyf is %g\n",pyf);

  dxf = sxf - pxf;
  dyf = syf - pyf;
  if (DEBUG) printf("# dxf is %g\n",dxf);
  if (DEBUG) printf("# dyf is %g\n",dyf);

  /* Shift PSF accounting for OverSamplingFactor */
  mxstatus++;
  dmxf = dxf*PsfS.PSF_MX;
  dmyf = dyf*PsfS.PSF_MY;
  status = ip_Shift2d_i4(
    PsfS.psf0,
    dmxf,
    dmyf,
    PsfS.psf
  );
  if (status) goto mx_error;

  pxp += dxf;
  pyp += dyf;
  if (DEBUG) printf("# pxp is now %g\n",pxp);
  if (DEBUG) printf("# pyp is now %g\n",pyp);
  pxj = pxp;
  pyj = pyp;
  if (DEBUG) printf("# pxj is now %d\n",pxj);
  if (DEBUG) printf("# pyj is now %d\n",pyj);
  pxf = pxp - pxj;
  pyf = pyp - pyj;
  if (DEBUG) printf("# pxf is now %g\n",pxf);
  if (DEBUG) printf("# pyf is now %g\n",pyf);

  /* Initialize the PSF array */
  for (iy=0; iy<PsfS.nyi; ++iy) {
    for (ix=0; ix<PsfS.nxi; ++ix) {
      PsfS.PSF->matrixd[iy][ix] = 0.0;
    }
  }

  /* PSF is psf rebinned by OverSamplingFactor */
  for (iy=0; iy<PsfS.PSF_NY; ++iy) {
    for (ix=0; ix<PsfS.PSF_NX; ++ix) {
      pixeld = 0.0;
      for (jy=0; jy<PsfS.PSF_MY; ++jy) {
	ky = (iy*PsfS.PSF_MY)+jy;
	for (jx=0; jx<PsfS.PSF_MX; ++jx) {
	  kx = (ix*PsfS.PSF_MX)+jx;
	  pixeld += PsfS.psf->matrixd[ky][kx];
	}
      }
      PsfS.PSF->matrixd[iy][ix] = pixeld;
    }
  }

  mxstatus++;
  status = inc_ImageS_Add_Star_i7 (
    ImageS.mxIrafImage->data, 
    ParS.se, 
    sxp,
    syp,
    PsfS.PSF,
    pxj-1,
    pyj-1
  );
  if (status) goto mx_error;

  /*
  ** =========================================================================
  */

  mxstatus++;
  status = inc_ParS_Put_i1 (&ParS);
  if (status) goto mx_error;

  inc_ParS_Print_v2 (&ParS, stdout);

ok:
  status = 0;
  goto bye;

mx_error:
mx_error1:
  mxp_tmpmsg_init_f0();
mx_error2:
  mxp_errmsg_append_f3 (mxfunc, mxstatus, MX.tmpmsg);
  mxp_tmpmsg_init_f0();
  mxp_errmsg_dump_f0();
  goto bye;

bye:
  inc_ImageS_Close_i1( &ImageS );
  inc_PsfS_Close_i1( &PsfS );
  fflush (stdout);
  return(status);

}
#undef DEBUG
/* end-of-file */