/*      FILE: fdels0d.c
 *   PURPOSE: Compute del(model)/del(I_0) (=PSF)
 *    AUTHOR: Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE: ANSI C
 *      DATE: 2001NOV09
 * COPYRIGHT: (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 *
 * COMMENT: *** SOON TO BE DIGITAL ***
 *
 *    MXNPARL=7  
 *  
 *  and
 *
 *    ad[0] is B0 (=B_0      constant value of background )
 *    ad[1] is BX (=B_X      X slope value of background )
 *    ad[2] is BY (=B_Y      Y slope value of background )
 *    ad[3] is S0 (=I_0      intensity of 2-d Gaussian )
 *    ad[4] is SX (=X_0      X coordinate of 2-d Gaussian )
 *    ad[5] is SY (=Y_0      Y coordiante of 2-d Gaussian )
 *    ad[6] is SS (=SIGMA_0  standard deviation of 2-d Gaussian )
 *
 *  where
 *
 *    Z = B_0 + + (X*B_X) + (Y*B_Y) + I_0*Gaussian( X, Y; X_0, Y_0, SIGMA_0)
 *
 *  is a Gaussian centered at (X_0,Y_0) on a tilted background
 *
 */

#include "mx.h"
#include "inc.h"


int fdels0d_i7(
  inc_PsfS *PsfS,
  struct mxip_image_s *image,
#ifdef nogo
  double **matrix,
#endif
  double **xmatrix,
  double **ymatrix,
  int nx,
  int ny,
  double ad[]
){
  char mxfunc[] = "fdels0d_i7";

  int
    xi,
    yi;
  const int res=5;
  double xd, yd, md, dzdi, dzdx, dzdy, dzds;
  int mxstatus=0;
  int status;

  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 td1;
  double td2;
  double td3;
  double td4;

  double pixeld;

  double dmxf;
  double dmyf;

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

  int debug = MX_TRUE;
#define DEBUG (debug==MX_TRUE)

  if (DEBUG) printf("\n*** fdels0d_i7 ***\n");

  sxp = ad[4];
  syp = ad[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 PsfS->PSF->matrixd */
  for (iy=0; iy<PsfS->nyi; ++iy) {
    for (ix=0; ix<PsfS->nxi; ++ix) {
      PsfS->PSF->matrixd[iy][ix] = 0.0;
    }
  }

  /* PsfS->PSF is PsfS->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;
    }
  }

  /* Initialize image->matrixd */
  for (yi=0; yi<ny; yi++) {
    for (xi=0; xi<nx; xi++) {
      image->matrixd[yi][xi] = 0.0;
    }
  }

  mxstatus++;
  status = inc_ImageS_Add_Star_i7 (
    image,
    1.000,
    ad[4],
    ad[5],
    PsfS->PSF, 
    pxj-1,
    pyj-1
  );
  if (status) goto mx_error;

#ifdef NOGO
if (DEBUG) printf("# yi,xi,image->matrixd[yi][xi],dzdi,td1-td2:\n");
  for (yi=0; yi<ny; yi++) {
    for (xi=0; xi<nx; xi++) {
      Gaussian2d_v12(
	ad[3], 
        ad[4], 
        ad[5], 
        ad[6], 
        res, 
        xmatrix[yi][xi],
        ymatrix[yi][xi],
	&md, 
        &dzdi, 
        &dzdx, 
        &dzdy, 
        &dzds
      );

      td1 = image->matrixd[yi][xi];
      td2 = dzdi;
      td3 = td1-td2;
      td4 = 44.4444444;
      if (td2!=0.) td4 = td1/td2;
      if (td2!=0.0) if (DEBUG) printf("%d %d %12.6e %12.6e %12.6e %12.6f\n",yi,xi,td1,td2,td3,td4);
      /*
       *      image->matrixd[yi][xi] = dzdi;
       */
    }
  }
#endif

ok:
  status = 0;
  goto bye;

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

bye:
  status=0;
error:
  return(status);
}
#undef DEBUG
/* end-of-file */