Index: /trunk/Ohana/src/opihi/cmd.astro/imfit.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 2937)
+++ /trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 2938)
@@ -1,19 +1,49 @@
 # include "astro.h"
+# define FFACTOR 200
+# define FSCALE 1.2
+
+float mrq2dchi (float *x, float *t, float *y, float *dy, int Npts, 
+		float *par, int Npar, 
+		float (funcs)(float, float, float *, int, float *));
 
 float fgaussTD (float, float, float *, int, float *);
 float pgaussTD (float, float, float *, int, float *);
+float sgaussTD (float, float, float *, int, float *);
 float vgaussTD (float, float, float *, int, float *);
-float qgaussTD (float, float, float *, int, float *);
+
+float PgaussTD (float, float, float *, int, float *);
+float SgaussTD (float, float, float *, int, float *);
+
+float tgaussTD (float, float, float *, int, float *);
 float rgaussTD (float, float, float *, int, float *);
-float sgaussTD (float, float, float *, int, float *);
 
 static float Npow;
+static float Npin;
+static int MissingVariable = FALSE;
+static float *fpar;
+static int Nfpar;
+
+int mrq2dlimits (float *pmin, float *pmax, int Npar);
+
+float get_variable_default (char *name, float dvalue) {
+
+  char *value;
+  float fvalue;
+
+  value = get_variable (name);
+  if (value == NULL) {
+    MissingVariable = TRUE;
+    return (dvalue);
+  }
+  fvalue = atof (value);
+  return (fvalue);
+}
 
 int imfit (int argc, char **argv) {
 
-  int i, j, N, Npts, Npar, Save, VERBOSE;
+  int i, j, N, Npts, Npar, Save, VERBOSE, ShapeVariation;
   int sx, sy, nx, ny, Nx, Ny;
   float chisq, ochisq, dchisq, Gain, RDnoise;
-  float *x, *y, *z, *dz, *par, *V;
+  float *x, *y, *z, *dz, *par, *V, *parmin, *parmax;
   char line[64];
   Buffer *buf;
@@ -26,4 +56,10 @@
   }
 
+  ShapeVariation = FALSE;
+  if ((N = get_argument (argc, argv, "-shapes"))) {
+    remove_argument (N, &argc, argv);
+    ShapeVariation = TRUE;
+  }
+
   /* Gain in e/DN */
   Gain = 1.0;
@@ -49,4 +85,5 @@
 
   /* set fitting function */
+  Nfpar = 0;
   Npar = 7;
   fitfunc = fgaussTD;
@@ -60,19 +97,29 @@
       fitfunc = pgaussTD;
     }
-    if (!strcmp(argv[N], "qgauss")) {
-      fitfunc = qgaussTD;
-      Npar = 9;
-    }
     if (!strcmp(argv[N], "vgauss")) {
       fitfunc = vgaussTD;
       Npar = 9;
     }
-    if (!strcmp(argv[N], "rgauss")) {
-      fitfunc = rgaussTD;
-      Npar = 11;
-    }
     if (!strcmp(argv[N], "sgauss")) {
       fitfunc = sgaussTD;
       Npar = 10;
+    }
+    if (!strcmp(argv[N], "rgauss")) {
+      fitfunc = rgaussTD;
+      Npar = 10;
+    }
+    if (!strcmp(argv[N], "tgauss")) {
+      fitfunc = tgaussTD;
+      Npar = 10;
+    }
+    if (!strcmp(argv[N], "Pgauss")) {
+      fitfunc = PgaussTD;
+      Npar = 4;
+      Nfpar = 3;
+    }
+    if (!strcmp(argv[N], "Sgauss")) {
+      fitfunc = SgaussTD;
+      Npar = 4;
+      Nfpar = 6;
     }
     if (fitfunc == NULL) {
@@ -102,6 +149,10 @@
   if (sy + ny >= Ny) goto range;
 
-  ALLOCATE (par, float, 2*Npar);
+  ALLOCATE (par, float, Npar);
   bzero (par, Npar*sizeof(float));
+  if (Nfpar > 0) {
+    ALLOCATE (fpar, float, Nfpar);
+    bzero (fpar, Nfpar*sizeof(float));
+  }    
 
   /* convert array z[x,y] to x[i], y[i], z[i] */
@@ -127,27 +178,42 @@
 
   /* make initial guesses on parameters */
-  par[0] = atof (get_variable ("Xg"));
-  par[1] = atof (get_variable ("Yg"));
-  par[2] = 2.35 * sqrt(2.0) / atof (get_variable ("SXg"));
-  par[3] = 2.35 * sqrt(2.0) / atof (get_variable ("SYg"));
-  par[4] = 0.0;
-  par[5] = atof (get_variable ("Zpk"));
-  par[6] = atof (get_variable ("Sg"));
-  if (Npar == 9) {
+  if (Npar >= 7) {
+    par[0] = get_variable_default ("Xg", 0);
+    par[1] = get_variable_default ("Yg", 0);
+    par[2] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 2.0);
+    par[3] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 2.0);
+    par[4] = 0.0;
+    par[5] = get_variable_default ("Zpk", 10000);
+    par[6] = get_variable_default ("Sg", 0.0);
+  }
+  if (Npar == 4) {
+    par[0] = get_variable_default ("Xg", 0);
+    par[1] = get_variable_default ("Yg", 0);
+    par[2] = get_variable_default ("Zpk", 10000);
+    par[3] = get_variable_default ("Sg", 0.0);
+  }
+  if (Npar == 9) { /** vgauss **/
     par[7] = par[8] = 1;
   }
-  if (Npar == 11) {
-    par[7] = 0.2*par[5];
-    par[8] = 0.25*par[2];
-    par[9] = 0.25*par[3];
-    par[10] = 0.0;
-  }
-  if (Npar == 10) {
-    Npow = atof (get_variable ("Npow"));
-    par[7] = 2.35 * sqrt(2.0) / atof (get_variable ("SXf"));
-    par[8] = 2.35 * sqrt(2.0) / atof (get_variable ("SYf"));
+  if (Npar == 10) { /** sgauss, tgauss, rgauss **/
+    Npow = get_variable_default ("Npow", 2.25);
+    Npin = get_variable_default ("Npin", 1.00);
+    par[7] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0);
+    par[8] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0);
     par[9] = 0.0;
   }
