Changeset 42078 for trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss.c

Timestamp:

Feb 28, 2022, 12:10:28 PM (4 years ago)

Author:

eugene

Message:

add Molweide projection; improvements to source fitting functions (imfit) for testing

File:

• trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss.c (modified) (2 diffs)

trunk/Ohana/src/opihi/cmd.astro/imfit-rgauss.c

@@ -9 +9 @@
 
   fitfunc = rgaussTD;
-  cleanup = rgaussCL;
-  Npar = 10;
-  Nfpar = 1;
+  imfit_cleanup = rgaussCL;
+  Npar =  8;
+  Nfpar = 0;
 
   /* allocate free and fixed parameters */
@@ -21 +21 @@
   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[2] = 2.35 / get_variable_default ("SXg", 2.0);
+  par[3] = 2.35 / get_variable_default ("SYg", 2.0);
+  par[4] = get_variable_default ("SXYg", 0);
   par[5] = get_variable_default ("Zpk", 10000);
-  par[6] = get_variable_default ("Sg", 0.0);
-  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[6] = get_variable_default ("Sg", 0.0); // sky
+  par[7] = get_variable_default ("Sr", 2.0);
 
-  fpar[0] = get_variable_default ("Npow", 2.25);
+//  fpar[0] = get_variable_default ("Npow", 2.25);
 
   sky = &par[6];
 }
 
-/* two components: (1 + z_1 + 0.5*z_1^2 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
+void rgaussCL () {
+  set_variable ("Xg",   par[0]);
+  set_variable ("Yg",   par[1]);
+  set_variable ("SXg",  2.35 / par[2]);
+  set_variable ("SYg",  2.35 / par[3]);
+  set_variable ("SXYg", par[4]);
+  set_variable ("Zpk",  par[5]);
+  set_variable ("Sg",   par[6]);
+  set_variable ("Sr",   par[7]);
+}
+
+/* rgauss: (1 + z + z^alpha)^(-1) -- x, y, sx, sy, sxy, I, sky, sr */
 opihi_flt rgaussTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) {
   OHANA_UNUSED_PARAM(Npar);
 
-  opihi_flt X, Y, px1, py1, px2, py2;
-  opihi_flt z1, z2, r, q1, q2, f;
+  opihi_flt X = x - par[0];
+  opihi_flt Y = y - par[1];
+  
+  opihi_flt px = par[2]*X;
+  opihi_flt py = par[3]*Y;
 
-  X = x - par[0];
-  Y = y - par[1];
-  
-  px1 = par[2]*X;
-  py1 = par[3]*Y;
-  px2 = par[7]*X;
-  py2 = par[8]*Y;
+  opihi_flt z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*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,fpar[0]));
-  f = par[5]*r + par[6];
-
-  q1 = par[5]*SQ(r)*(1 + z1);
-  q2 = par[5]*SQ(r)*fpar[0]*pow(z2,(fpar[0]-1));
+  opihi_flt p = pow(z,par[7] - 1.0);
+  opihi_flt r = 1.0 / (1 + z + z*p);
+  opihi_flt f = par[5]*r + par[6];
 
   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;
+    opihi_flt t = par[5]*SQ(r);
+    opihi_flt q = t*(1 + par[7]*p);
+
+    dpar[0] = +q*(par[2]*px + par[4]*Y);
+    dpar[1] = +q*(par[3]*py + par[4]*X);
+    dpar[2] = -q*px*X;
+    dpar[3] = -q*py*Y;
+    dpar[4] = -q*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;
+
+    // this model derivative is undefined at z = 0.0, but the limit is zero as z -> 0.0
+    dpar[7] = (z == 0.0) ? 0.0 : -t*log(z)*p*z;
   }
   return (f);
 }
 
-int rgaussCL () {
-  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]);
-  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]);
-}
+/******************************************************************************
+ * this file defines the RGAUSS source shape model (XXX need a better name!).  Note that these
+ * model functions are loaded by pmModelClass.c using 'include', and thus need no 'include'
+ * statements of their own.  The models use a psVector to represent the set of parameters, with
+ * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
+ * may thus vary depending on the specifics of the model.  All models which are used as a PSF
+ * representations share a few parameters, for which # define names are listed in pmModel.h:
+
+   power-law with fitted slope
+   1 / (1 + z + z^alpha)
+
+ * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
+ * PM_PAR_I0 1    - central intensity
+ * PM_PAR_XPOS 2  - X center of object
+ * PM_PAR_YPOS 3  - Y center of object
+ * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
+ * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
+ * PM_PAR_SXY 6   - X*Y term of elliptical contour
+ * PM_PAR_7   7   - power-law slope (alpha)
+ *****************************************************************************/
+