Index: /trunk/Ohana/src/opihi/cmd.astro/imfit.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 2875)
+++ /trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 2876)
@@ -5,4 +5,8 @@
 float vgaussTD (float, float, float *, int, float *);
 float qgaussTD (float, float, float *, int, float *);
+float rgaussTD (float, float, float *, int, float *);
+float sgaussTD (float, float, float *, int, float *);
+
+static float Npow;
 
 int imfit (int argc, char **argv) {
@@ -10,5 +14,5 @@
   int i, j, N, Npts, Npar, Save, VERBOSE;
   int sx, sy, nx, ny, Nx, Ny;
-  float chisq, ochisq, dchisq;
+  float chisq, ochisq, dchisq, Gain, RDnoise;
   float *x, *y, *z, *dz, *par, *V;
   char line[64];
@@ -22,4 +26,20 @@
   }
 
+  /* Gain in e/DN */
+  Gain = 1.0;
+  if ((N = get_argument (argc, argv, "-gain"))) {
+    remove_argument (N, &argc, argv);
+    Gain = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  /* RD noise in DN */
+  RDnoise = 0.0;
+  if ((N = get_argument (argc, argv, "-rdnoise"))) {
+    remove_argument (N, &argc, argv);
+    RDnoise = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   VERBOSE = FALSE;
   if ((N = get_argument (argc, argv, "-v"))) {
@@ -28,18 +48,37 @@
   }
 
+  /* set fitting function */
   Npar = 7;
   fitfunc = fgaussTD;
-  if ((N = get_argument (argc, argv, "-pgauss"))) {
-    remove_argument (N, &argc, argv);
-    fitfunc = pgaussTD;
-  }
-  if ((N = get_argument (argc, argv, "-qgauss"))) {
-    remove_argument (N, &argc, argv);
-    fitfunc = qgaussTD;
-  }
-  if ((N = get_argument (argc, argv, "-vgauss"))) {
-    remove_argument (N, &argc, argv);
-    fitfunc = vgaussTD;
-    Npar = 9;
+  if ((N = get_argument (argc, argv, "-func"))) {
+    fitfunc = NULL;
+    remove_argument (N, &argc, argv);
+    if (!strcmp(argv[N], "fgauss")) {
+      fitfunc = fgaussTD;
+    }
+    if (!strcmp(argv[N], "pgauss")) {
+      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 (fitfunc == NULL) {
+      fprintf (stderr, "unknown function %s\n", argv[N]);
+      return (FALSE);
+    }
+    remove_argument (N, &argc, argv);
   }
 
@@ -49,4 +88,5 @@
   }
 
+  /* non-optional arguments */
   if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
   sx = atof (argv[2]);
@@ -66,4 +106,5 @@
 
   /* convert array z[x,y] to x[i], y[i], z[i] */
+  N = 0;
   Npts = nx*ny;
   ALLOCATE (x,  float, 2*Npts);
@@ -71,18 +112,19 @@
   ALLOCATE (z,  float, 2*Npts);
   ALLOCATE (dz, float, 2*Npts);
-
-  N = 0;
   for (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++) {
+    for (i = 0; i < nx; i++, V++) {
+      dz[N] = (SQ(RDnoise) + *V/Gain);
+      if (dz[N] <= 0) continue;
+      dz[N] = 1.0 / dz[N];
       x[N] = i + sx;
       y[N] = j + sy;
       z[N] = *V;
-      dz[N] = 1.0 / *V;  /* temporarily no gain (or, gain = 1) */
-    }
-  }
+      N++;
+    }
+  }
+  Npts = N;
 
   /* make initial guesses on parameters */
-  /* fitting X^2 / (2*Sx), Y^2 / (2*Sy) */
   par[0] = atof (get_variable ("Xg"));
   par[1] = atof (get_variable ("Yg"));
@@ -95,4 +137,17 @@
     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"));
+    par[9] = 0.0;
+  }
+
 
   /* run fit routine */
@@ -100,5 +155,5 @@
   dchisq = ochisq;
 
-  for (i = 0; (i < 20) && ((dchisq <= 0.0) || (dchisq > 0.001*(Npts - Npar))); i++) {
+  for (i = 0; (i < 50) && ((dchisq <= 0.0) || (dchisq > 0.001*(Npts - Npar))); i++) {
     chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
     dchisq = ochisq - chisq;
@@ -107,5 +162,5 @@
   }  
 
-  /** create output image **/
+  /** create output image (keep in sky : par[6]?) **/
   if (Save) {
     Buffer *out;
@@ -137,5 +192,5 @@
 	Vo[(i+nx)+(j   )*2*nx] = vf;
 	Vo[(i   )+(j+ny)*2*nx] = vr - vf;
-	Vo[(i+nx)+(j+ny)*2*nx] = 0;
+	Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf);
       }
     }
@@ -146,4 +201,13 @@
   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);
@@ -152,4 +216,5 @@
   }
   if (VERBOSE) fprintf (stderr, "  %f\n", chisq/Npts);
+  set_variable ("ChiSq", (double) (chisq/Npts));
 
   free (par);