-
+  if (fitfunc == PgaussTD) {
+    fpar[0] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 15.0);
+    fpar[1] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 15.0);
+    fpar[2] = get_variable_default ("SXYg", 15.0);
+  }    
+  if (fitfunc == SgaussTD) {
+    fpar[0] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 15.0);
+    fpar[1] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 15.0);
+    fpar[2] = get_variable_default ("SXYg", 15.0);
+    fpar[3] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0);
+    fpar[4] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0);
+    fpar[5] = get_variable_default ("SXYf", 15.0);
+  }    
 
   /* run fit routine */
@@ -159,15 +225,76 @@
     dchisq = ochisq - chisq;
     ochisq = chisq;
-    if (VERBOSE) fprintf (stderr, "dchisq: %f, Ndof: %d\n", dchisq, Npts - Npar);
   }  
 
-  /** create output image (keep in sky : par[6]?) **/
+  if (ShapeVariation) {
+    /* find dChi/dSx and dChi/dSy given by increasing shape terms by 5% */
+    float tp1, tp2, chix, chiy;
+    chix = chiy = 0;
+    if (fitfunc == sgaussTD) {
+      tp1 = par[2];
+      tp2 = par[7];
+      par[2] = par[2]*1.05;
+      par[7] = par[7]*1.05;
+      chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
+      par[2] = tp1;
+      par[7] = tp2;
+
+      tp1 = par[3];
+      tp2 = par[8];
+      par[3] = par[3]*1.05;
+      par[8] = par[8]*1.05;
+      chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
+      par[3] = tp1;
+      par[8] = tp2;
+    }
+    if (fitfunc == pgaussTD) {
+      tp1 = par[2];
+      par[2] = par[2]*1.05;
+      chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
+      par[2] = tp1;
+
+      tp1 = par[3];
+      par[3] = par[3]*1.05;
+      chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
+      par[3] = tp1;
+    }
+    if (fitfunc == SgaussTD) {
+      tp1 = par[0];
+      tp2 = par[3];
+      par[0] = par[0]*1.05;
+      par[3] = par[3]*1.05;
+      chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
+      par[0] = tp1;
+      par[3] = tp2;
+
+      tp1 = par[1];
+      tp2 = par[4];
+      par[1] = par[1]*1.05;
+      par[4] = par[4]*1.05;
+      chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
+      par[1] = tp1;
+      par[4] = tp2;
+    }
+    if (fitfunc == PgaussTD) {
+      tp1 = par[0];
+      par[0] = par[0]*1.05;
+      chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
+      par[0] = tp1;
+
+      tp1 = par[1];
+      par[1] = par[1]*1.05;
+      chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
+      par[1] = tp1;
+    }
+    set_variable ("dChiX", chix/chisq);
+    set_variable ("dChiY", chiy/chisq);
+  }
+
+  /** create output image (keep in sky) **/
   if (Save) {
     Buffer *out;
-    float *Vi, *Vo, vr, vf, *dpar;
+    float *Vi, *Vo, vr, vf, sky;
 
     if ((out = SelectBuffer ("out",   ANYBUFFER, TRUE)) == NULL) return (FALSE);
-    ALLOCATE (dpar, float, 2*Npar);
-  
     free (out[0].header.buffer);
     free (out[0].matrix.buffer);
@@ -181,4 +308,9 @@
     out[0].header.Naxis[1] = 2*ny;
     CreateBuffer (out);
+    if (Npar == 4) {
+      sky = par[3];
+    } else {
+      sky = par[6];
+    }
 
     /* four panels: 1) raw image. 2) fit  3) raw - fit   4) ?? */
@@ -187,40 +319,62 @@
     for (j = 0; j < ny; j++) {
       for (i = 0; i < nx; i++) {
-	vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, dpar);
+	vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL);
 	vr = Vi[(i+sx)+(j+sy)*Nx];
 	Vo[(i   )+(j   )*2*nx] = vr;
 	Vo[(i+nx)+(j   )*2*nx] = vf;
-	Vo[(i   )+(j+ny)*2*nx] = vr - vf + par[6];
-	Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + par[6];
+	Vo[(i   )+(j+ny)*2*nx] = vr - vf + sky;
+	Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + sky;
       }
     }
-    free (dpar);
   }
 
   /* parameter renormalize */
-  par[2] = 2.35*sqrt(2.0) / par[2];
-  par[3] = 2.35*sqrt(2.0) / par[3];
-  if (Npar == 11) {
-    par[8] = 2.35*sqrt(2.0) / par[8];
-    par[9] = 2.35*sqrt(2.0) / par[9];
-  }
-  if (Npar == 10
-) {
-    par[7] = 2.35*sqrt(2.0) / par[7];
-    par[8] = 2.35*sqrt(2.0) / par[8];
-  }
-  for (i = 0; i < Npar; i++) {
-    sprintf (line, "PAR%d", i);
-    set_variable (line, (double) par[i]);
-    if (VERBOSE) fprintf (stderr, "%f ", par[i]);
-  }
-  if (VERBOSE) fprintf (stderr, "  %f\n", chisq/Npts);
-  set_variable ("ChiSq", (double) (chisq/Npts));
-
-  free (par);
+  if (Npar >= 7) {
+    set_variable ("Xg",  par[0]);
+    set_variable ("Yg",  par[1]);
+    set_variable ("SXg", 2.35 * sqrt(2.0) / par[2]);
+    set_variable ("SYg", 2.35 * sqrt(2.0) / par[3]);
+    set_variable ("SXYg", par[4]);
+    set_variable ("Zpk", par[5]);
+    set_variable ("Sg",  par[6]);
+  }
+  if (Npar == 9) {
+    set_variable ("SXf", par[7]);
+    set_variable ("SYf", par[8]);
+  }
+  if (Npar == 10) {
+    set_variable ("SXf", 2.35 * sqrt(2.0) / par[7]);
+    set_variable ("SYf", 2.35 * sqrt(2.0) / par[8]);
+    set_variable ("SXYf", par[9]);
+  }
+  if (fitfunc == PgaussTD) {
+    set_variable ("Xg",  par[0]);
+    set_variable ("Yg",  par[1]);
+    set_variable ("Zpk", par[2]);
+    set_variable ("Sg",  par[3]);
+    set_variable ("SXg", 2.35 * sqrt(2.0) / fpar[0]);
+    set_variable ("SYg", 2.35 * sqrt(2.0) / fpar[1]);
+    set_variable ("SXYg", fpar[2]);
+  }
+  if (fitfunc == SgaussTD) {
+    set_variable ("Xg",  par[0]);
+    set_variable ("Yg",  par[1]);
+    set_variable ("Zpk", par[2]);
+    set_variable ("Sg",  par[3]);
+    set_variable ("SXg", 2.35 * sqrt(2.0) / fpar[0]);
+    set_variable ("SYg", 2.35 * sqrt(2.0) / fpar[1]);
+    set_variable ("SXYg", fpar[2]);
+    set_variable ("SXf", 2.35 * sqrt(2.0) / fpar[3]);
+    set_variable ("SYf", 2.35 * sqrt(2.0) / fpar[4]);
+    set_variable ("SXYf", fpar[5]);
+  }    
+  set_variable ("ChiSq", chisq/(Npts - Npar));
+
   free (x);
   free (y);
   free (z);
   free (dz);
