Changeset 36680 for branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_EXP.c

Timestamp:

Apr 21, 2014, 5:42:34 AM (12 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ps2-tc3-20130727

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects (modified) (1 prop)
• psModules/src/objects/models/pmModel_EXP.c (modified) (12 diffs)

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects

@@ -12 +12 @@
 pmSourceIO_CMF_PS1_V1.v1.c
 pmSourceIO_CMF_PS1_V4.c
+pmSourceIO_CMF_PS1_V5.c
 pmSourceIO_CMF_PS1_SV1.c
 pmSourceIO_CMF_PS1_SV2.c

@@ -12 +12 @@
 pmSourceIO_CMF_PS1_V1.v1.c
 pmSourceIO_CMF_PS1_V4.c
+pmSourceIO_CMF_PS1_V5.c
 pmSourceIO_CMF_PS1_SV1.c
 pmSourceIO_CMF_PS1_SV2.c

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/models/pmModel_EXP.c

@@ -40 +40 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -45 +46 @@
 #include "pmPSFtry.h"
 #include "pmDetections.h"
+#include "pmModel_CentralPixel.h"
 
 #include "pmModel_EXP.h"
@@ -62 +64 @@
 // 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
 // values need to be pixel coords
+//
+
+// Notes on changing kappa value from 1.70056 to 1.678
+// I'm using a functional form f(x,y) = Io exp(-kappa (r / r_e)).  
+// The article by Graham & Driver (2005) uses a form Ie exp(-bn [(r / r_e) -1]) 
+// which is equal to Ie exp(-bn (r / r_e)) exp(bn).  
+// Thus, my Io = Ie exp(bn) and my kappa is their bn.
+// My value of kappa is 1.700, their value for bn is 1.678., so I am off by a small amount there (1.5%).  
+
+
+#define KAPPA_EXP 1.678
+#define OLD_KAPP_EXP 1.70056
+
 
 // Lax parameter limits
 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.05, 0.05, -1.0 };
-static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0 };
+static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0 };
 
 // Moderate parameter limits
@@ -78 +93 @@
 static float *paramsMinUse = paramsMinLax;
 static float *paramsMaxUse = paramsMaxLax;
-static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5};
+static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5};
 
 static bool limitsApply = true;         // Apply limits?
 
-# include "pmModel_SERSIC.CP.h"
+// # include "pmModel_SERSIC.CP.h"
+
+// the problems I'm having with the SERSIC-like functions are:
+// 1) making sure I have the right functional form so that PAR[SXX,etc] represent R_eff (half-light radius)
+// 2) getting the central pixel right
+// 3) getting the derivaties right.
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -101 +121 @@
     psAssert (z >= 0, "do not allow negative z values in model");
 
-    float index = 1.0;
-    float par7 = 0.5;
-    float bn = 1.9992*index - 0.3271;
-    float Io = exp(bn);
-
-    psF32 f2 = bn*sqrt(z);
-    psF32 f1 = Io*exp(-f2);
-
+    // for EXP, we can hard-wire kappa(1):
+    // float index = 1.0;
+    float kappa = KAPPA_EXP;
+
+    // sqrt(z) is r
+    float q = kappa*sqrt(z);
+    psF32 f0 = exp(-q);
+
+    assert (isfinite(q));
+
+    // only worry about the central 4 pixels at most
     psF32 radius = hypot(X, Y);
