/*      FILE: fdelsxd.c
 *   PURPOSE: Compute del(model)/del(S_X) [partial derivative w.r.t. X]
 *    AUTHOR: Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE: ANSI C
 *      DATE: 2001SEP10
 * COPYRIGHT: (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 *
 * COMMENT: *** DIGITAL ***
 *
 */

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

#define DEBUGT (MX_TRUE)
#define DEBUGF (MX_FALSE)
#define DEBUG (DEBUGT)

int fdelsxd_i5(
  double **delsx,
  double **psf,
  double s0,
  int nx,
  int ny
)
{
  const int
    xdim = nx,
    ydim = ny,
    xmin = 0,
    xmax = (xdim-1),
    ymin = 0,
    ymax = (ydim-1),
    xminp1 = (xmin+1),
    xminp2 = (xmin+2),
    xmaxm1 = (xmax-1),
    xmaxm2 = (xmax-2),
    xw = xmax - xmin + 1;
  int
    status,
    ix,
    iy;

  if (xw<5) {
    status = 1;
    goto bye;
  }

  /* initialize */
  for (iy=ymin; iy<=ymax; iy++) {
    for (ix=xmin; ix<=xmax; ix++) {
      delsx[iy][ix] = 0.0;
    }
  }

  /*
   * Work on the outer 2 pixels on each side using a 
   * 2 or 3 point formula (where appropriate).
   */
  for (iy=ymin; iy<=ymax; iy++) {
    delsx[iy][xmin] = psf[iy][xmin] - psf[iy][xminp1];
    delsx[iy][xmax] = 0.0;
    delsx[iy][xminp1] = (psf[iy][xmin] - psf[iy][xminp2])/2.0;
    delsx[iy][xmaxm1] = (psf[iy][xmaxm2] - psf[iy][xmax])/2.0;
  }

  /*
   * Work on the inner pixels using a 5-point formula
   */
  for (iy=ymin; iy<=ymax; iy++) {
    for (ix=xminp2; ix<=xmaxm2; ix++) {
      delsx[iy][ix] = (
        -psf[iy][ix-2] + 8.0*psf[iy][ix-1] - 8.0*psf[iy][ix+1] + psf[iy][ix+2]
      )/12.0;
    }
  }

  /* and multiply by the intenstiy (s0) */
  for (iy=ymin; iy<=ymax; iy++) {
    for (ix=xmin; ix<=xmax; ix++) {
      delsx[iy][ix] *= s0;
    }
  }

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