+  free (par);
+  if (Nfpar > 0) free (fpar);
 
   mrq2dfree (Npar);
@@ -232,5 +386,5 @@
 }
 
-/* pars: x, y, sx, sy, sxy, I, sky */
+/* real 2D gaussian -- x, y, sx, sy, sxy, I, sky */
 float fgaussTD (float x, float y, float *par, int Npar, float *dpar) {
 
@@ -249,17 +403,17 @@
   f = q + par[6];
 
-  dpar[0] = q*(2*px*par[2] + par[4]*Y);
-  dpar[1] = q*(2*py*par[3] + par[4]*X);
-  dpar[2] = -2*q*px*X;
-  dpar[3] = -2*q*py*Y;
-  dpar[4] = -q*X*Y;
-  dpar[5] = +r;
-  dpar[6] = +1;
-
-  return (f);
-}
-
-/***** fix the derivatives.  q should not be in function */
-/* pars: x, y, sx, sy, sxy, I, sky */
+  if (dpar != NULL) {
+    dpar[0] = q*(2*px*par[2] + par[4]*Y);
+    dpar[1] = q*(2*py*par[3] + par[4]*X);
+    dpar[2] = -2*q*px*X;
+    dpar[3] = -2*q*py*Y;
+    dpar[4] = -q*X*Y;
+    dpar[5] = +r;
+    dpar[6] = +1;
+  }
+  return (f);
+}
+
+/* pseudo 2D gaussian -- x, y, sx, sy, sxy, I, sky */
 float pgaussTD (float x, float y, float *par, int Npar, float *dpar) {
 
@@ -279,16 +433,46 @@
   /* note difference from gaussian: q = par[5]*r */
 
-  dpar[0] = q*(2*px*par[2] + par[4]*Y);
-  dpar[1] = q*(2*py*par[3] + par[4]*X);
-  dpar[2] = -2*q*px*X;
-  dpar[3] = -2*q*py*Y;
-  dpar[4] = -q*X*Y;
-  dpar[5] = +r;
-  dpar[6] = +1;
-
-  return (f);
-}
-
-/* pars: x, y, sx, sy, sxy, I, sky, f1, f2 */
+  if (dpar != NULL) {
+    dpar[0] = q*(2*px*par[2] + par[4]*Y);
+    dpar[1] = q*(2*py*par[3] + par[4]*X);
+    dpar[2] = -2*q*px*X;
+    dpar[3] = -2*q*py*Y;
+    dpar[4] = -q*X*Y;
+    dpar[5] = +r;
+    dpar[6] = +1;
+  }
+  return (f);
+}
+
+/* pseudo 2D gaussian -- x, y, (sx), (sy), (sxy), I, sky */
+float PgaussTD (float x, float y, float *par, int Npar, float *dpar) {
+
+  float X, Y, px, py;
+  float z, r, q, f;
+
+  /* par -> fpar: (2,0), (3,1), (4,2) */
+
+  X = x - par[0];
+  Y = y - par[1];
+  
+  px = fpar[0]*X;
+  py = fpar[1]*Y;
+
+  z = 0.5*SQ(px) + 0.5*SQ(py) + fpar[2]*X*Y;
+  r = 1.0 / (1 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */
+  f = par[2]*r + par[3];
+  q = par[2]*r*r*(1 + z + 0.5*z*z);
+  /* note difference from gaussian: q = par[5]*r */
+
+  if (dpar != NULL) {
+    dpar[0] = q*(2*px*fpar[0] + fpar[2]*Y);
+    dpar[1] = q*(2*py*fpar[1] + fpar[2]*X);
+    dpar[2] = +r;
+    dpar[3] = +1;
+  }
+  return (f);
+}
+
+/* pseudo 2D gaussian with floating 2nd and 3rd order terms -- x, y, sx, sy, sxy, I, sky, f1, f2 */
 float vgaussTD (float x, float y, float *par, int Npar, float *dpar) {
 
@@ -309,18 +493,321 @@
   /* note difference from gaussian: q = par[5]*r */
 
-  dpar[0] = q*(2*px*par[2] + par[4]*Y);
-  dpar[1] = q*(2*py*par[3] + par[4]*X);
-  dpar[2] = -2*q*px*X;
-  dpar[3] = -2*q*py*Y;
-  dpar[4] = -q*X*Y;
-  dpar[5] = +r;
-  dpar[6] = +1;
-  dpar[7] = -100*par[5]*r*r*k;
-  dpar[8] = -100*par[5]*r*r*par[7]*(z*z*z)/6;
-
-  return (f);
-}
-
+  if (dpar != NULL) {
+    dpar[0] = q*(2*px*par[2] + par[4]*Y);
+    dpar[1] = q*(2*py*par[3] + par[4]*X);
+    dpar[2] = -2*q*px*X;
+    dpar[3] = -2*q*py*Y;
+    dpar[4] = -q*X*Y;
+    dpar[5] = +r;
+    dpar[6] = +1;
+    dpar[7] = -100*par[5]*r*r*k;
+    dpar[8] = -100*par[5]*r*r*par[7]*(z*z*z)/6;
+  }
+  return (f);
+}
+
+/* two components: (1 + z_1 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
+float sgaussTD (float x, float y, float *par, int Npar, float *dpar) {
+
+  float X, Y, px1, py1, px2, py2;
+  float z1, z2, r, q1, q2, f, f1, f2;
+
+  X = x - par[0];
+  Y = y - par[1];
+  
+  px1 = par[2]*X;
+  py1 = par[3]*Y;
+  px2 = par[7]*X;
+  py2 = par[8]*Y;
+
+  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
+  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
+
+  r = 1.0 / (1 + z1 + pow(z2,Npow));
+  f = par[5]*r + par[6];
+
+  q1 = par[5]*SQ(r);
+  q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
+
+  if (dpar != NULL) {
+    dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
+    dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
+
+    /* these fudge factors impede the growth of par[2] beyond par[7] */
+    f1 = fabs(par[7]) / fabs(par[2]);
+    f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
+    dpar[2] = -2*q1*px1*X*f2;
+
+    f1 = fabs(par[8]) / fabs(par[3]);
+    f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
+    dpar[3] = -2*q1*py1*Y*f2;
+
+    dpar[4] = -q1*X*Y;
+    dpar[5] = +r;
+    dpar[6] = +1;
+    dpar[7] = -2*q2*px2*X;
+    dpar[8] = -2*q2*py2*Y;
+    dpar[9] = -q2*X*Y;
+  }
+  return (f);
+}
+
+/* two components: (1 + z_1 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
+float SgaussTD (float x, float y, float *par, int Npar, float *dpar) {
+
+  float X, Y, px1, py1, px2, py2;
+  float z1, z2, r, q1, q2, f, f1, f2;
+
+  X = x - par[0];
+  Y = y - par[1];
+  
+  px1 = fpar[0]*X;
+  py1 = fpar[1]*Y;
+  px2 = fpar[3]*X;
+  py2 = fpar[4]*Y;
+
+  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + fpar[2]*X*Y;
+  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + fpar[5]*X*Y;
+
+  r = 1.0 / (1 + z1 + pow(z2,Npow));
+  f = par[2]*r + par[3];
+
+  q1 = par[2]*SQ(r);
+  q2 = par[2]*SQ(r)*Npow*pow(z2,(Npow-1));
+
+  if (dpar != NULL) {
+    dpar[0] = q1*(2*px1*fpar[0] + fpar[2]*Y) + q2*(2*px2*fpar[3] + fpar[5]*Y);
+    dpar[1] = q1*(2*py1*fpar[1] + fpar[2]*X) + q2*(2*py2*fpar[4] + fpar[5]*X);
+    dpar[2] = +r;
+    dpar[3] = +1;
+  }
+  return (f);
+}
+
+/* two components: (1 + z_1 + 0.5*z_1^2 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
+float rgaussTD (float x, float y, float *par, int Npar, float *dpar) {
+
+  float X, Y, px1, py1, px2, py2;
+  float z1, z2, r, q1, q2, f;
+
+  X = x - par[0];
+  Y = y - par[1];
+  
+  px1 = par[2]*X;
+  py1 = par[3]*Y;
+  px2 = par[7]*X;
+  py2 = par[8]*Y;
+
+  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
+  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
+
+  r = 1.0 / (1 + z1 + 0.5*SQ(z1)+ pow(z2,Npow));
+  f = par[5]*r + par[6];
+
+  q1 = par[5]*SQ(r)*(1 + z1);
+  q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
+
+  if (dpar != NULL) {
+    dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
+    dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
+    dpar[2] = -2*q1*px1*X;
+    dpar[3] = -2*q1*py1*Y;
+    dpar[4] = -q1*X*Y;
+    dpar[5] = +r;
+    dpar[6] = +1;
+    dpar[7] = -2*q2*px2*X;
+    dpar[8] = -2*q2*py2*Y;
+    dpar[9] = -q2*X*Y;
+  }
+  return (f);
+}
+
+/* two components: (1 + z_1^M + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
+float tgaussTD (float x, float y, float *par, int Npar, float *dpar) {
+
+  float X, Y, px1, py1, px2, py2;
+  float z1, z2, r, q1, q2, f;
+
+  X = x - par[0];
+  Y = y - par[1];
+  
+  px1 = par[2]*X;
+  py1 = par[3]*Y;
+  px2 = par[7]*X;
+  py2 = par[8]*Y;
+
+  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
+  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
+
+  r = 1.0 / (1 + pow(z1,Npin) + pow(z2,Npow));
+  f = par[5]*r + par[6];
+
+  q1 = par[5]*SQ(r)*Npin*pow(z1,(Npin-1));
+  q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
+
+  if (dpar != NULL) {
+    dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
+    dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
+    dpar[2] = -2*q1*px1*X*2;
+    dpar[3] = -2*q1*py1*Y*2;
+    dpar[4] = -q1*X*Y;
+    dpar[5] = +r;
+    dpar[6] = +1;
+    dpar[7] = -2*q2*px2*X;
+    dpar[8] = -2*q2*py2*Y;
+    dpar[9] = -q2*X*Y;
+  }
+  return (f);
+}
+
+/* alternate for exp(-z) */
+/* E = 1.0 / (1 + f1*Z + f2*(0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
+/* E = 1.0 / (1 +    Z +    (0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
+
+# if (0)
 /* pars: x, y, sx, sy, sxy, I, sky */
