Index: trunk/Ohana/src/opihi/cmd.astro/imfit.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 2886)
+++ 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;
+  }
+  */
+