Index: trunk/Ohana/src/opihi/mana/Makefile
===================================================================
--- trunk/Ohana/src/opihi/mana/Makefile	(revision 40654)
+++ trunk/Ohana/src/opihi/mana/Makefile	(revision 40655)
@@ -26,6 +26,8 @@
 $(SRC)/deimos_getobj.$(ARCH).o \
 $(SRC)/deimos_fitobj.$(ARCH).o \
+$(SRC)/deimos_fitalt.$(ARCH).o \
 $(SRC)/deimos_mkalt.$(ARCH).o \
 $(SRC)/deimos_mkmodel.$(ARCH).o \
+$(SRC)/deimos_arclines.$(ARCH).o \
 $(SRC)/findrowpeaks.$(ARCH).o 
 
Index: trunk/Ohana/src/opihi/mana/deimos.c
===================================================================
--- trunk/Ohana/src/opihi/mana/deimos.c	(revision 40654)
+++ trunk/Ohana/src/opihi/mana/deimos.c	(revision 40655)
@@ -7,7 +7,10 @@
 int deimos_getobj (int argc, char **argv);
 int deimos_fitobj (int argc, char **argv);
+int deimos_fitalt (int argc, char **argv);
+int deimos_arclines (int argc, char **argv);
 
 static Command deimos_commands[] = {
   {1, "fitobj", deimos_fitobj, "fit for object parameters"},
+  {1, "fitalt", deimos_fitalt, "fit for object parameters using LMM"},
   {1, "getobj", deimos_getobj, "determine crude object parameters"},
   {1, "mkobj", deimos_mkobj, "make a full object image"},
@@ -15,4 +18,5 @@
   {1, "mkslit", deimos_mkslit, "make a slit image"},
   {1, "fitslit", deimos_fitslit, "fit slit image to observed slit flux"},
+  {1, "arclines", deimos_arclines, "detect arclines using LSF and STILT"},
 };
 
