/*      FILE: fdelsyd.c
 *   PURPOSE: Compute del(model)/del(S_Y) [partial derivative w.r.t. Y]
 *    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 fdelsyd_i5(
  double **delsy,
  double **psf,
  double s0, 
  int nx,
  int ny
)
{
  const int
    xdim = nx,
    ydim = ny,
    xmin = 0,
    xmax = (xdim-1),
    ymin = 0,
    ymax = (ydim-1),
    yminp1 = (ymin+1),
    yminp2 = (ymin+2),
    ymaxm1 = (ymax-1),
    ymaxm2 = (ymax-2),
    yw = ymax - ymin + 1;
  int
    status,
    ix,
    iy;

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

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

  /*
   * Work on the outer 2 pixels on each side using a 
   * 2 or 3 point formula (where appropriate).
   */
  for (ix=xmin; ix<=xmax; ix++) {
    delsy[ymin][ix] = psf[ymin][ix] - psf[yminp1][ix];
    delsy[ymax][ix] = 0.0;
    delsy[yminp1][ix] = (psf[ymin][ix] - psf[yminp2][ix])/2.0;
    delsy[ymaxm1][ix] = (psf[ymaxm2][ix] - psf[ymax][ix])/2.0;
  }

  /*
   * Work on the inner pixels using a 5-point formula
   */
  for (iy=yminp2; iy<=ymaxm2; iy++) {
    for (ix=xmin; ix<=xmax; ix++) {
      delsy[iy][ix] = (
        -psf[iy-2][ix] + 8.0*psf[iy-1][ix] - 8.0*psf[iy+1][ix] + psf[iy+2][ix]
      )/12.0;
    }
  }

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

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