Changeset 28781 for branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_SERSIC.c

Timestamp:

Jul 29, 2010, 2:35:46 PM (16 years ago)

Author:

eugene

Message:

attempting to get a better starting guess; attempting to downweight neighbor detections

File:

• branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_SERSIC.c (modified) (6 diffs)

branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_SERSIC.c

@@ -126 +126 @@
         dPAR[PM_PAR_SKY]  = +1.0;
         dPAR[PM_PAR_I0]   = +f1;
-        dPAR[PM_PAR_7]    = (z < 0.01) ? -1.5*z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -1.5*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));
@@ -133 +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;
@@ -224 +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;
@@ -236 +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);
 
@@ -242 +262 @@
     if (!isfinite(shape.sxy)) return false;
 
-    // the other parameters depend on the guess for PAR_7
-    if (!isfinite(PAR[PM_PAR_7])) {
-        PAR[PM_PAR_7]    = 0.25;
-        // PAR[PM_PAR_7]    = 0.125;
-        // PAR[PM_PAR_7]    = 0.5;
-    }    
-
-    float index = 0.5 / PAR[PM_PAR_7];
     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);
@@ -329 +345 @@
     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]);