Index: trunk/Ohana/src/opihi/mana/Makefile
===================================================================
--- trunk/Ohana/src/opihi/mana/Makefile	(revision 40621)
+++ trunk/Ohana/src/opihi/mana/Makefile	(revision 40622)
@@ -21,7 +21,11 @@
 $(SRC)/init.$(ARCH).o \
 $(SRC)/mana.$(ARCH).o \
+$(SRC)/deimos_mkobj.$(ARCH).o \
+$(SRC)/deimos_mkslit.$(ARCH).o \
+$(SRC)/deimos_fitslit.$(ARCH).o \
 $(SRC)/findrowpeaks.$(ARCH).o 
 
 cmds = \
+$(SRC)/deimos.$(ARCH).o \
 $(SRC)/rawstars.$(ARCH).o \
 $(SRC)/fitcontour.$(ARCH).o \
Index: trunk/Ohana/src/opihi/mana/deimos.c
===================================================================
--- trunk/Ohana/src/opihi/mana/deimos.c	(revision 40621)
+++ trunk/Ohana/src/opihi/mana/deimos.c	(revision 40622)
@@ -1,31 +1,35 @@
 # include "data.h"
 
-int deimos_list (int argc, char **argv);
+int deimos_mkslit (int argc, char **argv);
+int deimos_fitslit (int argc, char **argv);
+int deimos_mkobj (int argc, char **argv);
 
