Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile	(revision 40837)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile	(revision 40838)
@@ -28,9 +28,10 @@
 $(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)/deimos_fitprofile.$(ARCH).o \
 $(SRC)/findrowpeaks.$(ARCH).o 
+
+# $(SRC)/deimos_mkalt.$(ARCH).o
 
 cmds = \
Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_arclines.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_arclines.c	(revision 40837)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_arclines.c	(revision 40838)
@@ -53,8 +53,14 @@
   deimos_make_LSF_kernel (LSF, stilt, Nx);
 
+  // FILE *Fout = NULL;
+
   // loop over the rows
   for (int iy = 0; iy < Ny; iy++) {
     float Fsum = 0;
     float Ksum = 0;
+
+    // if (iy == 5802) Fout = fopen ("test.5802.dat", "w");
+    // if (iy == 5815) Fout = fopen ("test.5815.dat", "w");
+
     // cross-correlation of buffer values and kernel values
     for (int ix = 0; ix < Nx; ix++) {
@@ -67,6 +73,12 @@
 	Fsum += Fin[ix + (iy + ky)*Nx] * kernel[ix + ko*Nx];
 	Ksum += kernel[ix + ko*Nx];
+	// if ((iy == 5815) || (iy == 5802)) {
+	//   fprintf (Fout, "%d %d %d : %f %f\n", iy, ix, ky, Fin[ix + (iy + ky)*Nx], kernel[ix + ko*Nx]);
+	// }
       }
     }