+float fgalaxyTD (float x, float y, float *par, int Npar, float *dpar) {
+
+  float X, Y, Z, E, F, q, R, f, p2, p3;
+
+  X = x - par[0];
+  Y = y - par[1];
+  
+  p2 = X / par[2];
+  p3 = Y / par[3];
+
+  Z = sqrt (0.5*p2*X + X*Y*par[4] + 0.5*p3*Y);                 /* R */
+  E = 1.0 / (1 + Z);   
+
+  q = par[5] * E;
+  R = q*E;
+  F = 0.5 / Z;
+  
+  f = q + par[6];
+
+  dpar[0] = F*R*(p2 + par[4]*Y);
+  dpar[1] = F*R*(p3 + par[4]*X);
+  dpar[2] = F*0.5*R*p2*p2;
+  dpar[3] = F*0.5*R*p3*p3;
+  dpar[4] = -R*X*Y*F;
+    
+  dpar[5] = E;
+  dpar[6] = 1;
+  return (f);
+}
+
+/* pars: x, y, sx, sy, sxy, I, sky */
+float fbarTD (float x, float y, float *par, int Npar, float *dpar) {
+
+  float X, Y, Z, E, F, q, R, f, p2, p3;
+
+  X = x - par[0];
+  Y = y - par[1];
+  
+  p2 = X / par[2];
+  p3 = Y / par[3];
+
+  Z = 0.5*p2*X + X*Y*par[4] + 0.5*p3*Y;                 /* R */
+  E = 1.0 / (1 + Z*Z*Z);   
+
+  q = par[5] * E;
+  F = 3*Z*Z;
+  R = q*E*F;
+  
+  f = q + par[6];
+
+  dpar[0] = R*(p2 + par[4]*Y);
+  dpar[1] = R*(p3 + par[4]*X);
+  dpar[2] = 0.5*R*p2*p2;
+  dpar[3] = 0.5*R*p3*p3;
+  dpar[4] = -R*X*Y;
+    
+  dpar[5] = E;
+  dpar[6] = 1;
+  return (f);
+}
+
+/* convert from x,y to major,minor */
+void fix_ellipsegauss_pars (float *par, int Npar) {
+
+  float p2, p4, angle, t1, t2, tmp, area;
+
+  /* par[0], par[1] = Xo, Yo - stay the same */
+
+  p2 = 1/par[2];
+  p4 = 1/par[3];
+
+  angle = 0.5 * atan2 (-2*par[4], p4 - p2); 
+
+  tmp = sqrt (SQ(p2 - p4) + 4*SQ(par[4]));
+  t1 = (p2 + p4 + tmp) / 2;
+  t2 = t1 - tmp;
+
+  par[2] = 2.35482*sqrt(1/t2);
+  par[3] = 2.35482*sqrt(1/t1);
+  par[4] = DEG_RAD * angle;
+
+  area = 2*M_PI/sqrt(t1*t2);
+
+  par[5] *= area;
+
+}
+# endif
+
+/***  options for later
+
+  Subtract = FALSE;
+  if ((N = get_argument (argc, argv, "-sub"))) {
+    remove_argument (N, &argc, argv);
+    Subtract  = TRUE;
+  }
+
+  DFact = 1;
+  if ((N = get_argument (argc, argv, "-D"))) {
+    remove_argument (N, &argc, argv);
+    DFact  = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  fitfunc = fgaussTD;
+  if ((N = get_argument (argc, argv, "-gal"))) {
+    remove_argument (N, &argc, argv);
+    fitfunc = fgalaxyTD; 
+  }
+  if ((N = get_argument (argc, argv, "-bar"))) {
+    remove_argument (N, &argc, argv);
+    fitfunc = fbarTD; 
+  }
+
+
+  f1 = 1;
+  if ((c = get_variable ("BETA1")) != (char *) NULL) f1 = atof (c);
+
+  f2 = 1;
+  if ((c = get_variable ("BETA2")) != (char *) NULL) f2 = atof (c);
+
+  if (Subtract) {
+    tmpsky = par[6];
+    par[6] = 0;
+    for (N = j = 0; j < ny; j++) {
+      V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx; 
+      for (i = 0; i < nx; i++, V++, N++) {
+	dx = i + sx;
+	dy = j + sy;
+	*V -= fitfunc (dx, dy, par, Npar, (float *) NULL);
+      }
+    }
+    par[6] = tmpsky;
+  }
+
+***/
+
+# if (0)
+
+/* these two tests were not very succcessful.  the first did not model the shape well because 
+   it could not match the change in roundness with radius.  the second did not work because the 
+   parameters were degenerate (amplitude and slope of second component) */
+
+/* test: fixed, non-integer higher-order term -- x, y, sx, sy, sxy, I, sky, f1, f2 */
 float qgaussTD (float x, float y, float *par, int Npar, float *dpar) {
 
@@ -354,5 +841,5 @@
 }
 
