Index: trunk/Ohana/src/opihi/cmd.astro/imfit.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 3105)
+++ trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 3144)
@@ -1,55 +1,12 @@
-# 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 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 sersic (float, float, float *, int, float *);
-float serbulge (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);
-}
+# include "imfit.h"
 
 int imfit (int argc, char **argv) {
 
-  int i, j, N, Npts, Npar, Save, VERBOSE, ShapeVariation;
+  int i, j, N, Npts, Save, VERBOSE, ShapeVariation;
   int sx, sy, nx, ny, Nx, Ny;
-  float chisq, ochisq, dchisq, Gain, RDnoise;
-  float *x, *y, *z, *dz, *par, *V, *parmin, *parmax;
+  float chisq, ochisq, dchisq, Gain, RDnoise, SatThreshold;
+  float *x, *y, *z, *dz, *V, *parmin, *parmax;
   char line[64];
   Buffer *buf;
-  float (*fitfunc)(float, float, float *, int, float *);
 
   Save = FALSE;
@@ -63,4 +20,11 @@
     remove_argument (N, &argc, argv);
     ShapeVariation = TRUE;
+  }
+
+  SatThreshold = 0xffff;
+  if ((N = get_argument (argc, argv, "-sat"))) {
+    remove_argument (N, &argc, argv);
+    SatThreshold = atof (argv[N]);
+    remove_argument (N, &argc, argv);
   }
 
@@ -88,50 +52,17 @@
 
   /* set fitting function */
-  Nfpar = 0;
-  Npar = 7;
-  fitfunc = fgaussTD;
+  fgauss_setup ("fgauss");
   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], "vgauss")) {
-      fitfunc = vgaussTD;
-      Npar = 9;
-    }
-    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 (!strcmp(argv[N], "sersic")) {
-      fitfunc = sersic;
-      Npar = 8;
-    }
-    if (!strcmp(argv[N], "serbulge")) {
-      fitfunc = serbulge;
-      Npar = 12;
-    }
+    fgauss_setup (argv[N]);
+    pgauss_setup (argv[N]);
+    Pgauss_setup (argv[N]);
+    sgauss_setup (argv[N]);
+    Sgauss_setup (argv[N]);
+    qgauss_setup (argv[N]);
+    Qgauss_setup (argv[N]);
+    qfgauss_setup (argv[N]);
+    qrgauss_setup (argv[N]);
     if (fitfunc == NULL) {
       fprintf (stderr, "unknown function %s\n", argv[N]);
@@ -155,15 +86,9 @@
   Ny = buf[0].matrix.Naxis[1];
 
-  if (sx < 0) goto range;
-  if (sy < 0) goto range;
-  if (sx + nx >= Nx) goto range;
-  if (sy + ny >= Ny) goto range;
-
-  ALLOCATE (par, float, Npar);
-  bzero (par, Npar*sizeof(float));
-  if (Nfpar > 0) {
-    ALLOCATE (fpar, float, Nfpar);
-    bzero (fpar, Nfpar*sizeof(float));
-  }    
+  /* check if region is valid */
+  if (sx + 0.5*nx < 0) goto range;
+  if (sy + 0.5*ny < 0) goto range;
+  if (sx + 0.5*nx >= Nx) goto range;
+  if (sy + 0.5*ny >= Ny) goto range;
 
   /* convert array z[x,y] to x[i], y[i], z[i] */
@@ -175,8 +100,13 @@
   ALLOCATE (dz, float, 2*Npts);
   for (j = 0; j < ny; j++) {
+    if (j + sy < 0) continue;
+    if (j + sy >= Ny) continue;
     V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx; 
-    for (i = 0; i < nx; i++, V++) {
+    for (i = 0; i < nx; i++) {
+      if (i + sx < 0) continue;
+      if (i + sx >= Nx) continue;
+      if (*V > SatThreshold) goto next;
       dz[N] = (SQ(RDnoise) + *V/Gain);
-      if (dz[N] <= 0) continue;
+      if (dz[N] <= 0) goto next;
       dz[N] = 1.0 / dz[N];
       x[N] = i + sx;
@@ -184,66 +114,13 @@
       z[N] = *V;
       N++;
+    next:
+      V++;
     }
   }
   Npts = N;
 
-  /* make initial guesses on parameters */
-  if (Npar >= 7) { /* not Pgauss, Sgauss */
-    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) { /* Pgauss, Sgauss */
-    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 (fitfunc == vgaussTD) {
-    par[7] = par[8] = 1;
-  }
-  if (fitfunc == sersic) {
-    par[7] = get_variable_default ("Sr", 1.0);
-  }
-  if (fitfunc == serbulge) {
-    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] = get_variable_default ("SXYf", 0.0);
-    par[11] = get_variable_default ("Sr", 1.0);
-    /* arbitrary guess: equal distribution between bulge and disk */
-    par[10] = par[5]/2.0;
-    par[5]  = par[10];
-  }
-  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] = get_variable_default ("SXYf", 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", 0.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", 0.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", 0.0);
-  }    
-
   /* run fit routine */
   ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
   dchisq = ochisq;
-
   for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
     chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
@@ -252,74 +129,9 @@
   }  
   set_int_variable ("Niter",  i);