+
+    // if ((iy == 5815) || (iy == 5802)) fclose (Fout);
+
     coord->elements.Flt[iy] = iy;
     flux->elements.Flt[iy] = Fsum;
Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitalt.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitalt.c	(revision 40837)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitalt.c	(revision 40838)
@@ -1,13 +1,4 @@
 # 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
 
 /*
@@ -18,6 +9,6 @@
 
 // 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 float      deimos_get_chisq (float *buffer, float *vardata, float *model, int Nx, int Ny, int row, int *Npts);
+double            deimos_LMM_update (float *data, float *vardata, 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);
@@ -26,10 +17,11 @@
 
   // input parameters:
-  // * buffer      : 2D image of slit (full frame or cutout?)
-  // * slit_trace  : spline fit of slit central x pos vs y-coord
-  // * psf_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? 
+  // * buffer         : 2D image of slit (full frame or cutout?)
+  // * slit_trace_red : spline fit of slit central x pos vs y-coord
+  // * slit_trace_blu : spline fit of slit central x pos vs y-coord
+  // * psf_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:
@@ -40,17 +32,19 @@
   int N;
 
-  Spline *slit_trace = NULL;
+  Spline *slit_trace_red = NULL;
+  Spline *slit_trace_blu = NULL;
   Spline *psf_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
+
+  Vector *profile    = NULL;
+  Vector *psf        = NULL;
+
+  Vector *obj        = NULL;
+  Vector *sky        = NULL;
+  Vector *bck        = NULL;
+
+  Buffer *buffer     = NULL;
+
+  float stilt        = 0.0; // angle of the slit
+
   ohana_gaussdev_init ();
 
@@ -65,4 +59,20 @@
     remove_argument (N, &argc, argv);
     VERBOSE = TRUE;
+  }
+
+  int SKIPFIT = FALSE;
+  if ((N = get_argument (argc, argv, "-skipfit"))) {
+    remove_argument (N, &argc, argv);
+    SKIPFIT = TRUE;
+  }
+
+  // for a red vs blu spline, we need to specify the split point
+  // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
+  // this more generic (need to define the ranges somewhere)
+  int redlimit = 4096;
+  if ((N = get_argument (argc, argv, "-redlimit"))) {
+    remove_argument (N, &argc, argv);
+    redlimit = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
   }
 
@@ -108,9 +118,15 @@
 
   // Input parameters:
-  if ((N = get_argument (argc, argv, "-slit-trace"))) {
-    remove_argument (N, &argc, argv);
-    if ((slit_trace = FindSpline (argv[N])) == NULL) return (FALSE);
-    remove_argument (N, &argc, argv);
-  } else { goto usage; }
+  if ((N = get_argument (argc, argv, "-slit-trace-red"))) {
+    remove_argument (N, &argc, argv);
+    if ((slit_trace_red = FindSpline (argv[N])) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  } else { goto usage; }
+  if ((N = get_argument (argc, argv, "-slit-trace-blu"))) {
+    remove_argument (N, &argc, argv);
+    if ((slit_trace_blu = FindSpline (argv[N])) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  } else { goto usage; }
+
   if ((N = get_argument (argc, argv, "-psf-trace"))) {
     remove_argument (N, &argc, argv);
@@ -150,4 +166,20 @@
     remove_argument (N, &argc, argv);
   } else { goto usage; }
+
+  // save the model
+  Buffer *model = NULL;
+  if ((N = get_argument (argc, argv, "-model"))) {
+    remove_argument (N, &argc, argv);
+    if ((model = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+
+  // save the model
+  Buffer *varim = NULL;
+  if ((N = get_argument (argc, argv, "-variance"))) {
+    remove_argument (N, &argc, argv);
+    if ((varim = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
 
   if (argc != 3) goto usage;
@@ -171,4 +203,10 @@
   float *bufVal = (float *) buffer->matrix.buffer;
 
+  if (varim) {
+    if (Nx != varim[0].matrix.Naxis[0]) { gprint (GP_ERR, "inconsistent sizes for image and varim\n"); return FALSE; }
+    if (Ny != varim[0].matrix.Naxis[1]) { gprint (GP_ERR, "inconsistent sizes for image and varim\n"); return FALSE; }
+  }
+  float *varVal = varim ? (float *) varim->matrix.buffer : NULL;
+
   // obj, sky, bck must be consistent with data
   if (Ny != obj->Nelements) {
@@ -183,4 +221,8 @@
     gprint (GP_ERR, "inconsistent wavelength scales (backgnd)\n");
     return FALSE;
+  }
+
+  if (model) {
+    if (!CreateBuffer (model, Nx, Ny, -32, 1.0, 0.0)) { gprint (GP_ERR, "error generating output model buffer\n"); return FALSE; }
   }
 
@@ -230,5 +272,5 @@
     // 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++) {
+    for (int row = dRow; !SKIPFIT && (row < Ny - dRow - 1); row++) {
       
       // I am using a numerical equivalent of Levenberg-Marquardt, but only for a single
@@ -251,5 +293,5 @@
       // 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, slit_trace, psf_trace, Nx, Nrow, row_0);
+      float *model_ref = deimos_make_model (objCurr, skyCurr, bckCurr, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Nrow, row_0);
 
       // in the test functions below, dRow is the central pixel of the subset, e.g., obj[row] => objCurr[dRow]
@@ -259,16 +301,20 @@
       // 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, slit_trace, psf_trace, Nx, Nrow, row_0);
+      float *model_obj = deimos_make_model (objtest, skyCurr, bckCurr, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_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, slit_trace, psf_trace, Nx, Nrow, row_0);
+      float *model_sky = deimos_make_model (objCurr, skytest, bckCurr, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_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, slit_trace, psf_trace, Nx, Nrow, row_0);
+      float *model_bck = deimos_make_model (objCurr, skyCurr, bcktest, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_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);
+      float chisq = deimos_LMM_update (bufVal, varVal, model_ref, model_obj, model_sky, model_bck, &dObj, &dSky, &dBck, Nx, Nrow, row_0, lambda);
+
+      if (FALSE && (row > 10) && (row < 20)) {
+	fprintf (stderr, "%d %d : %f %f %f\n", iter, row, dObj, dSky, dBck);
+      }
 
       if (1) {
@@ -301,8 +347,15 @@
     }
 
-    float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, slit_trace, psf_trace, Nx, Ny, 0);
+    float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, slit_trace_red, slit_trace_blu, redlimit, psf_trace, Nx, Ny, 0);
 
     int Npts = 0;
-    float chisq_full  = deimos_get_chisq (bufVal, model_full, Nx, Ny, 0, &Npts);
+    float chisq_full  = deimos_get_chisq (bufVal, varVal, model_full, Nx, Ny, 0, &Npts);
+
+    if (model) {
+      free (model[0].matrix.buffer);
+      model[0].matrix.buffer = (char *) model_full;
+    } else {
+      free (model_full);
+    }
 
     fprintf (stderr, "** Full chisq: %f (%d) : %f)\n", chisq_full, Npts, chisq_full / (1.0 * Npts));
@@ -312,5 +365,5 @@
 
  usage:
-  gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -slit-trace spline -psf-trace spline -profile vector -psf vector -stilt (angle)\n");
+  gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -slit-trace-red (spline) -slit-trace-blu (spline) -psf-trace (spline) -profile vector -psf vector -stilt (angle)\n");
   return FALSE;
 }
@@ -351,5 +404,5 @@
 
 
-double deimos_LMM_update (float *data, float *model_ref,
+double deimos_LMM_update (float *data, float *vardata, 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) {
@@ -365,6 +418,4 @@
   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++) {
@@ -380,4 +431,8 @@
       if (!isfinite(model_bck[pix])) continue;
 
+      if (vardata && !isfinite(vardata[pix_data])) continue;
+      float Var = vardata ? vardata[pix_data] : 1.0;
+      float W = (fabs(Var) < 0.01) ? 100.0 : 1.0 / Var;
+
       float df = W * (data[pix_data] - model_ref[pix]);
 
@@ -431,5 +486,5 @@
 }
 
-static float deimos_get_chisq (float *buffer, float *model, int Nx, int Ny, int row, int *Npts) {
+static float deimos_get_chisq (float *buffer, float *vardata, float *model, int Nx, int Ny, int row, int *Npts) {
 
   int npts = 0;
@@ -445,5 +500,12 @@
       if (!isfinite(buffer[pix_buffer])) continue;
       if (!isfinite(model[pix_model])) continue;
-      chisq += SQ(buffer[pix_buffer] - model[pix_model]);
+
+      if (vardata && !isfinite(vardata[pix_buffer])) continue;
+      float Var = vardata ? vardata[pix_buffer] : 1.0;
+      float W = (fabs(Var) < 0.01) ? 100.0 : 1.0 / Var;
+
+      float dchisq = SQ(buffer[pix_buffer] - model[pix_model]) * W;
+      if (!isfinite(dchisq)) continue;
+      chisq += dchisq;
       npts ++;
     }
@@ -452,3 +514,2 @@
   return chisq;
 }
-# endif
Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkalt.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkalt.c	(revision 40837)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkalt.c	(revision 40838)
@@ -1,4 +1,6 @@
 # include "data.h"
 # include "deimos.h"
+
+XXX: deprecated
 
 // Use this to test the deimos_make_object function used by deimos_fitobj
Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkmodel.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkmodel.c	(revision 40837)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkmodel.c	(revision 40838)
@@ -86,5 +86,6 @@
 // generate a subimage of the full model for the range row - row + Ny
 // 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 *slit_trace, Spline *psf_trace, int Nx, int Ny, int row) {
+// the slit_trace_* and psf_trace refer to the full image starting at row
+float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, Spline *psf_trace, int Nx, int Ny, int row) {
 
   float *raw = deimos_make_straight_image (obj, sky, psf, psf_trace, Nx, Ny, row);
@@ -92,5 +93,5 @@
   deimos_apply_profile (profile, new, Nx, Ny);
   deimos_add_background (bck, new, Nx, Ny, row);
-  float *out = deimos_apply_trace (slit_trace, new, Nx, Ny, row);
+  float *out = deimos_apply_trace (slit_trace_red, slit_trace_blu, redlimit, new, Nx, Ny, row);
 
   free (raw);
@@ -113,5 +114,5 @@
 
   // Xref is a global supplied by the user (defaults to Nx/2 if < 0)
-  int Xoff = (int)(Xref - Npsf);
+  int XoffRef = (int)(Xref - Npsf);
 
   // Nx : width of output window
@@ -123,4 +124,5 @@
     
     // we are generating the image for just the row range row to row+Ny
+    // note that obj & sky are subset vectors for just this range of pixels
     opihi_flt objVy = obj[iy]; // if we shift to using the row offset: obj[iy + row]
     opihi_flt skyVy = sky[iy]; // if we shift to using the row offset: sky[iy + row]
@@ -129,11 +131,14 @@
     int dxi = 0;
     float dxf = 0.0, dxr = 1.0;
+    int Xoff = XoffRef;
     
     if (psf_trace) {
-      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
+      // the psf_trace is referenced against the full image, so we need to add row to iy:
+      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy + row);
       dxi = floor(dx);
       dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
       dxr = 1 - dxf;
-      Xref += dxi;
+      Xoff += dxi;
+      // fprintf (stderr, "%d: %f : %d\n", iy + row, dx, Xoff);
     }
 
@@ -153,11 +158,17 @@
 	if (doInterp) {
 	  if ((n > 0) && (n < NpsfFull)) {
-	    value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
+	    if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
+	      value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
+	    }
 	  }
 	  if (n == 0) {
-	    value += objVy*psfV[n]*dxr;
+	    if (isfinite(psfV[n])) {
+	      value += objVy*psfV[n]*dxr;
+	    }
 	  }
 	} else {
-	  value += objVy * psfV[n];
+	  if (isfinite(psfV[n])) {
+	    value += objVy * psfV[n];
+	  }
 	}
       }
@@ -243,5 +254,5 @@
 
 // shift pixels in x based on the trace
-float *deimos_apply_trace (Spline *slit_trace, float *input, int Nx, int Ny, int row) {
+float *deimos_apply_trace (Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, float *input, int Nx, int Ny, int row) {
 
   ALLOCATE_PTR (output, float, Nx*Ny);
@@ -251,5 +262,7 @@
     // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
     // NOTE: spline is evaluated at full dispersion coordinate
-    float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy + row);
+    int iy_full = iy + row;
+    Spline *slit_trace = (iy_full < redlimit) ? slit_trace_red : slit_trace_blu;
+    float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy_full);
       
     // extract the integer pixel offset and the fractional offset
Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkobj.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkobj.c	(revision 40837)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkobj.c	(revision 40838)
@@ -223,4 +223,5 @@
       dxr = 1 - dxf;
       Npof += dxi;
+      // fprintf (stderr, "%d: %f : %d\n", iy, dx, Npof);
       if (VERBOSE && (iy % 200 == 0)) { gprint (GP_ERR, "%d : %f : %d : %f %f\n", iy, dx, dxi, dxf, dxr); }
     }
@@ -241,11 +242,17 @@
 	  if (doInterp) {
 	    if ((n > 0) && (n < NpsfFull)) {
-	      value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
+	      if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
+		value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
+	      }
 	    }
 	    if (n == 0) {
-	      value += objVy*psfV[n]*dxr;
+	      if (isfinite(psfV[n])) {
+		value += objVy*psfV[n]*dxr;
+	      }
 	    }
 	  } else {
-	    value += objVy*psfV[n];
+	    if (isfinite(psfV[n])) {
+	      value += objVy*psfV[n];
+	    }
 	  }
 	}