Changeset 29060 for branches/sc_branches/trunkTest/psModules/src/objects/models/pmModel_SERSIC.c

Timestamp:

Aug 26, 2010, 9:18:39 AM (16 years ago)

Author:

Serge CHASTEL

Message:

Merging trunk in branch

Location:

branches/sc_branches/trunkTest

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects/models/pmModel_SERSIC.c (modified) (13 diffs)

branches/sc_branches/trunkTest/psModules/src/objects/models/pmModel_SERSIC.c

@@ -18 +18 @@
    * PM_PAR_7   7   - normalized sersic parameter
 
-   note that a standard sersic model uses exp(-K*(z^(1/n) - 1).  the additional elements (K,
+   note that a standard sersic model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
    the -1 offset) are absorbed in this model by the normalization, the exponent, and the
    radial scale.  We fit the elements in this form, then re-normalize them on output.
@@ -25 +25 @@
 #include <stdio.h>
 #include <pslib.h>
-
+#include "pmHDU.h"
+#include "pmFPA.h"
+
+#include "pmTrend2D.h"
+#include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmSpan.h"
+#include "pmFootprintSpans.h"
+#include "pmFootprint.h"
+#include "pmPeaks.h"
 #include "pmMoments.h"
-#include "pmPeaks.h"
+#include "pmModelFuncs.h"
+#include "pmModel.h"
+#include "pmModelUtils.h"
+#include "pmModelClass.h"
+#include "pmSourceMasks.h"
+#include "pmSourceExtendedPars.h"
+#include "pmSourceDiffStats.h"
 #include "pmSource.h"
-#include "pmModel.h"
+#include "pmSourceFitModel.h"
+#include "pmPSF.h"
+#include "pmPSFtry.h"
+#include "pmDetections.h"
+
 #include "pmModel_SERSIC.h"
 
+# define PM_MODEL_NPARAM          8
 # define PM_MODEL_FUNC            pmModelFunc_SERSIC
 # define PM_MODEL_FLUX            pmModelFlux_SERSIC
@@ -47 +67 @@
 
 // Lax parameter limits
-static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0, 0.05 };
+static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.05 };
 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
 
@@ -84 +104 @@
     assert (z >= 0);
 
-    psF32 f2 = pow(z,PAR[PM_PAR_7]);
-    psF32 f1 = exp(-f2);
+    float index = 0.5 / PAR[PM_PAR_7];
+    float bn = 1.9992*index - 0.3271;
+    float Io = exp(bn);
+
+    psF32 f2 = bn*pow(z,PAR[PM_PAR_7]);
+    psF32 f1 = Io*exp(-f2);
     psF32 z0 = PAR[PM_PAR_I0]*f1;
     psF32 f0 = PAR[PM_PAR_SKY] + z0;
@@ -98 +122 @@
 
         // gradient is infinite for z = 0; saturate at z = 0.01
-        psF32 z1 = (z < 0.01) ? z0*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0);
+        psF32 z1 = (z < 0.01) ? z0*bn*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*bn*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0);
 
         dPAR[PM_PAR_SKY]  = +1.0;
         dPAR[PM_PAR_I0]   = +f1;
-        dPAR[PM_PAR_7]    = (z == 0.0) ? 0.0 : -z0*f2*log(z);
+        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z);
+        dPAR[PM_PAR_7]   *= 3.0;
 
         assert (isfinite(z1));
@@ -109 +134 @@
         dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
         dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
-        dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
+        dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag?
         dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
-        dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
     }
@@ -127 +151 @@
         return true;
     }
-    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
+    psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds");
 
     // we need to calculate the limits for SXY specially
@@ -201 +225 @@
     psF32     *PAR  = model->params->data.F32;
 
+    // the other parameters depend on the guess for PAR_7
+    if (!isfinite(PAR[PM_PAR_7])) {
+        PAR[PM_PAR_7]    = 0.25;
+    }    
+    float index = 0.5 / PAR[PM_PAR_7];
+
+    // the scale-length is a function of the moments and the index:
+    // Rmajor_guess = Rmajor_moments * Scale * Zero
+    float Scale = 0.70 + 0.053 * PAR[PM_PAR_7];
+    float Zero  = 1.16 - 0.615 * PAR[PM_PAR_7];
+
     psEllipseMoments emoments;
     emoments.x2 = moments->Mxx;