-/* pars: x, y, sx, sy, sxy, I, sky */
+/* test: two component model: inner pseudo gaussian with outer z^1.75 x, y, sx, sy, sxy, I, sky */
 float rgaussTD (float x, float y, float *par, int Npar, float *dpar) {
 
@@ -394,180 +881,36 @@
 }
 
-/* pars: x, y, sx, sy, sxy, I, sky */
-float sgaussTD (float x, float y, float *par, int Npar, float *dpar) {
-
-  float X, Y, px1, py1, px2, py2;
-  float z1, z2, r, q1, q2, f;
-
-  X = x - par[0];
-  Y = y - par[1];
-  
-  px1 = par[2]*X;
-  py1 = par[3]*Y;
-  px2 = par[7]*X;
-  py2 = par[8]*Y;
-
-  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
-  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
-
-  r = 1.0 / (1 + z1 + pow(z2,Npow));
-  f = par[5]*r + par[6];
-
-  q1 = par[5]*SQ(r);
-  q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
-
-  dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
-  dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
-  dpar[2] = -2*q1*px1*X;
-  dpar[3] = -2*q1*py1*Y;
-  dpar[4] = -q1*X*Y;
-  dpar[5] = +r;
-  dpar[6] = +1;
-  dpar[7] = -2*q2*px2*X;
-  dpar[8] = -2*q2*py2*Y;
-  dpar[9] = -q2*X*Y;
-
-  return (f);
-}
-
-/* alternate for exp(-z) */
-/* E = 1.0 / (1 + f1*Z + f2*(0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
-/* E = 1.0 / (1 +    Z +    (0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
-
-# if (0)
-/* pars: x, y, sx, sy, sxy, I, sky */
-float fgalaxyTD (float x, float y, float *par, int Npar, float *dpar) {
-
-  float X, Y, Z, E, F, q, R, f, p2, p3;
-
-  X = x - par[0];
-  Y = y - par[1];
-  
-  p2 = X / par[2];
-  p3 = Y / par[3];
-
-  Z = sqrt (0.5*p2*X + X*Y*par[4] + 0.5*p3*Y);                 /* R */
-  E = 1.0 / (1 + Z);   
-
-  q = par[5] * E;
-  R = q*E;
-  F = 0.5 / Z;
-  
-  f = q + par[6];
-
-  dpar[0] = F*R*(p2 + par[4]*Y);
-  dpar[1] = F*R*(p3 + par[4]*X);
-  dpar[2] = F*0.5*R*p2*p2;
-  dpar[3] = F*0.5*R*p3*p3;
-  dpar[4] = -R*X*Y*F;
-    
-  dpar[5] = E;
-  dpar[6] = 1;
-  return (f);
-}
-
-/* pars: x, y, sx, sy, sxy, I, sky */
-float fbarTD (float x, float y, float *par, int Npar, float *dpar) {
-
-  float X, Y, Z, E, F, q, R, f, p2, p3;
-
-  X = x - par[0];
-  Y = y - par[1];
-  
-  p2 = X / par[2];
-  p3 = Y / par[3];
-
-  Z = 0.5*p2*X + X*Y*par[4] + 0.5*p3*Y;                 /* R */
-  E = 1.0 / (1 + Z*Z*Z);   
-
-  q = par[5] * E;
-  F = 3*Z*Z;
-  R = q*E*F;
-  
-  f = q + par[6];
-
-  dpar[0] = R*(p2 + par[4]*Y);
-  dpar[1] = R*(p3 + par[4]*X);
-  dpar[2] = 0.5*R*p2*p2;
-  dpar[3] = 0.5*R*p3*p3;
-  dpar[4] = -R*X*Y;
-    
-  dpar[5] = E;
-  dpar[6] = 1;
-  return (f);
-}
-
-/* convert from x,y to major,minor */
-void fix_ellipsegauss_pars (float *par, int Npar) {
-
-  float p2, p4, angle, t1, t2, tmp, area;
-
-  /* par[0], par[1] = Xo, Yo - stay the same */
-
-  p2 = 1/par[2];
-  p4 = 1/par[3];
-
-  angle = 0.5 * atan2 (-2*par[4], p4 - p2); 
-
-  tmp = sqrt (SQ(p2 - p4) + 4*SQ(par[4]));
-  t1 = (p2 + p4 + tmp) / 2;
-  t2 = t1 - tmp;
-
-  par[2] = 2.35482*sqrt(1/t2);
-  par[3] = 2.35482*sqrt(1/t1);
-  par[4] = DEG_RAD * angle;
-
-  area = 2*M_PI/sqrt(t1*t2);
-
-  par[5] *= area;
-
-}
 # endif
 