Index: trunk/Ohana/src/opihi/mana/deimos_arclines.c
===================================================================
--- trunk/Ohana/src/opihi/mana/deimos_arclines.c	(revision 40655)
+++ trunk/Ohana/src/opihi/mana/deimos_arclines.c	(revision 40655)
@@ -0,0 +1,155 @@
+# include "data.h"
+
+// global parameters of the slit kernel:
+static int   Nkernel_pix = 0; // total number of pixels
+static int   Nkernel_y2  = 0;  // center of the interpolation 
+static float *kernel     = NULL; 
+
+void deimos_make_LSF_kernel (Vector *LSF, float stilt, int Nx);
+
+int deimos_arclines (int argc, char **argv) {
+
+  // arclines (buffer) (LSF) (STILT) (coord) (flux)
+
+  int N;
+  
+  Vector *LSF     = NULL;
+  Vector *flux    = NULL;
+  Vector *coord   = NULL;
+  Buffer *buff    = NULL;
+
+  Buffer *kern    = NULL;
+  if ((N = get_argument (argc, argv, "-save-kern"))) {
+    remove_argument (N, &argc, argv);
+    if ((kern = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 6) {
+    gprint (GP_ERR, "USAGE: deimos arclines (buffer) (LSF) (STILT) (coord) (flux)\n");
+    gprint (GP_ERR, "  inputs:  buffer (observed 2D flux), LSF (line spread function vector), STILT (degrees)\n");
+    gprint (GP_ERR, "  outputs: coord (dispersion coordinate vector), flux (intensity at coord)\n");
+    gprint (GP_ERR, "  options: [-save-kern (buffer)]\n");
+    return FALSE;
+  }
+
+  // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
+  if ((buff  = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((LSF   = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  float stilt = atof (argv[3]);
+  if ((coord = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((flux  = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  // define the output window
+  int Nx = buff[0].matrix.Naxis[0];
+  int Ny = buff[0].matrix.Naxis[1];
+  
+  ResetVector (coord, OPIHI_FLT, Ny);
+  ResetVector (flux,  OPIHI_FLT, Ny);
+
+  float *Fin = (float *) buff[0].matrix.buffer;
+
+  // generate a convolution kernel 
+  deimos_make_LSF_kernel (LSF, stilt, Nx);
+
+  // loop over the rows
+  for (int iy = 0; iy < Ny; iy++) {
+    float Fsum = 0;
+    // cross-correlation of buffer values and kernel values
+    for (int ix = 0; ix < Nx; ix++) {
+      for (int ky = -Nkernel_y2; ky <= Nkernel_y2; ky++) {
+	int ko = ky + Nkernel_y2;
+	if ((iy + ky) < 0) continue;
+	if ((iy + ky) >= Ny) continue;
+	Fsum += Fin[ix + (iy + ky)*Nx] * kernel[ix + ko*Nx];
+      }
+    }
+    coord->elements.Flt[iy] = iy;
+    flux->elements.Flt[iy] = Fsum;
+  }
+
+  if (kern) {
+    ResetBuffer (kern, Nx, 2*Nkernel_y2+1, -32, 0.0, 1.0);
+    float *kB = (float *) kern[0].matrix.buffer;
+    for (int iy = 0; iy < 2*Nkernel_y2 + 1; iy++) {
+      for (int ix = 0; ix < 2*Nx; ix++) {
+	kB[ix + iy*Nx] = kernel[ix + iy*Nx];
+      }
+    }
+  } 
+  free (kernel);
+
+  return TRUE;
+}
+
+// kernel for slit tilt; no LSF for now
+void deimos_make_LSF_kernel (Vector *LSF, float stilt, int Nx) {
+
+  // first, generate the interpolation kernels.  For a slit tilt of theta, there is a
+  // displacement in the y-direction of dy = dx sin(theta) where dx is the x-coord
+  // relative to the slit window center (Nx / 2).  we thus need an image of size Nx,
+  // Nx*sin(theta)
+
+  opihi_flt *LSFv = LSF->elements.Flt;
+
+  int Ns = LSF->Nelements;
+  int Ns2 = Ns/2;
+  
+  float sin_stilt = sin(stilt*RAD_DEG);
+  int Ny = Ns + ceil(Nx * fabs(sin_stilt)) + 1; // one extra pixel buffer
+  if (Ny < 3) Ny = 3; // minimum of 3 pixels (or kernel assumption below fails)
+
+  if (Ny % 2 == 0) Ny ++; // force Nyk to be odd
+
+  Nkernel_y2 = Ny/2; // Nkernel_y2 is the center pixel in the interpolation direction
+
+  ALLOCATE (kernel, float, Nx*Ny);
+  Nkernel_pix = Nx*Ny;
+  for (int ix = 0; ix < Nkernel_pix; ix++) { kernel[ix] = 0.0; }
+
+  // fprintf (stderr, "allocating %d pixels (%d,%d)\n", Nkernel_pix, Nx, Ny);
+
+  // XXX use Xref not Nx/2?
+  int Nx2 = floor(Nx/2);
+
+  for (int ix = 0; ix < Nx; ix++) {
+    // displacement in y-dir due to slit tilt:
+    float dy = (ix - Nx2) * sin_stilt;
+    int dyi = floor(dy);
+    float dyf = dy - dyi;
+
+    // offset of the LSF relative to the kernel window
+    int Sy = Nkernel_y2 - Ns2 + dyi;
+
+    // if fractional offset is small, do not interpolate
+    int doInterp = fabs(dyf) < 1e-5 ? FALSE : TRUE;
+
+    // loop over the LSF direction and insert into kernel
+    for (int iy = 0; iy < Ny; iy++) {
+
+      // equivalent coord in the LSF:
+      int py = iy - Sy;
+      if (py < 0) continue;
+      if (py >= Ns) continue;
+
+      // a default value:
+      float vout = NAN;
+
+      if (doInterp) {
+	if ((py > 0) && (py < Ns - 1)) {
+	  vout = LSFv[py]*(1 - dyf) + LSFv[py-1]*dyf;
+	}
+	if (py == 0) {
+	  vout = LSFv[py]*dyf;
+	}
+	if (py == Ns - 1) {
+	  vout = LSFv[Ns - 1]*(1.0 - dyf);
+	}
+      } else {
+	vout = LSFv[py];
+      }
+
+      kernel[ix + iy*Nx] = vout;
+    }
+  }
+}
Index: trunk/Ohana/src/opihi/mana/deimos_fitalt.c
===================================================================
--- trunk/Ohana/src/opihi/mana/deimos_fitalt.c	(revision 40655)
+++ trunk/Ohana/src/opihi/mana/deimos_fitalt.c	(revision 40655)
@@ -0,0 +1,447 @@
+# include "data.h"
+# include "deimos.h"
+
+# if (0) 
+int deimos_fitalt (int argc, char **argv) {
+  OHANA_UNUSED_PARAM(argc);
+  OHANA_UNUSED_PARAM(argv);
+  return FALSE;
+}
+
+# else
+
+/*
+  this is starting to work OK.  some improvements to make
+  * use Gaussdev to sample
+  * do not scale down range (or user-set scale-down)
+ */
+
+// internal functions to fitobj
+static float      deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts);
+double            deimos_LMM_update (float *data, float *model_ref, float *model_obj, float *model_sky, float *model_bck, float *dObj, float *dSky, float *dBck, int Nx, int Ny, int row, float lambda);
+static opihi_flt *deimos_make_test (opihi_flt *guess, float sigma, int row, int Nrow, float *dValue);
+opihi_flt        *deimos_subvector (opihi_flt *fullVec, int row, int Nrow);
+
+int deimos_fitalt (int argc, char **argv) {
+
+  // input parameters:
+  // * buffer      : 2D image of slit (full frame or cutout?)
+  // * trace       : spline fit of slit central x pos vs y-coord
+  // * profile     : slit window profile (vector)
+  // * PSF         : point-spread function vector (flux normalized, x-dir)
+  // * stilt       : slit tilt response : 2D kernel? 
+
+  // in-out parameters:
+  // * obj  : vector of object flux vs y-coord (starting guess and result)
+  // * sky  : vector of local sky signal vs y-coord 
+  // * bck  : vector of extra-slit background flux vs y-coord
+
+  int N;
+
+  Spline *trace   = NULL;
+  Vector *profile = NULL;
+
+  Vector *obj     = NULL;
+  Vector *sky     = NULL;
+  Vector *bck     = NULL;
+
+  Vector *psf     = NULL;
+
+  Buffer *buffer  = NULL;
+
+  float stilt = 0.0; // angle of the slit
+  ohana_gaussdev_init ();
+
+  int DO_PLUS = FALSE;
+  if ((N = get_argument (argc, argv, "-plus"))) {
+    remove_argument (N, &argc, argv);
+    DO_PLUS = TRUE;
+  }
+
+  int VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = TRUE;
+  }
+
+  // a crude noise model valid for both the guess vectors and the
+  // data buffer:
+  float gain = 1.0;
+  if ((N = get_argument (argc, argv, "-gain"))) {
+    remove_argument (N, &argc, argv);
+    gain = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  float noise = 0.0;
+  if ((N = get_argument (argc, argv, "-noise"))) {
+    remove_argument (N, &argc, argv);
+    noise = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  float nsigma = 2.0;
+  if ((N = get_argument (argc, argv, "-nsigma"))) {
+    remove_argument (N, &argc, argv);
+    nsigma = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  float lambda = 1.0;
+  if ((N = get_argument (argc, argv, "-lambda"))) {
+    remove_argument (N, &argc, argv);
+    lambda = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  int Niter = 25;
+  if ((N = get_argument (argc, argv, "-iter"))) {
+    remove_argument (N, &argc, argv);
+    Niter = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  int Nrow = 5;
+  if ((N = get_argument (argc, argv, "-row"))) {
+    remove_argument (N, &argc, argv);
+    Nrow = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  // Input parameters:
+  if ((N = get_argument (argc, argv, "-trace"))) {
+    remove_argument (N, &argc, argv);
+    if ((trace = FindSpline (argv[N])) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  } else { goto usage; }
+  if ((N = get_argument (argc, argv, "-profile"))) {
+    remove_argument (N, &argc, argv);
+    if ((profile = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  } else { goto usage; }
+  if ((N = get_argument (argc, argv, "-psf"))) {
+    remove_argument (N, &argc, argv);
+    if ((psf = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  } else { goto usage; }
+  if ((N = get_argument (argc, argv, "-stilt"))) {
+    remove_argument (N, &argc, argv);
+    stilt = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  } else { goto usage; }
+
+  // In-Out parameters:
+  if ((N = get_argument (argc, argv, "-object"))) {
+    remove_argument (N, &argc, argv);
+    if ((obj = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  } else { goto usage; }
+  if ((N = get_argument (argc, argv, "-sky"))) {
+    remove_argument (N, &argc, argv);
+    if ((sky = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  } else { goto usage; }
+  if ((N = get_argument (argc, argv, "-backgnd"))) {
+    remove_argument (N, &argc, argv);
+    if ((bck = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  } else { goto usage; }
+
+  if (argc != 3) goto usage;
+
+  // buffer is the full image, Xref is reference coordinate of the profile in the buffer
+  if ((buffer = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  
+  // profile and PSF are defined with central reference pixel mapped to Xref (+ trace)
+  if (profile->Nelements % 2 == 0) {
+    gprint (GP_ERR, "slit profile vector must have an odd number of pixels\n");
+    return FALSE;
+  }
+  if (psf->Nelements % 2 == 0) {
+    gprint (GP_ERR, "PSF vector must have an odd number of pixels\n");
+    return FALSE;
+  }
+
+  // observed buffer size
+  int Nx = buffer[0].matrix.Naxis[0];
+  int Ny = buffer[0].matrix.Naxis[1];
+  float *bufVal = (float *) buffer->matrix.buffer;
+
+  // obj, sky, bck must be consistent with data
+  if (Ny != obj->Nelements) {
+    gprint (GP_ERR, "inconsistent wavelength scales (object)\n");
+    return FALSE;
+  }
+  if (Ny != sky->Nelements) {
+    gprint (GP_ERR, "inconsistent wavelength scales (sky)\n");
+    return FALSE;
+  }
+  if (Ny != bck->Nelements) {
+    gprint (GP_ERR, "inconsistent wavelength scales (backgnd)\n");
+    return FALSE;
+  }
+
+  // for the cross-dispersion reference pixel, use the user value or set to Nx/2 if < 0
+  deimos_set_cross_ref (atoi(argv[2]), Nx); 
+  deimos_make_kernel (stilt, Nx);
+
+  Vector  *objMin = NULL;
+  Vector  *skyMin = NULL;
+  Vector  *bckMin = NULL;
+
+  if ((objMin = SelectVector ("objMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (objMin, OPIHI_FLT, Ny);
+  if ((skyMin = SelectVector ("skyMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (skyMin, OPIHI_FLT, Ny);
+  if ((bckMin = SelectVector ("bckMin", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bckMin, OPIHI_FLT, Ny);
+
+  // the functions below accept an opihi_flt array for object, sky, background
+  opihi_flt *objVal = objMin->elements.Flt;
+  opihi_flt *skyVal = skyMin->elements.Flt;
+  opihi_flt *bckVal = bckMin->elements.Flt;
+
+  // we make a copy of the input guess vectors and minimize them
+  for (int iy = 0; iy < Ny; iy++) {
+    objVal[iy] = obj->elements.Flt[iy];
+    skyVal[iy] = sky->elements.Flt[iy];
+    bckVal[iy] = bck->elements.Flt[iy];
+  }
+
+  float noiseVar = SQ(noise);
+  ALLOCATE_PTR (objNoise, float, Ny);
+  ALLOCATE_PTR (skyNoise, float, Ny);
+  ALLOCATE_PTR (bckNoise, float, Ny);
+
+  // generate noise vectors
+  for (int iy = 0; iy < Ny; iy++) {
+    objNoise[iy] = nsigma*sqrt(hypot(fabs(objVal[iy]) / gain, noiseVar));
+    skyNoise[iy] = nsigma*sqrt(hypot(fabs(skyVal[iy]) / gain, noiseVar));
+    bckNoise[iy] = nsigma*sqrt(hypot(fabs(bckVal[iy]) / gain, noiseVar));
+  }
+
+  // obj, sky, bck are initial guesses.  Also use these values to define a range for the initial guesses.
+
+  // XXX worry about last segment (< Nrow rows)
+  for (int iter = 0; iter < Niter; iter++) {
+
+    int dRow = Nrow/2;
+
+    // use scaled-noise to define test ranges
+    // we are going to try each row using a window Nrow wide to measure chisq:
+    for (int row = dRow; row < Ny - dRow - 1; row++) {
+      
+      // I am using a numerical equivalent of Levenberg-Marquardt, but only for a single
+      // p_m[i] = (obj,sky,bck) triplet [m is the index of obj,sky,bck; i is the row
+      // index] at a time.  I start with a current best guess set of (obj,sky,bck).  I
+      // need to calculate the model at this location, and then at offset locations p_m +
+      // dp_m.  Note that I do not need to calculate chisq for each of these to generate
+      // the elements of the LMM equations:
+
+      int row_0 = row - dRow; // start row of the subset test range
+
+      // Nrow = 5, dRow = 2, Ny = 100, last start is Ny - 2 - 1 = 97
+      // last row_0 is 95
+
+      // generate subset vectors for the range **centered** on row
+      opihi_flt *objCurr = deimos_subvector (objVal, row_0, Nrow);
+      opihi_flt *skyCurr = deimos_subvector (skyVal, row_0, Nrow);
+      opihi_flt *bckCurr = deimos_subvector (bckVal, row_0, Nrow);
+
+      // current vector value:
+      // for deimos_make_model, row is the starting point of the subimage
+      float *model_ref = deimos_make_model (objCurr, skyCurr, bckCurr, psf, profile, trace, Nx, Nrow, row_0);
+
+      // in the test functions below, dRow is the central pixel of the subset, e.g., obj[row] => objCurr[dRow]
+
+      float dObj = 0.0, dSky = 0.0, dBck = 0.0;
+
+      // delta obj (save dObj)
+      opihi_flt *objtest = deimos_make_test (objCurr, objNoise[row], dRow, Nrow, &dObj);
+      float *model_obj = deimos_make_model (objtest, skyCurr, bckCurr, psf, profile, trace, Nx, Nrow, row_0);
+
+      // delta sky (save dSky)
+      opihi_flt *skytest = deimos_make_test (skyCurr, skyNoise[row], dRow, Nrow, &dSky);
+      float *model_sky = deimos_make_model (objCurr, skytest, bckCurr, psf, profile, trace, Nx, Nrow, row_0);
+
+      // delta bck (save dBck)
+      opihi_flt *bcktest = deimos_make_test (bckCurr, bckNoise[row], dRow, Nrow, &dBck);
+      float *model_bck = deimos_make_model (objCurr, skyCurr, bcktest, psf, profile, trace, Nx, Nrow, row_0);
+
+      // dObj, dSky, dBck carry the delta used derivatives in and the delta applied to the values out
+      float chisq = deimos_LMM_update (bufVal, model_ref, model_obj, model_sky, model_bck, &dObj, &dSky, &dBck, Nx, Nrow, row_0, lambda);
+
+      if (1) {
+	if (DO_PLUS) {
+	  objVal[row] += dObj;
+	  skyVal[row] += dSky;
+	  bckVal[row] += dBck;
+	} else {
+	  objVal[row] -= dObj;
+	  skyVal[row] -= dSky;
+	  bckVal[row] -= dBck;
+	}
+      }
+
+      free (objCurr);
+      free (skyCurr);
+      free (bckCurr);
+
+      free (model_ref);
+
+      free (model_obj);
+      free (model_sky);
+      free (model_bck);
+
+      free (objtest);
+      free (skytest);
+      free (bcktest);
+
+      if (VERBOSE) fprintf (stderr, "chisq: %f (%f, %f, %f)\n", chisq, dObj, dSky, dBck);
+    }
+
+    float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, trace, Nx, Ny, 0);
+
+    int Npts = 0;
+    float chisq_full  = deimos_get_chisq (bufVal, model_full, Nx, Ny, 0, &Npts);
+
+    fprintf (stderr, "** Full chisq: %f (%d) : %f)\n", chisq_full, Npts, chisq_full / (1.0 * Npts));
+  }
+
+  return TRUE;
+
+ usage:
+  gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n");
+  return FALSE;
+}
+
+/****************** fitobj Support Functions *******************/
+
+// we generate a new vector with a single element (row) modified based on sigma
+opihi_flt *deimos_subvector (opihi_flt *fullVec, int row, int Nrow) {
+
+  ALLOCATE_PTR (value, opihi_flt, Nrow);
+
+  for (int i = 0; i < Nrow; i++) {
+    value[i] = fullVec[i + row];
+  }
+
+  return value;
+}
+
+
+// we generate a new vector with a single element (row) modified based on sigma
+static opihi_flt *deimos_make_test (opihi_flt *guess, float sigma, int row, int Nrow, float *dValue) {
+
+  // if fullInput is TRUE,  the guess vector runs from 0 to Ny (full wavelength range)
+  // if fullInput is FALSE, the guess vector runs from row to row + Nrow (subset range)
+  // sigma is valid for the element at row
+
+  ALLOCATE_PTR (value, opihi_flt, Nrow);
+
+  for (int i = 0; i < Nrow; i++) {
+    value[i] = guess[i];
+  }
+
+  value[row] = guess[row] + 0.1*sigma;
+
+  *dValue = 0.1*sigma;
+  return value;
+}
+
+
+double deimos_LMM_update (float *data, float *model_ref,
+			  float *model_obj, float *model_sky, float *model_bck,
+			  float *dObj, float *dSky, float *dBck, int Nx, int Ny, int row, float lambda) {
+  
+  // we are going to loop over the images calculating the following value:
+
+  // W = 1 / sigma^2, leave as 1.0 for now
+  // chisq = sum W * (data - model_ref)^2 
+  // dF_X = sum W * (data - model_ref) * (model_ref - model_X) / dX
+  // d2F_XY = sum W * (model_ref - model_X) * (model_ref - model_Y) / (dX * dY)
+
+  double chisq = 0.0;
+  double dF_obj = 0.0, dF_sky = 0.0, dF_bck = 0.0;
+  double d2F_obj_obj = 0.0, d2F_obj_sky = 0.0, d2F_obj_bck = 0.0, d2F_sky_sky = 0.0, d2F_sky_bck = 0.0, d2F_bck_bck = 0.0;
+
+  float W = 1.0;
+
+  for (int iy = 0; iy < Ny; iy++) {
+    for (int ix = 0; ix < Nx; ix++) {
+      int pix = ix + iy*Nx;
+      int pix_data  = ix + (iy + row)*Nx;
+
+      // if any of the pixel values are NAN, skip the point:
+      if (!isfinite(data[pix_data])) continue;
+      if (!isfinite(model_ref[pix])) continue;
+      if (!isfinite(model_obj[pix])) continue;
+      if (!isfinite(model_sky[pix])) continue;
+      if (!isfinite(model_bck[pix])) continue;
+
+      float df = W * (data[pix_data] - model_ref[pix]);
+
+      chisq += SQ(df);
+
+      float df_obj = (model_ref[pix] - model_obj[pix]) / *dObj;
+      float df_sky = (model_ref[pix] - model_sky[pix]) / *dSky;
+      float df_bck = (model_ref[pix] - model_bck[pix]) / *dBck;
+
+      dF_obj += df*df_obj;
+      dF_sky += df*df_sky;
+      dF_bck += df*df_bck;
+
+      d2F_obj_obj += W * SQ(df_obj);
+      d2F_sky_sky += W * SQ(df_sky);
+      d2F_bck_bck += W * SQ(df_bck);
+
+      d2F_obj_sky += W * df_obj * df_sky;
+      d2F_obj_bck += W * df_obj * df_bck;
+      d2F_sky_bck += W * df_sky * df_bck;
+    }
+  }
+
+  // these are the elements of Ax = B
+  ALLOCATE_PTR (A, double *, 3);
+  ALLOCATE_PTR (B, double *, 3);
+  for (int i = 0; i < 3; i++) {
+    ALLOCATE (A[i], double, 3);
+    ALLOCATE (B[i], double, 1);
+  }
+
+  A[0][0] = d2F_obj_obj * (1.0 + lambda);
+  A[1][1] = d2F_sky_sky * (1.0 + lambda);
+  A[2][2] = d2F_bck_bck * (1.0 + lambda);
+
+  A[0][1] = A[1][0] = d2F_obj_sky;
+  A[1][2] = A[2][1] = d2F_sky_bck;
+  A[0][2] = A[2][0] = d2F_obj_bck;
+
+  B[0][0] = dF_obj;
+  B[1][0] = dF_sky;
+  B[2][0] = dF_bck;
+
+  dgaussjordan (A, B, 3, 1);
+
+  *dObj = B[0][0];
+  *dSky = B[1][0];
+  *dBck = B[2][0];
+
+  return chisq;
+}
+
+static float deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts) {
+
+  int npts = 0;
+  float chisq = 0;
+
+  for (int iy = 0; iy < Ny; iy++) {
+
+    for (int ix = 0; ix < Nx; ix++) {
+
+      int pix_buffer = ix + (iy + row)*Nx;
+      int pix_model  = ix +  iy*Nx;
+
+      if (!isfinite(buffer[pix_buffer])) continue;
+      if (!isfinite(model[pix_model])) continue;
+      chisq += SQ(buffer[pix_buffer] - model[pix_model]);
+      npts ++;
+    }
+  }
+  *Npts = npts;
+  return chisq;
+}
+# endif
Index: trunk/Ohana/src/opihi/mana/deimos_mkmodel.c
===================================================================
--- trunk/Ohana/src/opihi/mana/deimos_mkmodel.c	(revision 40654)
+++ trunk/Ohana/src/opihi/mana/deimos_mkmodel.c	(revision 40655)
@@ -37,5 +37,5 @@
   ALLOCATE (kernel, float, Nx*Ny);
   Nkernel_pix = Nx*Ny;
-  fprintf (stderr, "allocating %d pixels (%d,%d)\n", Nkernel_pix, Nx, Ny);
+  // fprintf (stderr, "allocating %d pixels (%d,%d)\n", Nkernel_pix, Nx, Ny);
 
   // XXX use Xref not Nx/2?
@@ -53,13 +53,13 @@
     int pix;
     pix = ix + (dyi + Nkernel_y2 + 0)*Nx;
-    fprintf (stderr, "%d, %d, %d -- ", ix, (dyi + Nkernel_y2 + 0), pix);
+    // fprintf (stderr, "%d, %d, %d -- ", ix, (dyi + Nkernel_y2 + 0), pix);
     myAssert (pix < Nkernel_pix, "oops (make_kernel 1)");
     kernel[pix] = 1.0 - dyf;
 
     pix = ix + (dyi + Nkernel_y2 + 1)*Nx;
-    fprintf (stderr, "%d, %d, %d -- ", ix, (dyi + Nkernel_y2 + 1), pix);
+    // fprintf (stderr, "%d, %d, %d -- ", ix, (dyi + Nkernel_y2 + 1), pix);
     myAssert (pix < Nkernel_pix, "oops (make_kernel 2)");
     kernel[pix] = dyf;
-    fprintf (stderr, ": %f %f\n", kernel[ix + (dyi + Nkernel_y2 + 0)*Nx], kernel[ix + (dyi + Nkernel_y2 + 1)*Nx]);
+    // fprintf (stderr, ": %f %f\n", kernel[ix + (dyi + Nkernel_y2 + 0)*Nx], kernel[ix + (dyi + Nkernel_y2 + 1)*Nx]);
   }
 }
@@ -85,5 +85,5 @@
 
 // generate a subimage of the full model for the range row - row + Ny
-// obj,sky,bck are vectors of the full length of the model (ie, not sub-vectors)
+// obj,sky,bck are subset vectors corresponding to the range row to row+Ny
 float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *trace, int Nx, int Ny, int row) {
 
Index: trunk/Ohana/src/opihi/mana/deimos_mkobj.c
===================================================================
--- trunk/Ohana/src/opihi/mana/deimos_mkobj.c	(revision 40654)
+++ trunk/Ohana/src/opihi/mana/deimos_mkobj.c	(revision 40655)
@@ -260,5 +260,5 @@
       myAssert (Nkpix_i < Nkpix, "oops");
       kernel[Nkpix_i] = dyf;
-      fprintf (stderr, "%d, %d: %f %f\n", ix, (dyi + Nyk2), kernel[ix + (dyi + Nyk2 + 0)*Nx], kernel[ix + (dyi + Nyk2 + 1)*Nx]);
+      // fprintf (stderr, "%d, %d: %f %f\n", ix, (dyi + Nyk2), kernel[ix + (dyi + Nyk2 + 0)*Nx], kernel[ix + (dyi + Nyk2 + 1)*Nx]);
     }
   }