Index: trunk/Ohana/src/opihi/cmd.astro/imfit.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 2938)
+++ trunk/Ohana/src/opihi/cmd.astro/imfit.c	(revision 3030)
@@ -17,4 +17,7 @@
 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;
@@ -122,4 +125,12 @@
       Npar = 4;
       Nfpar = 6;
+    }
+    if (!strcmp(argv[N], "sersic")) {
+      fitfunc = sersic;
+      Npar = 8;
+    }
+    if (!strcmp(argv[N], "serbulge")) {
+      fitfunc = serbulge;
+      Npar = 12;
     }
     if (fitfunc == NULL) {
@@ -178,5 +189,5 @@
 
   /* make initial guesses on parameters */
-  if (Npar >= 7) {
+  if (Npar >= 7) { /* not Pgauss, Sgauss */
     par[0] = get_variable_default ("Xg", 0);
     par[1] = get_variable_default ("Yg", 0);
@@ -187,5 +198,5 @@
     par[6] = get_variable_default ("Sg", 0.0);
   }
-  if (Npar == 4) {
+  if (Npar == 4) { /* Pgauss, Sgauss */
     par[0] = get_variable_default ("Xg", 0);
     par[1] = get_variable_default ("Yg", 0);
@@ -193,6 +204,20 @@
     par[3] = get_variable_default ("Sg", 0.0);
   }
-  if (Npar == 9) { /** vgauss **/
+  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 **/
@@ -201,18 +226,18 @@
     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;
+    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", 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", 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", 15.0);
+    fpar[5] = get_variable_default ("SXYf", 0.0);
   }    
 
@@ -339,7 +364,18 @@
     set_variable ("Sg",  par[6]);
   }
-  if (Npar == 9) {
+  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) {
@@ -371,4 +407,9 @@
   set_variable ("ChiSq", chisq/(Npts - Npar));
 
+  for (i = 0; i < Npar; i++) {
+    fprintf (stderr, "%g ", par[i]);
+  }
+  fprintf (stderr, "\n");
+
   free (x);
   free (y);
@@ -411,4 +452,84 @@
     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);