Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c	(revision 40660)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c	(revision 40661)
@@ -17,11 +17,21 @@
   Buffer *image;
 
-  if (argc != 3) {
-    gprint (GP_ERR, "USAGE: medimage add (name) (image)\n");
+  Buffer *var = NULL;
+  if ((N = get_argument (argc, argv, "-var"))) {
+    remove_argument (N, &argc, argv);
+    if ((var = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 3) {
+    gprint (GP_ERR, "USAGE: medimage add (name) (image) [-var var]\n");
     gprint (GP_ERR, "       add the given image to the set of images to be medianed\n");
+    gprint (GP_ERR, "       optionally supply variance image (for weighted calculations)\n");
     return FALSE;
   }
 
   if ((image = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+
+XXX: match dimensions wit var image
 
   MedImageType *median = FindMedImage (argv[1]);
@@ -48,8 +58,15 @@
   int Ninput = median->Ninput;
   median->Ninput ++;
-  REALLOCATE (median->buffers, float *, median->Ninput);
-
-  ALLOCATE (median->buffers[Ninput], float, median->Nx*median->Ny);
-  memcpy (median->buffers[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
+  REALLOCATE (median->sig, float *, median->Ninput);
+  REALLOCATE (median->var, float *, median->Ninput);
+
+  ALLOCATE (median->sig[Ninput], float, median->Nx*median->Ny);
+  memcpy (median->sig[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
+
+  median->var[Ninput] = NULL;
+  if (var) {
+    ALLOCATE (median->var[Ninput], float, median->Nx*median->Ny);
+    memcpy (median->var[Ninput], var->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
+  }
 
   return TRUE;
@@ -100,5 +117,5 @@
       int Npix = ix + Nx*iy;
       for (n = 0; n < Ninput; n++) {
-	float v = median->buffers[n][Npix];
+	float v = median->sig[n][Npix];
 	if (!isfinite(v)) continue;
 	value[N] = v;
@@ -106,4 +123,8 @@
       }
       if (N == 0) continue;
+
+    XXXX: add in calc options:
+      weighted mean
+      irls
 
       if (CALC_MEAN) {
Index: /branches/eam_branches/ohana.20190329/src/opihi/include/data.h
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/include/data.h	(revision 40660)
+++ /branches/eam_branches/ohana.20190329/src/opihi/include/data.h	(revision 40661)
@@ -48,5 +48,6 @@
   int Nx;
   int Ny;
-  float **buffers;
+  float **sig;
+  float **var;
 } MedImageType;
 
Index: /branches/eam_branches/ohana.20190329/src/opihi/lib.data/MedImageOps.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/lib.data/MedImageOps.c	(revision 40660)
+++ /branches/eam_branches/ohana.20190329/src/opihi/lib.data/MedImageOps.c	(revision 40661)
@@ -32,7 +32,7 @@
   free (medimage[0].name);
   for (i = 0; i < medimage[0].Ninput; i++) {
-    free (medimage[0].buffers[i]);
+    free (medimage[0].sig[i]);
   }
-  free (medimage[0].buffers);
+  free (medimage[0].sig);
   free (medimage);
 }
@@ -75,5 +75,5 @@
   medimage->Nx = Nx;
   medimage->Ny = Ny;
-  ALLOCATE (medimage->buffers, float *, 1);
+  ALLOCATE (medimage->sig, float *, 1);
 
   medimages[N] = medimage;
Index: /branches/eam_branches/ohana.20190329/src/opihi/lib.data/starfuncs.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/lib.data/starfuncs.c	(revision 40660)
+++ /branches/eam_branches/ohana.20190329/src/opihi/lib.data/starfuncs.c	(revision 40661)
@@ -75,4 +75,6 @@
   FWHMx = 2.355*sqrt (fabs(x2 / I - x*x));
   FWHMy = 2.355*sqrt (fabs(y2 / I - y*y));
+
+  fprintf (stderr, "Mxx, Myy: %f, %f\n", x2/I - x*x, y2/I - y*y);
   Sxy   = xy / I - x*y;
   mag = -2.5*log10(I);
Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkobj.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkobj.c	(revision 40660)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkobj.c	(revision 40661)
@@ -6,5 +6,6 @@
 
   // input parameters:
-  // * trace       : spline fit of slit central x pos vs y-coord
+  // * slit_trace  : spline fit of slit central x pos vs y-coord
+  // * psf_trace   : spline fit of psf central x pos vs y-coord
   // * profile     : slit window profile (vector)
   // * object      : vector of object flux vs y-coord 
@@ -18,21 +19,27 @@
 
   // if any of these are not defined, they will have assumed identity values
-  Spline *trace   = NULL;
-  Vector *profile = NULL;
-  Vector *object  = NULL;
-  Vector *sky     = NULL;
-  Vector *backgnd = NULL;
-
-  Vector *psf     = NULL;
-  Vector *lsf     = NULL;
+  Spline *slit_trace = NULL;
+  Spline *psf_trace  = NULL;
+  Vector *profile    = NULL;
+  Vector *object     = NULL;
+  Vector *sky        = NULL;
+  Vector *backgnd    = NULL;
+		     
+  Vector *psf        = NULL;
+  Vector *lsf        = NULL;
   // add tilt later
 
-  Buffer *output  = NULL;
+  Buffer *output     = NULL;
 
   float stilt = 0.0; // angle of the slit
 
-  if ((N = get_argument (argc, argv, "-trace"))) {
-    remove_argument (N, &argc, argv);
-    if ((trace = FindSpline (argv[N])) == NULL) return (FALSE);
+  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);
+  }
+  if ((N = get_argument (argc, argv, "-psf-trace"))) {
+    remove_argument (N, &argc, argv);
+    if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE);
     remove_argument (N, &argc, argv);
   }
@@ -82,5 +89,7 @@
 
   if (argc != 3) {
-    gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-Nwave N] [-object vector] [-sky vector] [-backgnd vector] [-trace spline] [-profile vector]\n");
+    gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-Nwave N] [-object vector] [-sky vector] [-backgnd vector] [-slit-trace spline] [-psf-trace spline] [-profile vector]\n");
+    gprint (GP_ERR, "  psf-trace : cross-dispersion shift *relative* to the slit trace\n");
+    gprint (GP_ERR, "  sky       : sky is relative to background outside slit\n");
     return FALSE;
   }
@@ -123,10 +132,10 @@
   // if we are appying a trace offset spline, we need to generate an output window which
   // is Nx + the full swing of the trace, then window back down
-  if (trace) {
+  if (slit_trace) {
     float dXmin = +1000;
     float dXmax = -1000;
     for (int iy = 0; iy < Ny; iy++) {
       // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
-      float dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
+      float dx = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
       dXmin = MIN (dx, dXmin);
       dXmax = MAX (dx, dXmax);
@@ -180,4 +189,9 @@
     int Npof = floor(Nx/2);
     if (psf) Npof -= Npsf;
+
+    if (psf_trace) {
+      float dx = -spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
+      Npof += dx;
+    }
 
     // flux = obj * PSF + sky
@@ -329,5 +343,5 @@
 
   // shift pixels in x based on the trace
-  if (trace) {
+  if (slit_trace) {
     ALLOCATE_PTR (outTraceBuffer, char, NxBase*Ny*sizeof(float));
     float *outTrace = (float *) outTraceBuffer;
@@ -337,5 +351,5 @@
 
       // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
-      float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
+      float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
       
       // extract the integer pixel offset and the fractional offset