-/***  options for later
-
-  Subtract = FALSE;
-  if ((N = get_argument (argc, argv, "-sub"))) {
-    remove_argument (N, &argc, argv);
-    Subtract  = TRUE;
-  }
-
-  DFact = 1;
-  if ((N = get_argument (argc, argv, "-D"))) {
-    remove_argument (N, &argc, argv);
-    DFact  = atof(argv[N]);
-    remove_argument (N, &argc, argv);
-  }
-
-  fitfunc = fgaussTD;
-  if ((N = get_argument (argc, argv, "-gal"))) {
-    remove_argument (N, &argc, argv);
-    fitfunc = fgalaxyTD; 
-  }
-  if ((N = get_argument (argc, argv, "-bar"))) {
-    remove_argument (N, &argc, argv);
-    fitfunc = fbarTD; 
-  }
-
-
-  f1 = 1;
-  if ((c = get_variable ("BETA1")) != (char *) NULL) f1 = atof (c);
-
-  f2 = 1;
-  if ((c = get_variable ("BETA2")) != (char *) NULL) f2 = atof (c);
-
-  if (Subtract) {
-    tmpsky = par[6];
-    par[6] = 0;
-    for (N = j = 0; j < ny; j++) {
-      V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx; 
-      for (i = 0; i < nx; i++, V++, N++) {
-	dx = i + sx;
-	dy = j + sy;
-	*V -= fitfunc (dx, dy, par, Npar, (float *) NULL);
-      }
-    }
-    par[6] = tmpsky;
-  }
-
-***/
+  /* forcing values to have a rational range
+  ALLOCATE (parmin, float, Npar);
+  ALLOCATE (parmax, float, Npar);
+  bzero (parmin, Npar*sizeof(float));
+  bzero (parmax, Npar*sizeof(float));
+  parmin[0] = parmin[1] = 0;
+  parmax[0] = buf[0].matrix.Naxis[0];
+  parmax[1] = buf[0].matrix.Naxis[1];
+
+  parmin[2] = parmin[3] = 0.01;
+  parmax[2] = parmax[3] = 100.0;
+  parmin[4] = -1000;
+  parmax[4] = -1000;
+  
+  parmin[5] = 1;
+  parmax[5] = 1e5;
+
+  parmin[6] = 0.0;
+  parmax[6] = 1e5;
+
+  if (Npar == 9) {
+    parmin[7] = parmin[8] = 0.01;
+    parmax[7] = parmax[8] = 10.0;
+  }
+  if (Npar == 10) {
+    parmin[7] = parmin[8] = 0.01;
+    parmax[7] = parmax[8] = 10.0;
+    parmin[9] = -1000;
+    parmax[9] = -1000;
+  }
+  */
+
Index: /trunk/Ohana/src/opihi/cmd.data/gridify.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/gridify.c	(revision 2937)
+++ /trunk/Ohana/src/opihi/cmd.data/gridify.c	(revision 2938)
@@ -3,5 +3,5 @@
 int gridify (int argc, char **argv) {
 
-  int i, Nx, Ny, Xb, Yb;
+  int i, Nx, Ny, Xb, Yb, Normalize, N;
   float Xmin, Xmax, dX, Ymin, Ymax, dY;
   float *buf, *val, *x, *y, *z;
@@ -9,4 +9,10 @@
   Buffer *bf;
   Vector *vx, *vy, *vz;
+
+  Normalize = TRUE;
+  if ((N = get_argument (argc, argv, "-raw"))) {
+    remove_argument (N, &argc, argv);
+    Normalize = FALSE;
+  }
 
   if (argc != 11) {
@@ -66,5 +72,5 @@
     if (Xb >= Nx) continue;
     if (Yb >= Ny) continue;
-    val[Xb + Yb*Nx] = *z;
+    val[Xb + Yb*Nx] += *z;
     Nval[Xb + Yb*Nx]++;
   }
@@ -72,9 +78,13 @@
   buf = (float *) bf[0].matrix.buffer;
   for (i = 0; i < Nx*Ny; i++) {
-    if (Nval[i] == 0) {
-      buf[i] = 0;
-      continue;
+    if (Normalize) {
+      if (Nval[i] == 0) {
+	buf[i] = 0;
+	continue;
+      }
+      buf[i] = val[i] / Nval[i];
+    } else {
+      buf[i] = val[i];
     }
-    buf[i] = val[i] / Nval[i];
   }
 
Index: /trunk/Ohana/src/opihi/cmd.data/imsmooth.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/imsmooth.c	(revision 2937)
+++ /trunk/Ohana/src/opihi/cmd.data/imsmooth.c	(revision 2938)
@@ -1,12 +1,18 @@
 # include "data.h"
-# define NSIGMA 3
 
 int imsmooth (int argc, char **argv) {
   
-  int i, j, n, Nx, Ny, Ns, Ngauss;
+  int i, j, n, N, Nx, Ny, Ns, Ngauss;
   float *vi, *vo, *gauss, *gaussnorm;
-  float g, s, sigma;
+  float g, s, sigma, Nsigma;
   Buffer *in;
   float *temp;
+
+  Nsigma = 3;
+  if ((N = get_argument (argc, argv, "-Nsigma"))) {
+    remove_argument (N, &argc, argv);
+    Nsigma = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
 
   if (argc != 3) {
@@ -23,5 +29,5 @@
 
   /* build a 1D gaussian */
-  Ns = (int) (NSIGMA*sigma + 0.5);
+  Ns = (int) (Nsigma*sigma + 0.5);
   Ngauss = 2*Ns + 1;
   ALLOCATE (gaussnorm, float, Ngauss);
Index: /trunk/Ohana/src/opihi/lib.data/mrq2dmin.c
===================================================================
--- /trunk/Ohana/src/opihi/lib.data/mrq2dmin.c	(revision 2937)
+++ /trunk/Ohana/src/opihi/lib.data/mrq2dmin.c	(revision 2938)
@@ -7,8 +7,13 @@
 */
 
+# define VERY_VERBOSE 0
+
 static float **alpha, **talpha;
 static float **beta, **tbeta;
 static float *partry, *dyda;
 static float ochisq, lambda;
+
+static float *parmin = NULL;
+static float *parmax = NULL;
 
 float mrq2dcof (float *x, float *t, float *y, float *dy, int Npts, 
@@ -41,6 +46,31 @@
       ta[k][j] = ta[j][k];
       
+# if (VERY_VERBOSE)
+  for (j = 0; j < Npar; j++) {
+    for (k = 0; k < Npar; k++) {
+      fprintf (stderr, "%9.3e  ", ta[j][k]);
+    }
+    fprintf (stderr, "    :   %9.3e  ", tb[j][0]);
+    fprintf (stderr, "\n");
+  }
+# endif
+
   return (chisq);
 
+}
+
+float mrq2dchi (float *x, float *t, float *y, float *dy, int Npts, 
+		float *par, int Npar, 
+		float (funcs)(float, float, float *, int, float *)) {
+
+  int k, j, i;
+  float ydiff, wt, chisq;
+
+  chisq = 0.0;
+  for (i = 0; i < Npts; i++) {
+    ydiff = funcs (x[i], t[i], par, Npar, dyda) - y[i];
+    chisq += SQ(ydiff) * dy[i];
+  }
+  return (chisq);
 }
 
@@ -65,5 +95,11 @@
   }
 
-  for (j = 0; j < Npar; j++) partry[j] = par[j] - tbeta[j][0];
+  for (j = 0; j < Npar; j++) {
+    partry[j] = par[j] - tbeta[j][0];
+    /*
+    if (parmin != NULL) partry[j] = MAX (parmin[j], partry[j]);
+    if (parmax != NULL) partry[j] = MIN (parmax[j], partry[j]);
+    */
+  }
 
   chisq = mrq2dcof (x, t, y, dy, Npts, partry, Npar, talpha, tbeta, funcs);
@@ -93,4 +129,17 @@
   return (chisq);
 
+}
+
+int mrq2dlimits (float *pmin, float *pmax, int Npar) {
+
+  int i;
+
+  ALLOCATE (parmin, float, Npar);
+  ALLOCATE (parmax, float, Npar);
+  for (i = 0; i < Npar; i++) {
+    parmin[i] = pmin[i];
+    parmax[i] = pmax[i];
+  }
+  return (TRUE);
 }
 
@@ -156,3 +205,6 @@
   free (dyda);
 
-}
+  if (parmin != NULL) free (parmin);
+  if (parmax != NULL) free (parmax);
+
+}
Index: /trunk/Ohana/src/opihi/mana/findpeaks.c
===================================================================
--- /trunk/Ohana/src/opihi/mana/findpeaks.c	(revision 2937)
+++ /trunk/Ohana/src/opihi/mana/findpeaks.c	(revision 2938)
@@ -7,7 +7,7 @@
   int Npeak, NPEAK, Npeaks;
   float *v;
-  int *peaks, *keep, *xp, *yp;
+  int *peaks, *keep, *xp, *yp, *zp;
   float threshold, vt, vo;
-  Vector *vecx, *vecy;
+  Vector *vecx, *vecy, *vecz;
   Buffer *buf;
 
@@ -19,4 +19,5 @@
   if ((vecx = SelectVector ("xp", ANYVECTOR, TRUE)) == NULL) return (FALSE);
   if ((vecy = SelectVector ("yp", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((vecz = SelectVector ("zp", ANYVECTOR, TRUE)) == NULL) return (FALSE);
 
   Nx = buf[0].matrix.Naxis[0];
@@ -27,4 +28,5 @@
   ALLOCATE (xp, int, NPEAK);
   ALLOCATE (yp, int, NPEAK);
+  ALLOCATE (zp, int, NPEAK);
 
   /* find peaks for each row */
@@ -37,8 +39,10 @@
       REALLOCATE (xp, int, NPEAK);
       REALLOCATE (yp, int, NPEAK);
+      REALLOCATE (zp, int, NPEAK);
     }
     for (i = 0; i < Npeaks; i++) {
       xp[Npeak + i] = peaks[i];
       yp[Npeak + i] = j;
+      zp[Npeak + i] = v[peaks[i] + j*Nx];
     }
     Npeak += Npeaks;
@@ -74,4 +78,5 @@
   REALLOCATE (vecx[0].elements, float, MAX (Npeak, 1));
   REALLOCATE (vecy[0].elements, float, MAX (Npeak, 1));
+  REALLOCATE (vecz[0].elements, float, MAX (Npeak, 1));
   /* eliminate non-local peaks */
   for (N = n = 0; n < Npeak; n++) {
@@ -79,10 +84,16 @@
     vecx[0].elements[N] = xp[n];
     vecy[0].elements[N] = yp[n];
+    vecz[0].elements[N] = zp[n];
     N ++;
   }
+  free (xp);
+  free (yp);
+  free (zp);
+  free (keep);
 
   REALLOCATE (vecx[0].elements, float, MAX (N, 1));
   REALLOCATE (vecy[0].elements, float, MAX (N, 1));
-  vecx[0].Nelements = vecy[0].Nelements = N;
+  REALLOCATE (vecz[0].elements, float, MAX (N, 1));
+  vecx[0].Nelements = vecy[0].Nelements = vecz[0].Nelements = N;
 
   return (TRUE);
@@ -93,2 +104,3 @@
 }
 
+
Index: /trunk/Ohana/src/opihi/mana/findrowpeaks.c
===================================================================
--- /trunk/Ohana/src/opihi/mana/findrowpeaks.c	(revision 2937)
+++ /trunk/Ohana/src/opihi/mana/findrowpeaks.c	(revision 2938)
@@ -3,16 +3,7 @@
 int *findrowpeaks (float *row, int Nrow, float threshold, int *npeaks) {
   
-  int i, rising;
+  int i;
   int Npeaks, NPEAKS;
   int *peaks;
-
-# if (0)
-  ALLOCATE (peaks, int, 3);
-  peaks[0] = 20;
-  peaks[1] = 30;
-  peaks[2] = 40;
-  *npeaks = 3;
-  return (peaks);
-# endif
 
   Npeaks = 0;
@@ -25,31 +16,29 @@
   }
     
-  rising = (row[1] > row[0]);
+  /* special case for first pixel in row */
+  if ((row[0] > row[1]) && (row[0] > threshold)) {
+    peaks[Npeaks] = 0;
+    Npeaks ++;
+  }    
 
-  for (i = 2; i < Nrow; i++) {
-    if (row[i] < threshold) {
-      rising = TRUE;
-      continue;
-    }
-    if (rising) {
-      if (row[i] < row[i-1]) {
-	rising = FALSE;
-	peaks[Npeaks] = i-1;
-	Npeaks ++;
-	if (Npeaks >= NPEAKS) {
-	  NPEAKS += 100;
-	  REALLOCATE (peaks, int, NPEAKS);
-	  fprintf (stderr, "#");
-	}
-      }
-      continue;
-    }
-    if (!rising) {
-      if (row[i] > row[i-1]) {
-	rising = TRUE;
-      }
-      continue;
+  for (i = 1; i < Nrow - 1; i++) {
+    if (row[i] < threshold) continue;
+    if (row[i] < row[i-1]) continue;
+    if (row[i] <= row[i+1]) continue;
+
+    peaks[Npeaks] = i;
+    Npeaks ++;
+    if (Npeaks >= NPEAKS) {
+      NPEAKS += 100;
+      REALLOCATE (peaks, int, NPEAKS);
     }
   }      
+
+  /* special case for last pixel in row */
+  if ((row[Nrow-1] >= row[Nrow-2]) && (row[Nrow-1] > threshold)) {
+    peaks[Npeaks] = Nrow-1;
+    Npeaks ++;
+  }    
+
   *npeaks = Npeaks;
   return (peaks);
Index: /trunk/Ohana/src/opihi/mana/rawstars.c
===================================================================
--- /trunk/Ohana/src/opihi/mana/rawstars.c	(revision 2937)
+++ /trunk/Ohana/src/opihi/mana/rawstars.c	(revision 2938)
@@ -7,5 +7,5 @@
   double Raper, Rinner, Router;
   Vector *xp, *yp;
-  Vector *xc, *yc, *sx, *sy, *sxy, *zs, *zp, *sk;
+  Vector *xc, *yc, *sx, *sy, *sxy, *zs, *zc, *sk;
   Buffer *buff;
 
@@ -44,10 +44,10 @@
   if ((xc = SelectVector ("xc", ANYVECTOR, TRUE)) == NULL) return (FALSE);
   if ((yc = SelectVector ("yc", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((zc = SelectVector ("zc", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((zs = SelectVector ("zs", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((sk = SelectVector ("sk", ANYVECTOR, TRUE)) == NULL) return (FALSE);
   if ((sx = SelectVector ("sx", ANYVECTOR, TRUE)) == NULL) return (FALSE);
   if ((sy = SelectVector ("sy", ANYVECTOR, TRUE)) == NULL) return (FALSE);
   if ((sxy = SelectVector ("sxy", ANYVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((zs = SelectVector ("zs", ANYVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((zp = SelectVector ("zp", ANYVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((sk = SelectVector ("sk", ANYVECTOR, TRUE)) == NULL) return (FALSE);
 
   Nx = buff[0].matrix.Naxis[0];
@@ -61,8 +61,8 @@
   REALLOCATE (sxy[0].elements, float, Np);
   REALLOCATE (zs[0].elements, float, Np);
-  REALLOCATE (zp[0].elements, float, Np);
+  REALLOCATE (zc[0].elements, float, Np);
   REALLOCATE (sk[0].elements, float, Np);
   xc[0].Nelements = yc[0].Nelements = sx[0].Nelements = Np;
-  sy[0].Nelements = zs[0].Nelements = zp[0].Nelements = Np;
+  sy[0].Nelements = zs[0].Nelements = zc[0].Nelements = Np;
   sxy[0].Nelements = sk[0].Nelements = Np;
 
@@ -83,5 +83,5 @@
 		    &sxy[0].elements[i], 
 		    &zs[0].elements[i], 
-		    &zp[0].elements[i],
+		    &zc[0].elements[i],
 		    &sk[0].elements[i]);
   }