-    if (radius < 1.0) {
-
-        // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center
-
-        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
-        psEllipseAxes axes;
-        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
-
-        // get the central pixel flux from the lookup table
-        float xPix = (axes.major - centralPixelXo) / centralPixeldX;
-        xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
-        float yPix = (index - centralPixelYo) / centralPixeldY;
-        yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
-
-        // the integral of a Sersic has an analytical form as follows:
-        float logGamma = lgamma(2.0*index);
-        float bnFactor = pow(bn, 2.0*index);
-        float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
-
-        // XXX interpolate to get the value
-        // XXX for the moment, just integerize
-        // XXX I need to multiply by the integrated flux to get the flux in the central pixel
-        float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
-        
-        float px1 = 1.0 / PAR[PM_PAR_SXX];
-        float py1 = 1.0 / PAR[PM_PAR_SYY];
-        float z10 = PS_SQR(px1);
-        float z01 = PS_SQR(py1);
-
-        // which pixels do we need for this interpolation?
-        // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
-        if ((X >= 0) && (Y >= 0)) {
-            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
-            float V00 = Vcenter;
-            float V10 = Io*exp(-bn*pow(z10,par7));
-            float V01 = Io*exp(-bn*pow(z01,par7));
-            float V11 = Io*exp(-bn*pow(z11,par7));
-            f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
-        }
-        if ((X < 0) && (Y >= 0)) {
-            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
-            float V00 = Io*exp(-bn*pow(z10,par7));
-            float V10 = Vcenter;
-            float V01 = Io*exp(-bn*pow(z11,par7));
-            float V11 = Io*exp(-bn*pow(z01,par7));
-            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
-        }
-        if ((X >= 0) && (Y < 0)) {
-            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
-            float V00 = Io*exp(-bn*pow(z01,par7));
-            float V10 = Io*exp(-bn*pow(z11,par7));
-            float V01 = Vcenter;
-            float V11 = Io*exp(-bn*pow(z10,par7));
-            f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
-        }
-        if ((X < 0) && (Y < 0)) {
-            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
-            float V00 = Io*exp(-bn*pow(z11,par7));
-            float V10 = Io*exp(-bn*pow(z10,par7));
-            float V01 = Io*exp(-bn*pow(z01,par7));
-            float V11 = Vcenter;
-            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
-        }
-    }   
-
-    psF32 z0 = PAR[PM_PAR_I0]*f1;
-    psF32 f0 = PAR[PM_PAR_SKY] + z0;
-
-    assert (isfinite(f2));
+    if (radius <= 1.5) {
+        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], 1.0, 51);
+    }
+    assert (isfinite(f0));
+
+    psF32 f1 = PAR[PM_PAR_I0]*f0;
+    psF32 f = PAR[PM_PAR_SKY] + f1;
+
     assert (isfinite(f1));
-    assert (isfinite(z0));
-    assert (isfinite(f0));
+    assert (isfinite(f));
 
     if (deriv != NULL) {
@@ -187 +148 @@
 
         dPAR[PM_PAR_SKY]  = +1.0;
-        dPAR[PM_PAR_I0]   = +f1;
-
-        // gradient is infinite for z = 0; saturate at z = 0.01
-        // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i]
-        psF32 z1 = (z < 0.01) ? 0.5*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z);
-
-        // XXX dampen SXX and SYY as in GAUSS?
-        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_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
-        dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
-    }
-    return (f0);
+        dPAR[PM_PAR_I0]   = +f0;
+
+        if (z > 0.01) {
+          float z1 = 0.5*f1*kappa/sqrt(z);
+          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]);
+          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]);
+          dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
+          dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
+          dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
+        } else {
+          // gradient -> 0 for z -> 0, but has undef form
+          float z1 = 0.5*f1*kappa;
+          dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]);
+          dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]);
+          dPAR[PM_PAR_SXX]  = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX];
+          dPAR[PM_PAR_SYY]  = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY];
+          dPAR[PM_PAR_SXY]  = -1.0*z1;
+        }
+    }
+    return (f);
 }
 
@@ -284 +251 @@
 bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
 {
+    // for the moment, we are going to require moments and KronFlux
+    if (!source->moments) return false;
+    pmMoments *moments = source->moments;
+
+    if (!isfinite(moments->KronFlux)) return false;
+    if (!isfinite(moments->Mrf)) return false;
+    if (moments->Mrf < 0.0) return false;
+
     psF32 *PAR  = model->params->data.F32;
 
@@ -289 +264 @@
     PAR[PM_PAR_SKY]  = 0.0;
 
-    // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) {
-      return false;
-    }
-
-    // set the model normalization
-    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
-      return false;
-    }
+    psEllipseMoments emoments;
+    emoments.x2 = moments->Mxx;
+    emoments.xy = moments->Mxy;
+    emoments.y2 = moments->Myy;
+
+    // force the axis ratio to be < 20.0
+    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
+
+    if (!isfinite(axes.major)) return false;
+    if (!isfinite(axes.minor)) return false;
+    if (!isfinite(axes.theta)) return false;
+
+    // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width 
+    float scale = moments->Mrf / axes.major;
+    axes.major *= scale;
+    axes.minor *= scale;
+
+    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], axes, true);
+
+    // psEllipseAxes axes;
+    // use the code in SetShape here to avoid doing this 2x
+    // pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
+
+    // float norm = pmSersicNorm (4);  // hardwire
+    float norm = 0.34578;
+    float normFlux = 2.0 * M_PI * axes.major * axes.minor * norm;
+    PAR[PM_PAR_I0] = moments->KronFlux / normFlux;
 
     // set the model position