-  
-
-  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, sky;
+    float *Vi, *Vo, vr, vf;
 
     if ((out = SelectBuffer ("out",   ANYBUFFER, TRUE)) == NULL) return (FALSE);
@@ -335,9 +147,4 @@
     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) ?? */
@@ -350,61 +157,13 @@
 	Vo[(i   )+(j   )*2*nx] = vr;
 	Vo[(i+nx)+(j   )*2*nx] = vf;
-	Vo[(i   )+(j+ny)*2*nx] = vr - vf + sky;
-	Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + sky;
+	Vo[(i   )+(j+ny)*2*nx] = vr - vf + *sky;
+	Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + *sky;
       }
     }
   }
 
-  /* parameter renormalize */
-  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 (fitfunc == vgaussTD) {
-    set_variable ("SXf", par[7]);
-    set_variable ("SYf", par[8]);
-  }
-  if (fitfunc == sersic) {
-    set_variable ("Sr", par[7]);
-  }
-  if (fitfunc == serbulge) {
-    set_variable ("Zb", par[5]);
-    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]);
-    set_variable ("Zd", par[10]);
-    set_variable ("Sr", par[11]);
-  }
-  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]);
-  }    
+  /* save parameters to opihi variables */
+  imfit_cleanup ();
+
   set_variable ("ChiSq", chisq/(Npts - Npar));
 
@@ -421,5 +180,5 @@
   free (dz);
   free (par);
-  if (Nfpar > 0) free (fpar);
+  free (fpar);
 
   mrq2dfree (Npar);
@@ -431,611 +190,2 @@
 }
 
-/* real 2D gaussian -- x, y, sx, sy, sxy, I, sky */
-float fgaussTD (float x, float y, float *par, int Npar, float *dpar) {
-
-  float X, Y, px, py;
-  float z, r, q, f;
-
-  X = x - par[0];
-  Y = y - par[1];
-  
-  px = par[2]*X;
-  py = par[3]*Y;
-
-  z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
-  r = exp (-z);
-  q = par[5]*r;
-  f = q + par[6];
-
-  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);
-}
-
-/* sersic galaxy model -- x, y, sx, sy, sxy, I, sky, n */
-/* exp (-b (r/r_e)^(1/n)) */
-float sersic (float x, float y, float *par, int Npar, float *dpar) {
-
-  float X, Y, px, py;
-  float z, r, t, q, f;
-
-  X = x - par[0];
-  Y = y - par[1];
-  
-  px = par[2]*X;
-  py = par[3]*Y;
-
-  z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
-  t = pow (z, par[7]);
-  r = exp (-t);
-  f = par[5]*r + par[6];
-  q = par[5]*r*par[7]*pow(z, par[7]-1);
-
-  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] = -q*log(z)*t;
-  }
-  return (f);
-}
-
-/*                                  0  1    2   3   4      5    6     7   8   9       10  11 */
-/* sersic galaxy model w/ bulge: -- x, y, (sx, sy, sxy)_1, I_1, sky, (sx, sy, sxy)_2, I_2, n */
-/* exp (-b (r/r_e)^(1/n)) + pgauss (r) */
-float serbulge (float x, float y, float *par, int Npar, float *dpar) {
-
-  float X, Y, px1, px2, py1, py2;
-  float z1, z2, r1, r2, t, 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;
-
-  /* bulge component */
-  r1 = 1.0 / (1 + z1 + 0.5*z1*z1*(1 + z1/3)); /* ~ exp (-Z) */
-
-  /* disk component */
-  t = pow (z2, par[11]);
-  r2 = exp (-t);
-
-  f = par[5]*r1 + par[10]*r2 + par[6];
-
-  q1 = par[5]*r1*r1*(1 + z1 + 0.5*z1*z1);
-  q2 = par[10]*r2*par[11]*pow(z2, par[11]-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] = +r1;
-    dpar[6] = +1;
-    dpar[7] = -2*q2*px2*X*50;
-    dpar[8] = -2*q2*py2*Y*50;
-    dpar[9] = -q2*X*Y*50;
-    dpar[10] = +r2*50;
-    dpar[11] = -q2*log(z2)*t*50;
-  }
-  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;
-
-  X = x - par[0];
-  Y = y - par[1];
-  
-  px = par[2]*X;
-  py = par[3]*Y;
-
-  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) */
-  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 */
-
-  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) {
-
-  float X, Y, px, py;
-  float z, r, q, f, k;
-
-  X = x - par[0];
-  Y = y - par[1];
-  
-  px = par[2]*X;
-  py = par[3]*Y;
-
-  z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
-  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 */
-
-  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) {
-
-  float X, Y, px, py;
-  float z, r, q, f, k;
-
-  X = x - par[0];
-  Y = y - par[1];
-  
-  px = par[2]*X;
-  py = par[3]*Y;
-
-  z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
-  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 = par[5]*r + 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;
-  dpar[7] = -10*par[5]*r*r*k;
-  dpar[8] = -10*par[5]*r*r*par[7]*k*1.75*log(z);
-
-  return (f);
-}
-
-/* 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) {
-
-  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);
-}
-
-# endif
-
-  /* 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;
-  }
-  */
-