@@ -213 +248 @@
     if (!isfinite(axes.theta)) return false;
 
+    // set a lower limit to avoid absurd solutions..
+    float Rmajor = PS_MAX(1.0, Scale * axes.major + Zero);
+    float Rminor = Rmajor * (axes.minor / axes.major);
+
+    // fprintf (stderr, "guess index: %f : %f, %f -> %f, %f\n", index, axes.major, axes.minor, Rmajor, Rminor);
+
+    axes.major = Rmajor;
+    axes.minor = Rminor;
     psEllipseShape shape = psEllipseAxesToShape (axes);
 
@@ -219 +262 @@
     if (!isfinite(shape.sxy)) return false;
 
+    float bn = 1.9992*index - 0.3271;
+    // float fR = 1.0 / (sqrt(2.0) * pow (bn, index));
+    float Io = exp(0.5*bn);
+
+    // XXX do we need this factor of sqrt(2)?
+    // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx);
+    // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy);
+
+    float Sxx = PS_MAX(0.5, shape.sx);
+    float Syy = PS_MAX(0.5, shape.sy);
+
     PAR[PM_PAR_SKY]  = 0.0;
-    PAR[PM_PAR_I0]   = peak->flux;
+    PAR[PM_PAR_I0]   = peak->flux / Io;
     PAR[PM_PAR_XPOS] = peak->xf;
     PAR[PM_PAR_YPOS] = peak->yf;
-    PAR[PM_PAR_SXX]  = PS_MAX(0.5, M_SQRT2*shape.sx);
-    PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
+    PAR[PM_PAR_SXX]  = Sxx;
+    PAR[PM_PAR_SYY]  = Syy;
     PAR[PM_PAR_SXY]  = shape.sxy;
-    PAR[PM_PAR_7]    = 0.5;
 
     return(true);
@@ -254 +307 @@
     float f1, f2;
     for (z = DZ; z < 50; z += DZ) {
-        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
+        // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
+        f1 = exp(-pow(z,PAR[PM_PAR_7]));
         z += DZ;
-        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
+        f2 = exp(-pow(z,PAR[PM_PAR_7]));
         norm += f0 + 4*f1 + f2;
         f0 = f2;
@@ -287 +341 @@
 
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 sigma = axes.major;
-
-    psF64 limit = flux / PAR[PM_PAR_I0];
-
-    psF64 z = pow (-log(limit), (1.0 / PAR[PM_PAR_7]));
-    psAssert (isfinite(z), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);
-
-    psF64 radius = sigma * sqrt (2.0 * z);
-    psAssert (isfinite(radius), "fix this code: radius should not be nan for %f, %f", PAR[PM_PAR_7], sigma);
-
-    if (isnan(radius))
-        psAbort("error in code: radius is NaN");
-
+
+    // f = Io exp(-z^n) -> z^n = ln(Io/f)
+    psF64 zn = log(PAR[PM_PAR_I0] / flux);
+    psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]);
+
+    // fprintf (stderr, "sersic model %f %f, n %f, radius: %f, zn: %f, f/Io: %f, major: %f\n", PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], PAR[PM_PAR_7], radius, zn, flux/PAR[PM_PAR_I0], axes.major);
+
+    psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);
     return (radius);
 }
@@ -448 +497 @@
     return;
 }
-
-# undef PM_MODEL_FUNC
-# undef PM_MODEL_FLUX
-# undef PM_MODEL_GUESS
-# undef PM_MODEL_LIMITS
-# undef PM_MODEL_RADIUS
-# undef PM_MODEL_FROM_PSF
-# undef PM_MODEL_PARAMS_FROM_PSF
-# undef PM_MODEL_FIT_STATUS
-# undef PM_MODEL_SET_LIMITS