Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile	(revision 40752)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile	(revision 40753)
@@ -31,4 +31,5 @@
 $(SRC)/deimos_mkmodel.$(ARCH).o \
 $(SRC)/deimos_arclines.$(ARCH).o \
+$(SRC)/deimos_fitprofile.$(ARCH).o \
 $(SRC)/findrowpeaks.$(ARCH).o 
 
Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos.c	(revision 40752)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos.c	(revision 40753)
@@ -10,4 +10,5 @@
 int deimos_fitalt (int argc, char **argv);
 int deimos_arclines (int argc, char **argv);
+int deimos_fitprofile (int argc, char **argv);
 
 static Command deimos_commands[] = {
@@ -21,4 +22,5 @@
   {1, "fitslit", deimos_fitslit, "fit slit image to observed slit flux"},
   {1, "arclines", deimos_arclines, "detect arclines using LSF and STILT"},
+  {1, "fitprofile", deimos_fitprofile, "detect arclines using LSF and STILT"},
 };
 
Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitprofile.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitprofile.c	(revision 40753)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitprofile.c	(revision 40753)
@@ -0,0 +1,233 @@
+# include "data.h"
+
+typedef opihi_flt FitFunc (opihi_flt, opihi_flt *, int, opihi_flt *);
+typedef enum {MODE_UPPER, MODE_LOWER} SigmoidMode;
+
+void sigmoid_fit (Vector *xprofile, Vector *fprofile, Vector *oprofile, double *par, int Npar, SigmoidMode mode, double limit);
+
+opihi_flt fsigmoidDeimos (opihi_flt, opihi_flt *, int, opihi_flt *);
+opihi_flt rsigmoidDeimos (opihi_flt, opihi_flt *, int, opihi_flt *);
+
+// use these globals to carry into sigmoid_fit
+int VERBOSE = FALSE;
+int QUIET   = FALSE;
+
+int deimos_fitprofile (int argc, char **argv) {
+
+  // fit simple slit window profile model to vector
+  // I am fitting a sigmoid at both edges
+
+  // input parameters:
+  // * xprofile    : extracted slit window coordinates (vector)
+  // * xprofile    : extracted slit window profile (vector)
+  // * Xs, Xe      : scale coordinates of window boundaries
+
+  // output parameters:
+  // * sigmoid_left_center, etc : coordinates of window boundaries
+
+  int N;
+  Vector *xprofile = NULL;
+  Vector *fprofile = NULL;
+
+  QUIET = FALSE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    QUIET = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-quiet"))) {
+    QUIET = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    VERBOSE = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 5) goto usage;
+  if ((xprofile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((fprofile = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if (xprofile->Nelements != fprofile->Nelements) goto usage;
+
+  // supplied guess for boundaries
+  double Xs = atof(argv[3]);
+  double Xe = atof(argv[4]);
+
+  // output vectors of fits
+  Vector *Lprofile = SelectVector ("sigmoid_left", ANYVECTOR, TRUE);
+  if (Lprofile == NULL) return (FALSE);
+  MatchVector (Lprofile, xprofile, OPIHI_FLT);
+  Vector *Rprofile = SelectVector ("sigmoid_right", ANYVECTOR, TRUE);
+  if (Rprofile == NULL) return (FALSE);
+  MatchVector (Rprofile, xprofile, OPIHI_FLT);
+
+  // first, use supplied Xs & Xe guess to make a guess for BckL, Sky, BckR
+  
+  float Sky_V = 0.0, BckL_V = 0.0, BckR_V = 0.0;
+  int Sky_N = 0, BckL_N = 0, BckR_N = 0;
+  for (int i = 0; i < xprofile->Nelements; i++) {
+    if (xprofile->elements.Flt[i] < Xs) {
+      BckL_V += fprofile->elements.Flt[i];
+      BckL_N ++;
+      continue;
+    }
+    if (xprofile->elements.Flt[i] > Xe) {
+      BckR_V += fprofile->elements.Flt[i];
+      BckR_N ++;
+      continue;
+    }
+    Sky_V += fprofile->elements.Flt[i];
+    Sky_N ++;
+  }
+  float Sky = Sky_V / Sky_N;
+  float BckL = BckL_V / BckL_N;
+  float BckR = BckR_V / BckR_N;
+
+  // temporary storage vectors 
+  opihi_flt par[4];
+
+  float Xmid = 0.5*(Xs + Xe);
+
+  int Npar = 4;
+  par[0] = Xs;
+  par[1] = 1;
+  par[2] = Sky;
+  par[3] = BckL;
+
+  sigmoid_fit (xprofile, fprofile, Lprofile, par, Npar, MODE_LOWER, Xmid);
+
+  par[0] = Xe;
+  par[1] = 1;
+  par[2] = Sky;
+  par[3] = BckR;
+
+  sigmoid_fit (xprofile, fprofile, Rprofile, par, Npar, MODE_UPPER, Xmid);
+ 
+  return TRUE;
+  
+usage:
+  gprint (GP_ERR, "USAGE: deimos fitprofile xprof fprof Xs Xe\n");
+  return FALSE;
+}
+
+/* pars: x_o, sigma, top, bottom */
+opihi_flt rsigmoidDeimos (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
+  OHANA_UNUSED_PARAM(Npar);
+
+  double z = (x - par[0])/par[1];
+  double r = exp (z);
+  double s = 1 + r;
+  double q = 1 / s;
+  double f = par[2]*q + par[3];
+
+  // df/dp0 = df/dq * dq/dr * dr/dz * dz/dp0
+  // df/dq  = par[2]
+  // dq/ds  = -s^-2 = -q/s
+  // ds/dr  = 1
+  // dr/dz  = r
+  // dz/dp0 = -1 / par[1]
+  // dz/dp1 = -z / par[1]
+
+  // df/dp0 = par[2]*(-q/s)*r*(-1/par[1]) = par[2]*q*r/(s*par[1])
+  // df/dp1 = par[2]*(-q/s)*r*(-z/par[1]) = par[2]*q*r*z/(s*par[1])
+
+  if (dpar) {
+    dpar[0] = par[2]*q*r/(s*par[1]); 
+    dpar[1] = dpar[0]*z;
+    dpar[2] = q;
+    dpar[3] = 1;
+  }
+  
+  return (f);
+}
+
+/* pars: x_o, sigma, top, bottom */
+opihi_flt fsigmoidDeimos (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
+  OHANA_UNUSED_PARAM(Npar);
+
+  double z = (x - par[0])/par[1];
+  double r = exp (-z);
+  double s = 1 + r;
+  double q = 1 / s;
+  double f = par[2]*q + par[3];
+
+  // df/dp0 = df/dq * dq/dr * dr/dz * dz/dp0
+  // df/dq  = par[2]
+  // dq/ds  = -s^-2 = -q/s
+  // ds/dr  = 1
+  // dr/dz  = -r
+  // dz/dp0 = -1 / par[1]
+  // dz/dp1 = -z / par[1]
+
+  // df/dp0 = -par[2]*(-q/s)*r*(-1/par[1]) = par[2]*q*r/(s*par[1])
+  // df/dp1 = -par[2]*(-q/s)*r*(-z/par[1]) = par[2]*q*r*z/(s*par[1])
+
+  if (dpar) {
+    dpar[0] = -par[2]*q*r/(s*par[1]); 
+    dpar[1] = dpar[0]*z;
+    dpar[2] = q;
+    dpar[3] = 1;
+  }
+  
+  return (f);
+}
+
+void sigmoid_fit (Vector *xprofile, Vector *fprofile, Vector *oprofile, double *par, int Npar, SigmoidMode mode, double limit) {
+
+  int Nelements = xprofile->Nelements;
+
+  ALLOCATE_PTR (xtemp, opihi_flt, Nelements);
+  ALLOCATE_PTR (ytemp, opihi_flt, Nelements);
+  ALLOCATE_PTR (wtemp, opihi_flt, Nelements);
+
+  // extract the points on the right side and fit
+  // XXX error model is pure poisson (allow external errors)
+
+  int Npts = 0;
+  for (int i = 0; i < Nelements; i++) {
+    int keep = (mode == MODE_UPPER) ^ (xprofile->elements.Flt[i] < limit);
+    if (keep) {
+      xtemp[Npts] = xprofile->elements.Flt[i];
+      ytemp[Npts] = fprofile->elements.Flt[i];
+      wtemp[Npts] = 1.0 / (fabs(fprofile->elements.Flt[i] + 0.01));
+      Npts ++;
+    }
+  }
+
+  FitFunc *FUNC = (mode == MODE_LOWER) ? fsigmoidDeimos : rsigmoidDeimos;
+
+  double ochisq = mrqinit (xtemp, ytemp, wtemp, Npts, par, Npar, FUNC, VERBOSE);
+  double dchisq = ochisq + 2*Npts;
+
+  int Niter;
+  for (Niter = 0; (Niter < 20) && ((dchisq > 0.1*(Npts - Npar)) || (dchisq <= 0.0)); Niter++) {
+    double chisq = mrqmin (xtemp, ytemp, wtemp, Npts, par, Npar, FUNC, VERBOSE);
+    dchisq = ochisq - chisq;
+    ochisq = chisq;
+    if (VERBOSE) gprint (GP_ERR, "dchisq: %f, Ndof: %d\n", dchisq, Npts - Npar);
+  }  
+  if (VERBOSE) gprint (GP_ERR, "%d iterations\n", Niter); 
+
+  if (!QUIET) gprint (GP_LOG, "sigmoid @ %f (%f) : top %f bottom %f\n", par[0], par[1], par[2], par[3]);
+  for (int i = 0; i < Nelements; i++) {
+    oprofile[0].elements.Flt[i] = FUNC (xprofile[0].elements.Flt[i], par, Npar, NULL);
+  }
+
+  if (mode == MODE_LOWER) {
+    set_variable ("sigmoid_left_center", par[0]);
+    set_variable ("sigmoid_left_sigma",  par[1]);
+    set_variable ("sigmoid_left_top",    par[2]);
+    set_variable ("sigmoid_left_bottom", par[3]);
+  } else {
+    set_variable ("sigmoid_right_center", par[0]);
+    set_variable ("sigmoid_right_sigma",  par[1]);
+    set_variable ("sigmoid_right_top",    par[2]);
+    set_variable ("sigmoid_right_bottom", par[3]);
+  }
+  mrqfree (Npar);
+
+  FREE (xtemp);
+  FREE (ytemp);
+  FREE (wtemp);
+
+}