-static Command book_commands[] = {
-  {1, "list",     book_list,     "list books"},
+static Command deimos_commands[] = {
+  {1, "mkobj", deimos_mkobj, "make a full object image"},
+  {1, "mkslit", deimos_mkslit, "make a slit image"},
+  {1, "fitslit", deimos_fitslit, "fit slit image to observed slit flux"},
 };
 
-int book_command (int argc, char **argv) {
+int deimos (int argc, char **argv) {
 
   int i, N, status;
 
   if (argc < 2) {
-    gprint (GP_ERR, "USAGE: book (command)\n");
-    gprint (GP_ERR, "    book list                                  : list books\n");
+    gprint (GP_ERR, "USAGE: deimos (command)\n");
+    gprint (GP_ERR, "    deimos mkslit : make slit image\n");
     return (FALSE);
   }
 
-  N = sizeof (book_commands) / sizeof (Command);
+  N = sizeof (deimos_commands) / sizeof (Command);
 
-  /* find the book sub-command which matches */
+  /* find the deimos sub-command which matches */
   for (i = 0; i < N; i++) {
-    if (!strcmp (book_commands[i].name, argv[1])) {
-      status = (*book_commands[i].func) (argc - 1, argv + 1);
+    if (!strcmp (deimos_commands[i].name, argv[1])) {
+      status = (*deimos_commands[i].func) (argc - 1, argv + 1);
       return (status);
     }
   }
 
-  gprint (GP_ERR, "unknown book command %s\n", argv[1]);
+  gprint (GP_ERR, "unknown deimos command %s\n", argv[1]);
   return (FALSE);
 }
Index: trunk/Ohana/src/opihi/mana/deimos_commands.c
===================================================================
--- trunk/Ohana/src/opihi/mana/deimos_commands.c	(revision 40621)
+++ trunk/Ohana/src/opihi/mana/deimos_commands.c	(revision 40622)
@@ -2,9 +2,4 @@
 
 int deimos_mkslit (int argc, char **argv) {
-  OHANA_UNUSED_PARAM(argv);
-  if (argc != 5) {
-    gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace) (flux) (buffer)\n");
-    return FALSE;
-  }
 
   // generate model observed flux given a uniformly illuminated slit
@@ -16,32 +11,111 @@
 
   Vector *profile = NULL;
-  Vector *flux = NULL;
-  Spline *trace = NULL;
-  Buffer *buff = NULL;
+  Vector *flux    = NULL;
+  Vector *sky     = NULL;
+  Spline *trace   = NULL;
+  Buffer *buff    = NULL;
+
+  // user-specified flux value (otherwise assumed to be 1)
+  if ((N = get_argument (argc, argv, "-flux"))) {
+    remove_argument (N, &argc, argv);
+    if ((flux = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+
+  // user-specified sky value (otherwise assumed to be 0)
+  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);
+  }
+
+  // either flux or fluxbins must be specified to define output size
+  int fluxbins = 0;
+  if ((N = get_argument (argc, argv, "-fluxbins"))) {
+    remove_argument (N, &argc, argv);
+    Nrows = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
+    return FALSE;
+  }
+
+  if (!flux && !fluxbins) {
+    gprint (GP_ERR, "either flux or fluxbins must be specified to define output size\n");
+    return FALSE;
+  }
 
   // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
   if ((profile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
   if ((trace   = FindSpline (argv[2])) == NULL) return (FALSE);
-  if ((flux    = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((buff    = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
-  
+  if ((buff    = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+
   // define the output window
   int Nx = 2*profile[0].Nelements;
-  int Ny = flux[0].Nelements;
+  int Ny = flux ? flux[0].Nelements : fluxbins;
   
   ResetBuffer (buff, Nx, Ny, -32, 0.0, 1.0);
-  
+
+  // the profile is registered to the midpoint of the vector : (int) N / 2
+  // the output window is registered to the midpoint of the window : (int) Nx / 2
+
+  int NxMidProfile = profile->Nelements / 2;
+  int NxMidWindow  = Nx / 2;
+
+  float *out = (float *) buff[0].matrix.buffer;
 
   // loop over the y positions
-  for (iy = 0; iy < flux->Nelements; iy++) {
+  for (int iy = 0; iy < Ny; iy++) {
 
-    // evaluate the trace spline at this y-coord to find the x-coord of the profile center
-    opihi_flt dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, 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);
     
-    // interpolate the profile with offset of dx
-    for (ix = 0; ix < profile->Nelements; ix++) {
-      // interpolate from profile->elements.Flt[ix], interpolate from profile->elements.Flt[ix+1] to ix+dx
+    // extract the integer pixel offset and the fractional offset
+    int dxi = floor(dx); // -1.7 -> -2, -0.5 -> -1, +0.5 -> 0, +1.7 -> 1
+    float dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
+
+    // starting point in output window is :
+    // XXX handle case if profile is wider than window
+    int Sx = (int)(NxMidWindow - NxMidProfile) + dxi;
+
+    // if fractional offset is small, do not interpolate
+    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
+
+    // loop over the relevant pixels in the window:
+    for (int ix = 0; ix < Nx; ix++) {
+
+      // equivalent coord in the profile:
+      int px = ix - Sx;
+      if (px < 0) continue;
+      if (px >= profile->Nelements) continue;
+
+      float vout = NAN;
+
+      if (doInterp) {
+	if ((px > 0) && (px < profile->Nelements - 1)) {
+	  // if ((px == 5) && (iy % 50 == 0)) {
+	  //   fprintf (stderr, "%d : %d %d : %d, %f\n", iy, Sx, px, dxi, dxf);
+	  // }
+	  // XXX something is not quite right in my math: I think (1-dx), dx are reversed
+	  // but if I use what I expect the interpolation is backwards...
+	  vout = profile->elements.Flt[px]*(1 - dxf) + profile->elements.Flt[px-1]*dxf;
+	}
+	if (px == 0) {
+	  vout = profile->elements.Flt[px]*dxf;
+	}
+	if (px == profile->Nelements - 1) {
+	  vout = profile->elements.Flt[profile->Nelements - 1]*(1.0 - dxf);
+	}
+      } else {
+	vout = profile->elements.Flt[px];
+      }
+
+      if (flux) vout *= flux->elements.Flt[iy];
+      if (sky) vout += sky->elements.Flt[iy];
+
+      out[ix + iy*Nx] = vout;
     }
-
   }
 
Index: trunk/Ohana/src/opihi/mana/deimos_fitslit.c
===================================================================
--- trunk/Ohana/src/opihi/mana/deimos_fitslit.c	(revision 40622)
+++ trunk/Ohana/src/opihi/mana/deimos_fitslit.c	(revision 40622)
@@ -0,0 +1,73 @@
+# include "data.h"
+
+int deimos_fitslit (int argc, char **argv) {
+
+  // fitslit (window) (slit) (flux) (sky)
+
+  // int N;
+
+  Vector *flux    = NULL;
+  Vector *sky     = NULL;
+
+  Buffer *wind    = NULL;
+  Buffer *slit    = NULL;
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: deimos fitslit (window) (slit) (flux) (sky)\n");
+    gprint (GP_ERR, "  inputs:  window (observed 2D flux), slit (model 2D flux)\n");
+    gprint (GP_ERR, "  outputs: flux (best-fit 1D flux), sky (best-fit 1D background)\n");
+    return FALSE;
+  }
+
+  // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
+  if ((wind = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((slit = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((flux = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((sky  = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  // XXX confirm slit and wind have same dimensions
+
+  // define the output window
+  int Nx = wind[0].matrix.Naxis[0];
+  int Ny = wind[0].matrix.Naxis[1];
+  
+  if (Nx != slit[0].matrix.Naxis[0]) return FALSE;
+  if (Ny != slit[0].matrix.Naxis[1]) return FALSE;
+
+  ResetVector (flux, OPIHI_FLT, Ny);
+  ResetVector (sky,  OPIHI_FLT, Ny);
+
+  float *Fwind = (float *) wind[0].matrix.buffer;
+  float *Fprof = (float *) slit[0].matrix.buffer;
+
+  // loop over the rows
+  for (int iy = 0; iy < Ny; iy++) {
+
+    // calculate elements of the chi-square for this row:
+
+    double R = 0.0, P = 0.0, P2 = 0.0, f = 0.0, fp = 0.0;
+    double wt = 1.0;
+
+    for (int ix = 0; ix < Nx; ix++) {
+      
+      float fxy = Fwind[ix + iy*Nx];
+      float pxy = Fprof[ix + iy*Nx];
+
+      R  += wt; 
+      P  += pxy*wt;
+      P2 += pxy*pxy*wt;
+      f  += fxy*wt;
+      fp += fxy*pxy*wt;
+    }
+
+    double det = R*P2 - P*P;
+    double S = (P2*f - P*fp) / det;
+    double F = (R*fp - P*f) / det;
+
+    flux->elements.Flt[iy] = F;
+    sky->elements.Flt[iy] = S;
+  }
+
+  return TRUE;
+}
+
Index: trunk/Ohana/src/opihi/mana/deimos_mkobj.c
===================================================================
--- trunk/Ohana/src/opihi/mana/deimos_mkobj.c	(revision 40622)
+++ trunk/Ohana/src/opihi/mana/deimos_mkobj.c	(revision 40622)
@@ -0,0 +1,176 @@
+# include "data.h"
+
+int deimos_mkobj (int argc, char **argv) {
+
+  // generate model observed flux for an object in a slit with background
+
+  // input parameters:
+  //   trace       : spline fit of slit central x pos vs y-coord
+  //   profile     : slit window profile (vector)
+  //   object      : vector of object flux vs y-coord
+  //   sky         : vector of local sky signal vs y-coord
+  //   background  : vector of extra-slit background flux vs y-coord
+  //   PSF         : point-spread function vector (flux normalized, x-dir)
+  //   LSF         : sline-spread function vector (flux normalized, y-dir)
+  //   stilt       : slit tilt response : 2D kernel? 
+
+  int N;
+
+  // 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;
+  // add tilt later
+
+  Buffer *output  = NULL;
+
+  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);
+  }
+  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);
+  }
+  if ((N = get_argument (argc, argv, "-object"))) {
+    remove_argument (N, &argc, argv);
+    if ((object = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+  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);
+  }
+  if ((N = get_argument (argc, argv, "-backgnd"))) {
+    remove_argument (N, &argc, argv);
+    if ((backgnd = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+  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);
+  }
+  if ((N = get_argument (argc, argv, "-lsf"))) {
+    remove_argument (N, &argc, argv);
+    if ((lsf = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+
+  // either flux or fluxbins must be specified to define output size
+  int Nwave = 0;
+  if ((N = get_argument (argc, argv, "-Nwave"))) {
+    remove_argument (N, &argc, argv);
+    Nwave = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-Nwave N] [-object vector] [-sky vector] [-backgnd vector] [-trace spline] [-profile vector]\n");
+    return FALSE;
+  }
+
+  int waveVector = (object != NULL) || (sky != NULL) || (backgnd != NULL);
+
+  if (!waveVector && !Nwave) {
+    gprint (GP_ERR, "either flux or fluxbins must be specified to define output size\n");
+    return FALSE;
+  }
+  if (waveVector && Nwave) {
+    gprint (GP_ERR, "-Nwave incompatible with -object, -sky, -backgnd\n");
+    return FALSE;
+  }
+
+  // all supplied vectors must be consistent 
+  int Ny = object ? object->Nelements : Nwave;
+  if (sky && Ny) {
+    if (Ny != sky->Nelements) {
+      gprint (GP_ERR, "inconsistent wavelength scales (sky)\n");
+      return FALSE;
+    }
+  } else {
+    Ny = sky ? sky->Nelements : Ny;
+  }
+  if (backgnd && Ny) {
+    if (Ny != backgnd->Nelements) {
+      gprint (GP_ERR, "inconsistent wavelength scales (backgnd)\n");
+      return FALSE;
+    }
+  } else {
+      Ny = backgnd ? backgnd->Nelements : Ny;
+  }
+
+  if ((output = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+  int Nx = atoi(argv[2]);
+  
+  ResetBuffer (output, Nx, Ny, -32, 0.0, 1.0);
+
+  opihi_flt *objV = (opihi_flt *) (object ? object->elements.Flt : NULL);
+  opihi_flt *skyV = (opihi_flt *) (sky    ? sky->elements.Flt : NULL);
+  opihi_flt *psfV = (opihi_flt *) (psf    ? psf->elements.Flt : NULL);
+
+  float *out = (float *) output[0].matrix.buffer;
+
+  // psf vector must have odd number of pixels:
+  int Npsf = 0;
+  int Npof = floor(Nx/2);
+  if (psf) {
+    if (psf->Nelements % 2 == 0) {
+      gprint (GP_ERR, "psf vector must have an odd number of pixels\n");
+      return FALSE;
+    }
+    Npsf = psf->Nelements / 2;
+    Npof -= Npsf;
+  }
+
+  // loop over the y positions
+  for (int iy = 0; iy < Ny; iy++) {
+
+    opihi_flt objVy = objV ? objV[iy] : 1.0; // if object is not supplied, flux defaults to 1.0
+    opihi_flt skyVy = skyV ? skyV[iy] : 0.0; // if sky is not supplied, flux defaults to 0.0
+      
+    // flux = obj * PSF + sky
+    for (int ix = 0; ix < Nx; ix++) {
+      
+      opihi_flt value = skyVy;
+
+      if (psfV) {
+	int n = ix - Npof; // n is the pixel in the PSF corresponding to the ix pixel in the output image
+	if (n < 0) continue;
+	if (n >= Npsf) continue;
+	
+	value += objVy * psfV[n];
+      } else {
+	if (ix == Nx / 2) value += objVy;
+      }	  
+      out[ix + iy*Nx] = value;
+    }
+  }
+  return TRUE;
+}
+
+  /* 
+    // convolve the object flux at this y-coord with PSF and add sky
+    for (int ix = 0; i < Nx; i++) {
+      opihi_flt g = s = 0;
+      for (int n = -Npsf; n <= Npsf; n++) {
+	if (ix + n < 0) continue;
+	if (ix + n >= Nx) continue;
+	opihi_flt value = objV ? objV[ix+n] : 1.0;
+	s += psfV[n]*value;
+	g += psfV[n];
+      }
+      out[ix + iy*Nx] = Npsf ? s / g : (objV ? objV[ix+n] : 1.0);
+      if (skyV) {
+	out[ix + iy*Nx] += skyV[iy]
+    }
+  */
+  
Index: trunk/Ohana/src/opihi/mana/deimos_mkslit.c
===================================================================
--- trunk/Ohana/src/opihi/mana/deimos_mkslit.c	(revision 40622)
+++ trunk/Ohana/src/opihi/mana/deimos_mkslit.c	(revision 40622)
@@ -0,0 +1,129 @@
+# include "data.h"
+
+int deimos_mkslit (int argc, char **argv) {
+
+  // generate model observed flux given a uniformly illuminated slit
+
+  // input parameters:
+  //   slit profile response : vector of fractional flux vs x-coord
+  //   trace                 : spline fit of slit central x pos vs y-coord
+  //   flux                  : vector of input signal vs y-coord
+
+  int N;
+
+  Vector *profile = NULL;
+  Vector *flux    = NULL;
+  Vector *sky     = NULL;
+  Spline *trace   = NULL;
+  Buffer *buff    = NULL;
+
+  // user-specified flux value (otherwise assumed to be 1)
+  if ((N = get_argument (argc, argv, "-flux"))) {
+    remove_argument (N, &argc, argv);
+    if ((flux = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+
+  // user-specified sky value (otherwise assumed to be 0)
+  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);
+  }
+
+  // either flux or fluxbins must be specified to define output size
+  int fluxbins = 0;
+  if ((N = get_argument (argc, argv, "-fluxbins"))) {
+    remove_argument (N, &argc, argv);
+    fluxbins = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
+    return FALSE;
+  }
+
+  if (!flux && !fluxbins) {
+    gprint (GP_ERR, "either flux or fluxbins must be specified to define output size\n");
+    return FALSE;
+  }
+
+  // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
+  if ((profile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((trace   = FindSpline (argv[2])) == NULL) return (FALSE);
+  if ((buff    = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+
+  // define the output window
+  int Nx = 2*profile[0].Nelements;
+  int Ny = flux ? flux[0].Nelements : fluxbins;
+  
+  ResetBuffer (buff, Nx, Ny, -32, 0.0, 1.0);
+
+  // the profile is registered to the midpoint of the vector : (int) N / 2
+  // the output window is registered to the midpoint of the window : (int) Nx / 2
+
+  int NxMidProfile = profile->Nelements / 2;
+  int NxMidWindow  = Nx / 2;
+
+  float *out = (float *) buff[0].matrix.buffer;
+
+  // loop over the y positions
+  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);
+    
+    // extract the integer pixel offset and the fractional offset
+    int dxi = floor(dx); // -1.7 -> -2, -0.5 -> -1, +0.5 -> 0, +1.7 -> 1
+    float dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
+
+    // starting point in output window is :
+    // XXX handle case if profile is wider than window
+    int Sx = (int)(NxMidWindow - NxMidProfile) + dxi;
+
+    // if fractional offset is small, do not interpolate
+    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
+
+    // loop over the relevant pixels in the window:
+    for (int ix = 0; ix < Nx; ix++) {
+
+      // set the sky in the entire window 
+      out[ix + iy*Nx] = sky ? sky->elements.Flt[iy] : 0.0;
+
+      // equivalent coord in the profile:
+      int px = ix - Sx;
+      if (px < 0) continue;
+      if (px >= profile->Nelements) continue;
+
+      // a default value:
+      float vout = NAN;
+
+      if (doInterp) {
+	if ((px > 0) && (px < profile->Nelements - 1)) {
+	  // if ((px == 5) && (iy % 50 == 0)) {
+	  //   fprintf (stderr, "%d : %d %d : %d, %f\n", iy, Sx, px, dxi, dxf);
+	  // }
+	  // XXX something is not quite right in my math: I think (1-dx), dx are reversed
+	  // but if I use what I expect the interpolation is backwards...
+	  vout = profile->elements.Flt[px]*(1 - dxf) + profile->elements.Flt[px-1]*dxf;
+	}
+	if (px == 0) {
+	  vout = profile->elements.Flt[px]*dxf;
+	}
+	if (px == profile->Nelements - 1) {
+	  vout = profile->elements.Flt[profile->Nelements - 1]*(1.0 - dxf);
+	}
+      } else {
+	vout = profile->elements.Flt[px];
+      }
+
+      if (flux) vout *= flux->elements.Flt[iy];
+
+      out[ix + iy*Nx] += vout;
+    }
+  }
+
+  return TRUE;
+}
+
Index: trunk/Ohana/src/opihi/mana/init.c
===================================================================
--- trunk/Ohana/src/opihi/mana/init.c	(revision 40621)
+++ trunk/Ohana/src/opihi/mana/init.c	(revision 40622)
@@ -1,4 +1,5 @@
 # include "opihi.h"
 
+int deimos	    PROTO((int, char **));
 int findpeaks	    PROTO((int, char **));
 int fitcontour	    PROTO((int, char **));
@@ -8,4 +9,5 @@
 
 static Command cmds[] = {  
+  {1, "deimos",      deimos,       "deimos multislit spectrograph tools"},
   {1, "findpeaks",   findpeaks,    "find image peaks"},
   {1, "fitcontour",  fitcontour,   "fit ellipse contour"},