Changeset 36021 for branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c

Timestamp:

Aug 25, 2013, 12:08:56 PM (13 years ago)

Author:

eugene

Message:

convert SERSIC to use pmModelCP_SersicSubpix; clean up DEV and EXP versions; add pmPCMMakeModel, pmSersicKappa, pmSersicNorm functions

File:

• branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c (modified) (5 diffs)

branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c

@@ -63 +63 @@
 # define PM_MODEL_SET_LIMITS      pmModelSetLimits_DEV
 
-// f = exp(-z^0.125) 
+// f = exp(-kappa*r^(1/index)) 
+// f = exp(-kappa*z^(0.5/index)) 
+// index = 4, 0.5/index = 0.125
 # define ALPHA 0.125 
-// # define ALPHA 0.25 
 
 // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
@@ -86 +87 @@
 static float *paramsMinUse = paramsMinLax;
 static float *paramsMaxUse = paramsMaxLax;
-static float betaUse[] = { 2, 3e6, 5, 5, 3.0, 3.0, 0.5 };
+static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5 };
 
 static bool limitsApply = true;         // Apply limits?
@@ -113 +114 @@
     // r = sqrt(z)
     float q = kappa*pow(z,ALPHA);
-    psF32 f0 = exp(-q);
-
-    psF32 f1 = PAR[PM_PAR_I0]*f0;
-    psF32 f = PAR[PM_PAR_SKY] + f1;
+    float f0 = exp(-q);
 
     assert (isfinite(q));
     assert (isfinite(f0));
-    assert (isfinite(f1));
-    assert (isfinite(f));
-
-    // only worry about the central 4 pixels at most
-    // If I use DELTA = 0.2, I'm way off for the total flux
-    // If I use DELTA = 0.02, I'm totally good (but I am under on the total flux for R = 30 by 0.2 mags -- aperture failure)
-    // For DELTA = 0.02 & Rmin/Rmaj = 0.25, I'm over flux by 0.15 mags (due to the central pixel)
-    psF32 radius = hypot(X, Y);
+
+    // only worry about the central pixels at most
+    float radius = hypot(X, Y);
     if (radius <= 1.5) {
         // Nsub ~ 10*index^2 + 1
@@ -135 +128 @@
         Nsub = PS_MAX (Nsub, 11);
         f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 4.0, Nsub);
-        f1 = PAR[PM_PAR_I0]*f0;
-        f = PAR[PM_PAR_SKY] + f1;
     }   
+
+    float f1 = PAR[PM_PAR_I0]*f0;
+    float f = PAR[PM_PAR_SKY] + f1;
+
+    assert (isfinite(f1));
+    assert (isfinite(f));
 
     if (deriv != NULL) {
@@ -304 +301 @@
     // f = Io exp(-kappa*z^n) -> z^n = ln(Io/f) / kappa
     psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
-    psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA);
+    psF64 radius = axes.major * pow(zn, 0.5 / ALPHA);
 
     psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)",