@@ -210,8 +275,7 @@
   z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
   r = 1.0 / (1 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */
-  q = par[5]*r;
-  k = q*r*(1 + z + 0.5*z*z);
+  f = par[5]*r + par[6];
+  q = par[5]*r*r*(1 + z + 0.5*z*z);
   /* note difference from gaussian: q = par[5]*r */
-  f = q + par[6];
 
   dpar[0] = q*(2*px*par[2] + par[4]*Y);
@@ -230,5 +294,5 @@
 
   float X, Y, px, py;
-  float z, r, q, f;
+  float z, r, q, f, k;
 
   X = x - par[0];
@@ -239,8 +303,9 @@
 
   z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
-  r = 1.0 / (1 + z + par[7]*0.5*z*z*(1 + par[8]*z/3)); /* ~ exp (-Z) */
-  q = par[5]*r*r*(1 + z + 0.5*z*z);
+  k = 0.5*z*z*(1 + par[8]*z/3);
+  r = 1.0 / (1 + z + par[7]*k); /* ~ exp (-Z) */
+  f = par[5]*r + par[6];
+  q = par[5]*r*r*(1 + par[7]*z*(1 + par[8]*z/2));
   /* note difference from gaussian: q = par[5]*r */
-  f = q + par[6];
 
   dpar[0] = q*(2*px*par[2] + par[4]*Y);
@@ -251,6 +316,6 @@
   dpar[5] = +r;
   dpar[6] = +1;
-  dpar[7] = +par[5]*r*z*z/2;
-  dpar[8] = +dpar[7]*z/3;
+  dpar[7] = -100*par[5]*r*r*k;
+  dpar[8] = -100*par[5]*r*r*par[7]*(z*z*z)/6;
 
   return (f);
@@ -261,5 +326,5 @@
 
   float X, Y, px, py;
-  float z, r, q, f;
+  float z, r, q, f, k;
 
   X = x - par[0];
@@ -270,8 +335,9 @@
 
   z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
-  r = 1.0 / (1 + z + par[7]*pow(z,par[8])); /* ~ exp (-Z) */
-  q = par[5]*r*r*par[7]*par[8]*pow(z,par[8]-1.0);
+  k = pow(z,1.75*par[8]);
+  r = 1.0 / (1 + z + par[7]*k); /* ~ exp (-Z) */
+  q = par[5]*r*r*(1 + 1.75*par[7]*par[8]*pow(z,1.75*par[8]-1));
   /* note difference from gaussian: q = par[5]*r */
-  f = q + par[6];
+  f = par[5]*r + par[6];
 
   dpar[0] = q*(2*px*par[2] + par[4]*Y);
@@ -282,4 +348,83 @@
   dpar[5] = +r;
   dpar[6] = +1;
+  dpar[7] = -10*par[5]*r*r*k;
+  dpar[8] = -10*par[5]*r*r*par[7]*k*1.75*log(z);
+
+  return (f);
+}
+
+/* pars: x, y, sx, sy, sxy, I, sky */
+float rgaussTD (float x, float y, float *par, int Npar, float *dpar) {
+
+  float X, Y, px1, py1, px2, py2;
+  float z1, z2, r1, r2, q1, q2, f;
+
+  X = x - par[0];
+  Y = y - par[1];
+  
+  px1 = par[2]*X;
+  py1 = par[3]*Y;
+  px2 = par[8]*X;
+  py2 = par[9]*Y;
+
+  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
+  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[10]*X*Y;
+
+  r1 = 1.0 / (1 + z1 + 0.5*SQ(z1)*(1 + z1/3)); /* ~ exp (-Z) */
+  r2 = 1.0 / (1 + pow(z2,1.75));
+
+  f = par[5]*r1 + par[6] + par[7]*r2;
+
+  q1 = par[5]*SQ(r1)*(1 + z1 + 0.5*SQ(z1));
+  q2 = par[7]*SQ(r2)*(1.75*pow(z2,0.75));
+
+  dpar[	0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[8] + par[10]*Y);
+  dpar[	1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[9] + par[10]*X);
+  dpar[	2] = -2*q1*px1*X;
+  dpar[	3] = -2*q1*py1*Y;
+  dpar[	4] = -q1*X*Y;
+  dpar[	5] = +r1;
+  dpar[	6] = +1;
+  dpar[	7] = +r2*2;
+  dpar[	8] = -2*q2*px2*X*2;
+  dpar[	9] = -2*q2*py2*Y*2;
+  dpar[10] = -q2*X*Y;
+
+  return (f);
+}
+
+/* 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);
Index: /trunk/Ohana/src/opihi/cmd.data/ungridify.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/ungridify.c	(revision 2875)
+++ /trunk/Ohana/src/opihi/cmd.data/ungridify.c	(revision 2876)
@@ -5,5 +5,5 @@
   int i, j, n;
   int Nx, Ny, NX, NY;
-  float Xmin, Xmax, Ymin, Ymax;
+  int Xmin, Xmax, Ymin, Ymax;
   float *v, *x, *y, *z;
   Buffer *bf;
@@ -49,4 +49,7 @@
     }
   }
+  if (n != Nx*Ny) {
+    fprintf (stderr, "error in ungrid: %d vs %d (%d x %d)\n", n, Nx*Ny, Nx, Ny);
+  }
   vx[0].Nelements = vy[0].Nelements = vz[0].Nelements = n;
   return (TRUE);