@@ -306 +299 @@
     return(true);
 }
-
 // An exponential model is equivalent to a Sersic with index = 1.0
 psF64 PM_MODEL_FLUX (const psVector *params)
@@ -314 +306 @@
     psEllipseAxes axes;
     pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
-    float AspectRatio = axes.minor / axes.major;
-
-    float index = 1.0;
-    float bn = 1.9992*index - 0.3271;
-
-    // the integral of a Sersic has an analytical form as follows:
-    float logGamma = lgamma(2.0*index);
-    float bnFactor = pow(bn, 2.0*index);
-    float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
-    
-    psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio;
-
-    return(Flux);
+
+    // static value for EXP:
+    float norm = 0.34578; // \int exp(-kappa*sqrt(z)) r dr
+
+    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
+
+    return(flux);
 }
 
@@ -345 +331 @@
     pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
-    // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f)
-    psF64 zn = log(PAR[PM_PAR_I0] / flux);
-    psF64 radius = axes.major * sqrt (2.0) * zn;
+    // static value for EXP:
+    float kappa = KAPPA_EXP;
+
+    // f = Io exp(-kappa*sqrt(z)) -> sqrt(z) = ln(Io/f) / kappa
+    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
+    psF64 radius = axes.major * zn;
 
     psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
@@ -501 +490 @@
     return;
 }
+
+# if (0)
+void bilin_inter_function () {
+        // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
+        psEllipseAxes axes;
+        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
+
+        // get the central pixel flux from the lookup table
+        float xPix = (axes.major - centralPixelXo) / centralPixeldX;
+        xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1);
+        float yPix = (index - centralPixelYo) / centralPixeldY;
+        yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1);
+
+        // the integral of a Sersic has an analytical form as follows:
+        float logGamma = lgamma(2.0*index);
+        float bnFactor = pow(bn, 2.0*index);
+        float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor;
+
+        // XXX interpolate to get the value
+        // XXX for the moment, just integerize
+        // XXX I need to multiply by the integrated flux to get the flux in the central pixel
+        float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm;
+        
+        float px1 = 1.0 / PAR[PM_PAR_SXX];
+        float py1 = 1.0 / PAR[PM_PAR_SYY];
+        float z10 = PS_SQR(px1);
+        float z01 = PS_SQR(py1);
+
+        // which pixels do we need for this interpolation?
+        // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...)
+        if ((X >= 0) && (Y >= 0)) {
+            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
+            float V00 = Vcenter;
+            float V10 = Io*exp(-bn*pow(z10,par7));
+            float V01 = Io*exp(-bn*pow(z01,par7));
+            float V11 = Io*exp(-bn*pow(z11,par7));
+            f1 = interpolatePixels(V00, V10, V01, V11, X, Y);
+        }
+        if ((X < 0) && (Y >= 0)) {
+            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
+            float V00 = Io*exp(-bn*pow(z10,par7));
+            float V10 = Vcenter;
+            float V01 = Io*exp(-bn*pow(z11,par7));
+            float V11 = Io*exp(-bn*pow(z01,par7));
+            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y);
+        }
+        if ((X >= 0) && (Y < 0)) {
+            float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative
+            float V00 = Io*exp(-bn*pow(z01,par7));
+            float V10 = Io*exp(-bn*pow(z11,par7));
+            float V01 = Vcenter;
+            float V11 = Io*exp(-bn*pow(z10,par7));
+            f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y));
+        }
+        if ((X < 0) && (Y < 0)) {
+            float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive
+            float V00 = Io*exp(-bn*pow(z11,par7));
+            float V10 = Io*exp(-bn*pow(z10,par7));
+            float V01 = Io*exp(-bn*pow(z01,par7));
+            float V11 = Vcenter;
+            f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y));
+        }
+}
+# endif

@@ -12 +12 @@
 pmSourceIO_CMF_PS1_V1.v1.c
 pmSourceIO_CMF_PS1_V4.c
+pmSourceIO_CMF_PS1_V5.c
 pmSourceIO_CMF_PS1_SV1.c
 pmSourceIO_CMF_PS1_SV2.c