Changeset 36680 for branches/eam_branches/ps2-tc3-20130727/psModules/src/objects

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/Makefile.am (modified) (9 diffs)
• psModules/src/objects/mksource.pl (modified) (1 diff)
• psModules/src/objects/models/pmModel_DEV.c (modified) (12 diffs)
• psModules/src/objects/models/pmModel_EXP.c (modified) (12 diffs)
• psModules/src/objects/models/pmModel_GAUSS.c (modified) (2 diffs)
• psModules/src/objects/models/pmModel_PGAUSS.c (modified) (2 diffs)
• psModules/src/objects/models/pmModel_PS1_V1.c (modified) (2 diffs)
• psModules/src/objects/models/pmModel_QGAUSS.c (modified) (2 diffs)
• psModules/src/objects/models/pmModel_RGAUSS.c (modified) (2 diffs)
• psModules/src/objects/models/pmModel_SERSIC.c (modified) (9 diffs)
• psModules/src/objects/models/pmModel_TRAIL.c (modified) (2 diffs)
• psModules/src/objects/pmDetEff.c (modified) (1 diff)
• psModules/src/objects/pmFootprintCullPeaks.c (modified) (4 diffs)
• psModules/src/objects/pmGrowthCurve.c (modified) (1 diff)
• psModules/src/objects/pmGrowthCurveGenerate.c (modified) (1 diff)
• psModules/src/objects/pmModel.c (modified) (1 diff)
• psModules/src/objects/pmModelClass.c (modified) (3 diffs)
• psModules/src/objects/pmModelClass.h (modified) (1 diff)
• psModules/src/objects/pmModelFuncs.h (modified) (1 diff)
• psModules/src/objects/pmModelUtils.c (modified) (4 diffs)
• psModules/src/objects/pmModel_CentralPixel.c (copied) (copied from trunk/psModules/src/objects/pmModel_CentralPixel.c )
• psModules/src/objects/pmModel_CentralPixel.h (copied) (copied from trunk/psModules/src/objects/pmModel_CentralPixel.h )
• psModules/src/objects/pmMoments.h (modified) (1 diff)
• psModules/src/objects/pmPCM_MinimizeChisq.c (modified) (12 diffs)
• psModules/src/objects/pmPCMdata.c (modified) (14 diffs)
• psModules/src/objects/pmPCMdata.h (modified) (3 diffs)
• psModules/src/objects/pmPSF.c (modified) (1 diff)
• psModules/src/objects/pmPSF_IO.c (modified) (1 diff)
• psModules/src/objects/pmPSFtry.c (modified) (1 diff)
• psModules/src/objects/pmPSFtryFitEXT.c (modified) (1 diff)
• psModules/src/objects/pmPSFtryFitPSF.c (modified) (3 diffs)
• psModules/src/objects/pmPSFtryMakePSF.c (modified) (1 diff)
• psModules/src/objects/pmPSFtryMetric.c (modified) (1 diff)
• psModules/src/objects/pmPSFtryModel.c (modified) (3 diffs)
• psModules/src/objects/pmPhotObj.c (modified) (1 diff)
• psModules/src/objects/pmSource.c (modified) (6 diffs)
• psModules/src/objects/pmSource.h (modified) (2 diffs)
• psModules/src/objects/pmSourceContour.c (modified) (1 diff)
• psModules/src/objects/pmSourceExtendedPars.c (modified) (1 diff)
• psModules/src/objects/pmSourceExtendedPars.h (modified) (2 diffs)
• psModules/src/objects/pmSourceFitModel.c (modified) (5 diffs)
• psModules/src/objects/pmSourceFitModel.h (modified) (2 diffs)
• psModules/src/objects/pmSourceFitPCM.c (modified) (5 diffs)
• psModules/src/objects/pmSourceFitSet.c (modified) (3 diffs)
• psModules/src/objects/pmSourceGroups.c (modified) (1 diff)
• psModules/src/objects/pmSourceIO.c (modified) (25 diffs)
• psModules/src/objects/pmSourceIO.h (modified) (3 diffs)
• psModules/src/objects/pmSourceIO_CFF.c (copied) (copied from trunk/psModules/src/objects/pmSourceIO_CFF.c )
• psModules/src/objects/pmSourceIO_CMF.c.in (modified) (19 diffs)
• psModules/src/objects/pmSourceIO_CMP.c (modified) (1 diff)
• psModules/src/objects/pmSourceIO_MatchedRefs.c (modified) (1 diff)
• psModules/src/objects/pmSourceIO_OBJ.c (modified) (1 diff)
• psModules/src/objects/pmSourceIO_PS1_CAL_0.c (modified) (2 diffs)
• psModules/src/objects/pmSourceIO_PS1_DEV_0.c (modified) (2 diffs)
• psModules/src/objects/pmSourceIO_PS1_DEV_1.c (modified) (2 diffs)
• psModules/src/objects/pmSourceIO_RAW.c (modified) (1 diff)
• psModules/src/objects/pmSourceIO_SMPDATA.c (modified) (2 diffs)
• psModules/src/objects/pmSourceIO_SX.c (modified) (1 diff)
• psModules/src/objects/pmSourceLensing.c (copied) (copied from trunk/psModules/src/objects/pmSourceLensing.c )
• psModules/src/objects/pmSourceLensing.h (copied) (copied from trunk/psModules/src/objects/pmSourceLensing.h )
• psModules/src/objects/pmSourceMasks.h (modified) (1 diff)
• psModules/src/objects/pmSourceMatch.c (modified) (1 diff)
• psModules/src/objects/pmSourceMoments.c (modified) (19 diffs)
• psModules/src/objects/pmSourceOutputs.c (modified) (1 diff)
• psModules/src/objects/pmSourcePhotometry.c (modified) (4 diffs)
• psModules/src/objects/pmSourcePlotApResid.c (modified) (1 diff)
• psModules/src/objects/pmSourcePlotMoments.c (modified) (1 diff)
• psModules/src/objects/pmSourcePlotPSFModel.c (modified) (1 diff)
• psModules/src/objects/pmSourceSky.c (modified) (1 diff)
• psModules/src/objects/pmSourceUtils.c (modified) (1 diff)
• psModules/src/objects/pmSourceVisual.c (modified) (1 diff)

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/Makefile.am

@@ -20 +20 @@
         pmModelClass.c \
         pmModelUtils.c \
+        pmModel_CentralPixel.c \
         pmSource.c \
         pmPhotObj.c \
@@ -39 +40 @@
         pmSourceIO_SX.c \
         pmSourceIO_CMP.c \
+        pmSourceIO_CFF.c \
         pmSourceIO_SMPDATA.c \
         pmSourceIO_PS1_DEV_0.c \
@@ -47 +49 @@
         pmSourceIO_CMF_PS1_V3.c \
         pmSourceIO_CMF_PS1_V4.c \
+        pmSourceIO_CMF_PS1_V5.c \
         pmSourceIO_CMF_PS1_SV1.c \
         pmSourceIO_CMF_PS1_SV2.c \
@@ -74 +77 @@
         pmGrowthCurve.c \
         pmSourceMatch.c \
+        pmSourceLensing.c \
         pmDetEff.c \
         pmSourceGroups.c \
@@ -97 +101 @@
         pmModelClass.h \
         pmModelUtils.h \
+        pmModel_CentralPixel.h \
         pmSource.h \
         pmPhotObj.h \
@@ -110 +115 @@
         pmSourceOutputs.h \
         pmSourceIO.h \
-        pmSourceSatstar.h \ 
+        pmSourceSatstar.h \
         pmSourcePlots.h \
         pmSourceVisual.h \
@@ -122 +127 @@
         pmGrowthCurveGenerate.h \
         pmSourceMatch.h \
+        pmSourceLensing.h \
         pmDetEff.h \
         pmSourceGroups.h \
@@ -142 +148 @@
 pmSourceIO_CMF_PS1_V3.c \
 pmSourceIO_CMF_PS1_V4.c \
+pmSourceIO_CMF_PS1_V5.c \
 pmSourceIO_CMF_PS1_DV1.c \
 pmSourceIO_CMF_PS1_DV2.c \
@@ -160 +167 @@
         mksource.pl pmSourceIO_CMF.c.in PS1_V4 pmSourceIO_CMF_PS1_V4.c
 
+pmSourceIO_CMF_PS1_V5.c : pmSourceIO_CMF.c.in mksource.pl
+        mksource.pl pmSourceIO_CMF.c.in PS1_V5 pmSourceIO_CMF_PS1_V5.c
+
 pmSourceIO_CMF_PS1_DV1.c : pmSourceIO_CMF.c.in mksource.pl
         mksource.pl pmSourceIO_CMF.c.in PS1_DV1 pmSourceIO_CMF_PS1_DV1.c

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/mksource.pl

@@ -19 +19 @@
                "PS1_V3", 3,
                "PS1_V4", 4,
+               "PS1_V5", 5,
     );
 %cmfmodes_dv = ("PS1_DV1", 1,

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

@@ -16 +16 @@
    * 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   - normalized dev parameter
 
    note that a standard dev model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
@@ -44 +43 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -49 +49 @@
 #include "pmPSFtry.h"
 #include "pmDetections.h"
+#include "pmModel_CentralPixel.h"
 
 #include "pmModel_DEV.h"
@@ -63 +64 @@
 # 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)
@@ -73 +75 @@
 // Lax parameter limits
 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -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
@@ -86 +88 @@
 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"
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -109 +109 @@
     psAssert (z >= 0, "do not allow negative z values in model");
 
-    float index = 0.5 / ALPHA;
-    float par7 = ALPHA;
-    float bn = 1.9992*index - 0.3271;
-    float Io = exp(bn);
-
-    psF32 f2 = bn*pow(z,ALPHA);
-    psF32 f1 = Io*exp(-f2);
-
-    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));
-        }
+    // for DEV, we can hard-wire kappa(4):
+    // float index = 4.0;
+    float kappa = 7.670628;
+
+    // r = sqrt(z)
+    float q = kappa*pow(z,ALPHA);
+    float f0 = exp(-q);
+
+    assert (isfinite(q));
+
+    // only worry about the central pixels at most
+    float radius = hypot(X, Y);
+    if (radius <= 1.5) {
+        // Nsub ~ 10*index^2 + 1
+        psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_DEV"));
+        int Nsub = 2 * ((int)(25 / axes.minor)) + 1;
+        Nsub = PS_MIN (Nsub, 121);
+        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);
     }   
-
-    psF32 z0 = PAR[PM_PAR_I0]*f1;
-    psF32 f0 = PAR[PM_PAR_SKY] + z0;
-
-    assert (isfinite(f2));
+    assert (isfinite(f0));
+
+    float f1 = PAR[PM_PAR_I0]*f0;
+    float f = PAR[PM_PAR_SKY] + f1;
+
     assert (isfinite(f1));
-    assert (isfinite(z0));
-    assert (isfinite(f0));
+    assert (isfinite(f));
 
     if (deriv != NULL) {
@@ -195 +141 @@
 
         dPAR[PM_PAR_SKY]  = +1.0;
-        dPAR[PM_PAR_I0]   = +2.0*f1; // XXX extra damping..
-
-        // gradient is infinite for z = 0; saturate at z = 0.01
-        psF32 z1 = (z < 0.01) ? z0*bn*ALPHA*pow(0.01,ALPHA - 1.0) : z0*bn*ALPHA*pow(z,ALPHA - 1.0);
-
-        assert (isfinite(z1));
-
-        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 = f1*kappa*ALPHA*pow(z,ALPHA-1.0);
+          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 = f1*kappa*ALPHA*pow(z,ALPHA);
+          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);
 }
 
@@ -292 +244 @@
 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;
 
@@ -297 +257 @@
     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;
-    }
-
-    // the normalization is modified by the slope
-    float index = 0.5 / ALPHA;
-    float bn = 1.9992*index - 0.3271;
-    float Io = exp(0.5*bn);
-
-    // set the model normalization
-    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
-      return false;
-    }
-    PAR[PM_PAR_I0] /= Io;
+    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 
+    // the factor of 2.3 comes from Table 1 of Graham and Driver (2005)
+    float scale = moments->Mrf / axes.major / 2.3;
+    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.00168012;
+    float normFlux = 2.0 * M_PI * axes.major * axes.minor * norm;
+    PAR[PM_PAR_I0] = moments->KronFlux / normFlux;
 
     // set the model position
@@ -328 +301 @@
     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 = 4.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);
+
+    float norm = 0.00168012;
+    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
+
+    return(flux);
 }
 
@@ -359 +324 @@
     pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
-    // 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 / ALPHA);
+    // static value for DEV:
+    float kappa = 7.670628;
+
+    // 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 * 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)", 

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

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

@@ -40 +40 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -61 +62 @@
 // Lax parameter limits
 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -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

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

@@ -40 +40 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -61 +62 @@
 // Lax parameter limits
 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -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

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

@@ -42 +42 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -70 +71 @@
 // Lax parameter limits
 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
-static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
+static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
 
 // Moderate parameter limits
 // Tolerate a small divot (k < 0)
 static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 };
-static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
+static float paramsMaxModerate[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
 
 // Strict parameter limits
 // k = PAR_7 < 0 is very undesirable (big divot in the middle)
 static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
-static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
+static float paramsMaxStrict[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
 
 // Parameter limits to use

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

@@ -42 +42 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -70 +71 @@
 // Lax parameter limits
 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
-static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
+static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
 
 // Moderate parameter limits
 // Tolerate a small divot (k < 0)
 static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 };
-static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
+static float paramsMaxModerate[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
 
 // Strict parameter limits
 // k = PAR_7 < 0 is very undesirable (big divot in the middle)
 static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
-static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
+static float paramsMaxStrict[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 };
 
 // Parameter limits to use

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

@@ -41 +41 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -66 +67 @@
 // Lax parameter limits
 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 };
-static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 };
+static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 4.0 };
 
 // Moderate parameter limits

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

@@ -20 +20 @@
    * note that a Sersic model is usually defined in terms of R_e, the half-light radius.  This
      construction does not include a factor of 2 in the X^2 term, etc, like for a Gaussian.
-     Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using setting:
+     Conversion from SXX, SYY, SXY to R_major, R_minor, theta can be done by using:
      shape.sx = SXX / sqrt(2), shape.sy = SYY / sqrt(2), shape.sxy = SXY, then calling
      psEllipseShapeToAxes, and multiplying the values of axes.major, axes.minor by sqrt(2)
@@ -50 +50 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -55 +56 @@
 #include "pmPSFtry.h"
 #include "pmDetections.h"
+#include "pmModel_CentralPixel.h"
 
 #include "pmModel_SERSIC.h"
@@ -74 +76 @@
 
 // Lax parameter limits
-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 };
+static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.001, 0.001, -1.0, 0.0625 };
+static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 1.0 };
 
 // Moderate parameter limits
@@ -88 +90 @@
 static float *paramsMinUse = paramsMinLax;
 static float *paramsMaxUse = paramsMaxLax;
-static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
+static float betaUse[] = { 2, 3e6, 5, 5, 10.0, 10.0, 0.5, 1.0};
 
 static bool limitsApply = true;         // Apply limits?
 
-# include "pmModel_SERSIC.CP.h"
+// # include "pmModel_SERSIC.CP.h"
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -111 +113 @@
     psAssert (z >= 0, "do not allow negative z values in model");
 
-    float index = 0.5 / PAR[PM_PAR_7];
-    float par7 = PAR[PM_PAR_7];
-    float bn = 1.9992*index - 0.3271;
-    float Io = exp(bn);
-
-    psF32 f2 = bn*pow(z,par7);
-    psF32 f1 = Io*exp(-f2);
+    float Sindex = 0.5 / PAR[PM_PAR_7];
+    float kappa = pmSersicKappa (Sindex);
+
+    float q = kappa*pow(z,PAR[PM_PAR_7]);
+    psF32 f0 = exp(-q);
+
+    assert (isfinite(q));
 
     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;
-
-    if (!isfinite(z0)) {
-        fprintf(stderr, "z0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, f1);
+    if (radius <= 1.5) {
+        // Nsub ~ 10*index^2 + 1
+        psEllipseAxes axes = pmPSF_ModelToAxes(PAR, pmModelClassGetType ("PS_MODEL_SERSIC"));
+        int Nsub = 2 * ((int)(6.0*Sindex / axes.minor)) + 1;
+        Nsub = PS_MIN (Nsub, 121);
+        Nsub = PS_MAX (Nsub, 11);
+        f0 = pmModelCP_SersicSubpix (X, Y, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], Sindex, Nsub);
+    }
+    if (!isfinite(f0)) {
+        fprintf(stderr, "f0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, q);
         fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4],
             PAR[5], PAR[6], PAR[7]);
     }
-
-    assert (isfinite(f2));
+    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) {
@@ -203 +147 @@
 
         dPAR[PM_PAR_SKY]  = +1.0;
-        dPAR[PM_PAR_I0]   = +f1;
-
-        // gradient is infinite for z = 0; saturate at z = 0.01
-        psF32 z1 = (z < 0.01) ? z0*bn*par7*pow(0.01,par7 - 1.0) : z0*bn*par7*pow(z,par7 - 1.0);
-
-        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,par7)*log(0.01) : -z0*f2*log(z);
-        dPAR[PM_PAR_7]   *= 3.0;
-
-        assert (isfinite(z1));
+        dPAR[PM_PAR_I0]   = +f0;
+
+        if (z > 0.01) {
+          float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]-1.0);
+          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;
+          dPAR[PM_PAR_7]    = -1.0*f1*q*log(z);
+        } else {
+          // gradient -> 0 for z -> 0, but has undef form
+          float z1 = f1*kappa*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7]);
+          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;
+          // dPAR[PM_PAR_7]    = -1.0*f1*q*log(z + 0.0001);
+          dPAR[PM_PAR_7]    = -1.0*f1*q*log(z + 0.0001); // factor of 16 to reduce the gain
+        }
         assert (isfinite(dPAR[PM_PAR_7]));
-
-        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]; // 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;
-    }
-    return (f0);
+    }
+    return (f);
 }
 
@@ -370 +320 @@
     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 = 0.5 / PAR[PM_PAR_7];
-    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);
+
+    float Sindex = 0.5 / PAR[PM_PAR_7];
+    float norm = pmSersicNorm (Sindex);
+
+    float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm;
+
+    return(flux);
 }
 
@@ -401 +345 @@
     pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
+    float Sindex = 0.5 / PAR[PM_PAR_7];
+    float kappa = pmSersicKappa (Sindex);
+
     // 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]);
+    psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa;
+    psF64 radius = axes.major * pow(zn, Sindex);
 
     psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f", 

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

@@ -40 +40 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -61 +62 @@
 // Lax parameter limits 
 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -1.0e2, -1.0e2,   0.5, -3.3, -0.5 };
-static float paramsMaxLax[] = {  1.0e5, 1.0e+8, +1.0e4, +1.0e4, 150.0, +3.3 , 5.0 };
+static float paramsMaxLax[] = {  1.0e5, 1.00+9, +1.0e5, +1.0e5, 150.0, +3.3 , 5.0 };
 
 // Moderate parameter limits

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

@@ -24 +24 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -25 +25 @@
 bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp);
 
- /*
-  * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
-  * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
-  * to reach a still higher peak --- and if that coll's more (less?) than nsigma DN below your
-  * starting point, discard the peak.
-  */
+/*
+ * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
+ * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
+ * to reach a still higher peak --- and if that coll's more (less?) than nsigma DN below your
+ * starting point, discard the peak.
+ */
 
 # define IN_PEAK 1
@@ -48 +48 @@
 
     if (fp->peaks == NULL || fp->peaks->n < 2) { // nothing to do
-        return PS_ERR_NONE;
+        return PS_ERR_NONE;
     }
 
@@ -91 +91 @@
 
         // max flux is above threshold for brightest peak
-      pmPeak *maxPeak = NULL;
-      for (int i = 0; i < fp->peaks->n; i++) {
-        pmPeak *testPeak = fp->peaks->data[i];
-        float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux;
+        pmPeak *maxPeak = NULL;
+        for (int i = 0; i < fp->peaks->n; i++) {
+            pmPeak *testPeak = fp->peaks->data[i];
+            float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux;
         
-        if (isfinite(this_peak)) {
-          maxPeak = fp->peaks->data[i];
-          break;
-        }
-      }
-      psAssert(maxPeak,"maxPeak was not set in these peaks");
-      //      = fp->peaks->data[0];
+            if (isfinite(this_peak)) {
+                maxPeak = fp->peaks->data[i];
+                break;
+            }
+        }
+        psAssert(maxPeak,"maxPeak was not set in these peaks");
+        //      = fp->peaks->data[0];
         float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux;
 
@@ -130 +130 @@
         }
 #if (0)
-        if (threshbounds->data.F32[threshbounds->n-1] > maxFlux) {
-            psWarning ("upper limit: %f does not include max flux: %f",
-                    threshbounds->data.F32[threshbounds->n-1], maxFlux);
-        }
+        if (threshbounds->data.F32[threshbounds->n-1] > maxFlux) {
+            psWarning ("upper limit: %f does not include max flux: %f",
+                       threshbounds->data.F32[threshbounds->n-1], maxFlux);
+        }
 #endif
         psHistogram *threshist = psHistogramAllocGeneric(threshbounds);

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

@@ -37 +37 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -217 +217 @@
     // the options allow us to modify various aspects of the model
     if (mode & PM_MODEL_OP_NORM) {
+        // if we are including the sky, renormalizing should force use to normalized down the sky flux
+        params->data.F32[PM_PAR_SKY] /= params->data.F32[PM_PAR_I0];
         params->data.F32[PM_PAR_I0] = 1.0;
     }
     if (!(mode & PM_MODEL_OP_SKY)) {
         params->data.F32[PM_PAR_SKY] = 0.0;
-    }
+    } 
     if (mode & PM_MODEL_OP_CENTER) {
         params->data.F32[PM_PAR_XPOS] = image->col0 + 0.5*image->numCols;

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

@@ -66 +66 @@
 
 static pmModelClass *models = NULL;
+static psVector *modelClassLookupTable = NULL;  // translation between model types in header and here
 static int Nmodels = 0;
 
@@ -135 +136 @@
     models = NULL;
     Nmodels = 0;
+    psFree(modelClassLookupTable);
+    modelClassLookupTable = NULL;
     return;
 }
@@ -193 +196 @@
 }
 
+
+bool pmModelClassWriteHeader(psMetadata *header)
+{
+    psMetadataAddS32(header, PS_LIST_TAIL, "MTNUM", PS_META_REPLACE, "number of model types", Nmodels);
+    for (int i = 0; i < Nmodels; i++) {
+        char modelNameKey[16];
+        char modelValKey[16];
+        sprintf(modelNameKey, "MTNAM%02d", i);
+        sprintf(modelValKey,  "MTVAL%02d", i);
+        psMetadataAddStr(header, PS_LIST_TAIL, modelNameKey, PS_META_REPLACE, "", models[i].name);
+        psMetadataAddS32(header, PS_LIST_TAIL, modelValKey, PS_META_REPLACE, "", i);
+    }
+
+    return true;
+}
+
+bool pmModelClassReadHeader(psMetadata *header) {
+    psFree(modelClassLookupTable);
+
+    bool status;
+    int numHeaderModels = psMetadataLookupS32(&status, header, "MTNUM");
+    if (!status) {
+        return false;
+    }
+
+    psVector *inputTypes = psVectorAlloc(numHeaderModels, PS_TYPE_S32);
+    psVector *localTypes = psVectorAlloc(numHeaderModels, PS_TYPE_S32);
+    int max_val = -1;
+    for (int i = 0; i < numHeaderModels; i++) {
+        char modelNameKey[16];
+        char modelValKey[16];
+        sprintf(modelNameKey, "MTNAM%02d", i);
+        sprintf(modelValKey,  "MTVAL%02d", i);
+        psString thisName = psMetadataLookupStr(&status, header, modelNameKey);
+        int thisVal = psMetadataLookupS32(&status, header, modelValKey);
+        if (thisVal > max_val) {
+            max_val = thisVal;
+        }
+        inputTypes->data.S32[i] = thisVal;
+        localTypes->data.S32[i] = pmModelClassGetType(thisName);
+    }
+    if (max_val < 0) {
+        psFree(inputTypes);
+        psFree(localTypes);
+        return false;
+    }
+
+    modelClassLookupTable = psVectorAlloc(max_val + 1, PS_TYPE_S32);
+    psVectorInit(modelClassLookupTable, -1);
+
+    for (int i = 0; i < numHeaderModels; i++) {
+        int thisVal = inputTypes->data.S32[i];
+        int localVal = localTypes->data.S32[i];
+        modelClassLookupTable->data.S32[thisVal] = localVal;
+    }
+    psFree(inputTypes);
+    psFree(localTypes);
+
+    return true;
+}
+
+pmModelType pmModelClassGetLocalType(pmModelType inputType) {
+    pmModelType localType = -1;
+
+    if (modelClassLookupTable) {
+        if (inputType >= 0 && inputType < modelClassLookupTable->n) {
+            localType = modelClassLookupTable->data.S32[inputType];
+        }
+    } else {
+        // no lookup table defined
+        // for backwards compatability if inputType refers to a defined model, return it
+        if (inputType >= 0 && pmModelClassGetName(inputType)) {
+            localType = inputType;
+        }
+    }
+
+    return localType;
+}

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelClass.h

@@ -76 +76 @@
 void pmModelClassSetLimits(pmModelLimitsType type);
 
+// write keywords to header definining the model type values used by this program
+bool pmModelClassWriteHeader(psMetadata *header);
+// create a lookup table for translating input model type values to local model type values
+bool pmModelClassReadHeader(psMetadata *header);
+// translate input model type value to local value
+pmModelType pmModelClassGetLocalType(pmModelType inputType);
 
 /// @}

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmModelFuncs.h

@@ -36 +36 @@
 
 typedef enum {
-    PM_MODEL_STATUS_NONE         = 0x00, ///< model fit not yet attempted, no other info
-    PM_MODEL_STATUS_FITTED       = 0x01, ///< model fit completed
-    PM_MODEL_STATUS_NONCONVERGE  = 0x02, ///< model fit did not converge
-    PM_MODEL_STATUS_OFFIMAGE     = 0x04, ///< model fit drove out of range
-    PM_MODEL_STATUS_BADARGS      = 0x08, ///< model fit called with invalid args
-    PM_MODEL_STATUS_LIMITS       = 0x10, ///< model parameters hit limits
-    PM_MODEL_STATUS_WEAK_FIT     = 0x20, ///< model fit met loose tolerance, but not tight tolerance
+    PM_MODEL_STATUS_NONE           = 0x000, ///< model fit not yet attempted, no other info
+    PM_MODEL_STATUS_FITTED         = 0x001, ///< model fit completed
+    PM_MODEL_STATUS_NONCONVERGE    = 0x002, ///< model fit did not converge
+    PM_MODEL_STATUS_OFFIMAGE       = 0x004, ///< model fit drove out of range
+    PM_MODEL_STATUS_BADARGS        = 0x008, ///< model fit called with invalid args
+    PM_MODEL_STATUS_LIMITS         = 0x010, ///< model parameters hit limits
+    PM_MODEL_STATUS_WEAK_FIT       = 0x020, ///< model fit met loose tolerance, but not tight tolerance
+    PM_MODEL_STATUS_NAN_CHISQ      = 0x040, ///< model fit failed with a NAN chisq 
+    PM_MODEL_SERSIC_PCM_FAIL_GUESS = 0x080, ///< sersic model fit failed on the initial moments-based guess
+    PM_MODEL_SERSIC_PCM_FAIL_GRID  = 0x100, ///< sersic model fit failed on the grid search
+    PM_MODEL_PCM_FAIL_GUESS        = 0x200, ///< non-sersic model fit failed on the initial moments-based guess
+    PM_MODEL_BEST_FIT              = 0x400, ///< this model was the best fit and was subtracted
 } pmModelStatus;
 

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

@@ -39 +39 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -129 +130 @@
 bool pmModelAxesToParams (float *Sxx, float *Sxy, float *Syy, psEllipseAxes axes, bool useReff)  {
 
+    // restrict axex to 0.5 here not below 
+    if (axes.minor < 0.2) axes.minor = 0.2;
+    if (axes.major < 0.2) axes.major = 0.2;
+
     psEllipseShape shape = psEllipseAxesToShape (axes);
 
@@ -137 +142 @@
     // set the shape parameters
     if (useReff) {
-        *Sxx  = PS_MAX(0.5, shape.sx);
-        *Syy  = PS_MAX(0.5, shape.sy);
+        // *Sxx  = PS_MAX(0.5, shape.sx);
+        // *Syy  = PS_MAX(0.5, shape.sy);
+        *Sxx  = shape.sx;
+        *Syy  = shape.sy;
         *Sxy  = shape.sxy * 2.0;
     } else {
-        *Sxx  = PS_MAX(0.5, M_SQRT2*shape.sx);
-        *Syy  = PS_MAX(0.5, M_SQRT2*shape.sy);
+        // *Sxx  = PS_MAX(0.5, M_SQRT2*shape.sx);
+        // *Syy  = PS_MAX(0.5, M_SQRT2*shape.sy);
+        *Sxx  = M_SQRT2*shape.sx;
+        *Syy  = M_SQRT2*shape.sy;
         *Sxy  = shape.sxy;
     }
@@ -190 +199 @@
 
     if (!isfinite(axes.major)) return false;
+    if (axes.major == 0) return false;
     if (!isfinite(axes.minor)) return false;
     if (!isfinite(axes.theta)) return false;
+    if (axes.major == 0) return false;
 
     // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width 

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmMoments.h

@@ -44 +44 @@
     float Myyyy;   ///< fourth moment
 
+  // float wSum;    ///< window-weighted sum (NOT needed by lensing)
+
     float Sum;    ///< Pixel sum above sky (background).
     float Peak;   ///< Peak counts above sky.

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

@@ -38 +38 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
 #include "pmPCMdata.h"
+
+# define SAVE_IMAGES 0
+# if (SAVE_IMAGES) 
+int psphotSaveImage (psMetadata *header, psImage *image, char *filename);
+# endif
 
 # define FACILITY "psModules.objects"
@@ -91 +97 @@
     psF32 lambda = 0.001;
     psF32 dLinear = 0.0;
-    psF32 nu = 2.0;
+    psF32 nu = 3.0;
 
 # if (USE_FFT && PRE_CONVOLVE)
@@ -130 +136 @@
         }
 
+        if (min->isInteractive) {
+            fprintf (stderr, "%d : ", min->iter);
+            for (int ti = 0; ti < params->n; ti++) {
+                fprintf (stderr, "%f  ", params->data.F32[ti]);
+            }
+            fprintf (stderr, " : %f\n", min->value);
+        }
+
+        char key[10]; // used for interactive responses
+        bool testValue = false;
+
         // set a new guess for Alpha, Beta, Params
         if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) {
+            if (false && min->isInteractive) {
+                fprintf (stdout, "guess failed (singular matrix or NaN values), continue? [Y,n] ");
+                if (!fgets(key, 8, stdin)) {
+                    psWarning("Unable to read option");
+                }
+                switch (key[0]) {
+                  case 'n':
+                  case 'N':
+                    done = true;
+                    break;
+                  case 'y':
+                  case 'Y':
+                  case '\n':
+                    lambda *= 10.0;
+                    continue;
+                  default:
+                    lambda *= 10.0;
+                    continue;
+                }
+                if (done) break;
+            }
             min->iter ++;
             if (min->iter >=  min->maxIter) break;
@@ -138 +176 @@
         }
 
+        if (false && min->isInteractive) {
+            p_psVectorPrint(psTraceGetDestination(), Params, "current parameters: ");
+            fprintf (stdout, "last chisq : %f\n", min->value);
+            bool getOptions = true;
+            while (getOptions) {
+                fprintf (stdout, "options: (m)odify, (g)o, (q)uit: ");
+                if (!fgets(key, 8, stdin)) {
+                    psWarning("Unable to read option");
+                }
+                switch (key[0]) {
+                  case 'm':
+                  case 'M':
+                    testValue = TRUE;
+                    fprintf (stdout, "enter (Npar) (value): ");
+                    int Npar = 0;
+                    float value= 0;
+                    int Nscan = fscanf (stdin, "%d %f", &Npar, &value);
+                    if (Nscan != 2) {
+                      fprintf (stderr, "scan failure\n");
+                    }
+                    Params->data.F32[Npar] = value;
+                    break;
+                  case 'g':
+                  case 'G':
+                  case '\n':
+                    getOptions = false;
+                    break;
+                  default:
+                    done = true;
+                    break;
+                }
+                fprintf (stderr, "foo\n");
+            }
+            if (done) break;
+        }
+            
         // dump some useful info if trace is defined
         if (psTraceGetLevel(FACILITY) >= 6) {
@@ -202 +276 @@
         // XXX : Madsen gives suggestion for better use of rho
         // rho is positive if the new chisq is smaller
-        if (rho >= -1e-6) {
+        if (testValue || (rho >= -1e-6)) {
             min->value = Chisq;
             alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
@@ -215 +289 @@
           case 0:
             if (rho >= -1e-6) {
-                lambda *= 0.25;
+                lambda *= 0.1;
             } else {
                 lambda *= 10.0;
@@ -234 +308 @@
             if (rho > 0.0) {
                 lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0)));
-                nu = 2.0;
+                nu = 3.0;
             } else {
                 lambda *= nu;
-                nu *= 2.0;
+                nu *= 3.0;
             }
             break;
@@ -408 +482 @@
 # else
     if (pcm->use1Dgauss) {
-        // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
-        // * the model flux is not masked
-        // * threading takes place above this level
-        pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
-        psImageSmooth_PreAlloc_F32 (pcm->modelConvFlux, pcm->smdata);
-        // psImageSmooth (pcm->modelConvFlux, pcm->sigma, pcm->nsigma);
+
+        if (USE_1D_CACHE) {
+            // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
+            // * the model flux is not masked
+            // * threading takes place above this level
+            pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
+            psImageSmoothCache_F32 (pcm->modelConvFlux, pcm->smdata);
+        } else {
+            pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
+            psImageSmooth2dCache_F32 (pcm->modelConvFlux, pcm->smdata2d);
+        }
     } else {
         psImageConvolveKernel (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psfFFT);
@@ -428 +507 @@
 # else
         if (pcm->use1Dgauss) {
-            // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
-            // * the model flux is not masked
-            // * threading takes place above this level
-            dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
-            psImageSmooth_PreAlloc_F32 (dmodelConv, pcm->smdata);
-            // psImageSmooth (dmodelConv, pcm->sigma, pcm->nsigma);
+            if (USE_1D_CACHE) {
+                // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
+                // * the model flux is not masked
+                // * threading takes place above this level
+                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
+                psImageSmoothCache_F32 (dmodelConv, pcm->smdata);
+            } else {
+                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
+                psImageSmooth2dCache_F32 (dmodelConv, pcm->smdata2d);
+            }
         } else {
             psImageConvolveKernel (dmodelConv, dmodel, NULL, 0, pcm->psfFFT);
@@ -449 +532 @@
 
         if (pcm->use1Dgauss) {
-            // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
-            // * the model flux is not masked
-            // * threading takes place above this level
-            dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
-            psImageSmooth_PreAlloc_F32 (dmodelConv, pcm->smdata);
-            // psImageSmooth (dmodelConv, pcm->sigma, pcm->nsigma);
+            if (USE_1D_CACHE) {
+                // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
+                // * the model flux is not masked
+                // * threading takes place above this level
+                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
+                psImageSmoothCache_F32 (dmodelConv, pcm->smdata);
+            } else {
+                dmodelConv = psImageCopy (dmodelConv, dmodel, dmodel->type.type);
+                psImageSmooth2dCache_F32 (dmodelConv, pcm->smdata2d);
+            }
         } else {
             psImageConvolveFFT (dmodelConv, dmodel, NULL, 0, pcm->psf);
@@ -474 +561 @@
     // XXX TEST : SAVE IMAGES
 # if (SAVE_IMAGES)
-    psphotSaveImage (NULL, pcm->psf->image, "psf.fits");
-    psphotSaveImage (NULL, pcm->modelFlux, "model.fits");
-    psphotSaveImage (NULL, pcm->modelConvFlux, "modelConv.fits");
-    psphotSaveImage (NULL, source->pixels, "obj.fits");
-    psphotSaveImage (NULL, source->maskObj, "mask.fits");
-    psphotSaveImage (NULL, source->variance, "variance.fits");
+    static int Npass = 0;
+    char name[128]; 
+    if (!pcm->use1Dgauss) {
+      snprintf (name, 128, "psf.%03d.fits", Npass); psphotSaveImage (NULL, pcm->psf->image, name);
+    }
+    snprintf (name, 128, "mod.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelFlux, name);
+    snprintf (name, 128, "cnv.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelConvFlux, name);
+    snprintf (name, 128, "obj.%03d.fits", Npass); psphotSaveImage (NULL, source->pixels, name);
+    snprintf (name, 128, "msk.%03d.fits", Npass); psphotSaveImage (NULL, source->maskObj, name);
+    snprintf (name, 128, "var.%03d.fits", Npass); psphotSaveImage (NULL, source->variance, name);
+    for (int n = 0; n < pcm->dmodelsFlux->n; n++) {
+        psImage *dmodelConv = pcm->dmodelsConvFlux->data[n];
+        if (!dmodelConv) continue;
+        snprintf (name, 128, "dpar.%01d.%03d.fits", n, Npass); psphotSaveImage (NULL, dmodelConv, name);
+    }
+    Npass ++;
 # endif
 
@@ -506 +603 @@
 
             float ymodel  = pcm->modelConvFlux->data.F32[i][j];
-            float yweight = 1.0 / source->variance->data.F32[i][j];
+
+            // XXXX note this point here:::
+            float yweight = pcm->poissonErrors ? 1.0 / source->variance->data.F32[i][j] : 1.0;
             float delta = ymodel - source->pixels->data.F32[i][j];
 

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

@@ -38 +38 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -43 +44 @@
 
 # define USE_DELTA_PSF 0
-# define USE_1D_GAUSS 1
 
 static void pmPCMdataFree (pmPCMdata *pcm) {
@@ -58 +58 @@
     psFree (pcm->psfFFT);
     psFree (pcm->constraint);
+
     psFree (pcm->smdata); // pre-allocated data for psImageSmooth_PreAlloc
+    psFree (pcm->smdata2d); // pre-allocated data for psImageSmooth_PreAlloc
     return;
 }
@@ -88 +90 @@
     }
 
+    pcm->smdata = NULL;
+    pcm->smdata2d = NULL;
+
     pcm->modelConv = NULL;
     pcm->psf = NULL;
@@ -94 +99 @@
     pcm->nDOF = 0;
 
+    pcm->poissonErrors = true;
+
     // full convolution with the PSF is expensive.  if we have to save time, we can do a 1D
     // convolution with a Gaussian approximation to the kernel
     pcm->use1Dgauss = false;
-    pcm->nsigma = 3.0; 
+    pcm->nsigma = NAN; // this is set to something defined by the user
     pcm->sigma = 1.0; // this should be set to something sensible when the psf is known
 
@@ -173 +180 @@
 }
 
+int pmPCMsetParams (psMinConstraint *constraint, pmSourceFitMode mode) {
+
+    // set parameter mask based on fitting mode
+    int nParams = 0;
+    int nParAll = constraint->paramMask->n;
+
+    switch (mode) {
+      case PM_SOURCE_FIT_NORM:
+        // fits only source normalization (Io)
+        nParams = 1;
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        break;
+
+      case PM_SOURCE_FIT_PSF:
+        // fits only x,y,Io
+        nParams = 3;
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
+        break;
+
+      case PM_SOURCE_FIT_EXT:
+        // fits all params except sky
+        nParams = nParAll - 1;
+        psVectorInit (constraint->paramMask, 0);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
+        break;
+
+      case PM_SOURCE_FIT_EXT_AND_SKY:
+        // fits all params including sky
+        nParams = nParAll;
+        psVectorInit (constraint->paramMask, 0);
+        break;
+
+      case PM_SOURCE_FIT_SHAPE:
+        // fits shape (Sxx, Sxy, Syy) and Io
+        nParams = 5;
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SXX] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SXY] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SYY] = 0;
+        break;
+
+      case PM_SOURCE_FIT_INDEX:
+        // fits only Io, index (PAR7) -- only Io for models with < 8 params
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        if (nParAll == 7) {
+            nParams = 1;
+        } else {
+            nParams = 2;
+            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
+        }
+        break;
+
+      case PM_SOURCE_FIT_NO_INDEX:
+        // fits all but index (PAR7) including sky
+        psVectorInit (constraint->paramMask, 0);
+        if (nParAll == 7) {
+            nParams = nParAll;
+        } else {
+            nParams = nParAll - 1;
+            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
+        }
+        break;
+      default:
+        psAbort("invalid fitting mode");
+    }
+    return nParams;
+}
+
+static int modelType_GAUSS = -1;
+static int modelType_PS1_V1 = -1;
+
+// generate a Gaussian smoothing kernel for supplied sigma.  sigma here does not need to match
+// that used to allocate the structure, but it is recommended
+bool psImageSmoothCacheKernel_PS1_V1 (psImageSmoothCacheData *smdata, float sigma, float kappa) {
+    // check for NULL structure elements?
+
+    int size = smdata->Nrange;
+
+    psFree (smdata->kernel);
+    smdata->kernel = psVectorAlloc(2 * smdata->Nrange + 1, PS_TYPE_F32);
+
+    double sum = 0.0;                   // Sum of Gaussian, for normalization
+    double factor = 1.0 / (sigma * M_SQRT2);    // Multiplier for i -> z
+
+    // PS1_V1 is a power-law with fitted linear term:
+    // 1 / (1 + kappa z + z^1.666)  where z = (r/sigma)^2
+
+    // generate the kernel (not normalized)
+    for (int i = -size, j = 0; i <= size; i++, j++) {
+        float z = PS_SQR(i * factor);
+        sum += smdata->kernel->data.F32[j] = 1.0 / (1 + kappa * z + pow(z,1.666));
+    }
+
+    // renormalize kernel to integral of 1.0
+    for (int i = 0; i < 2 * size + 1; i++) {
+        smdata->kernel->data.F32[i] /= sum;
+    }
+
+    return true;
+}
+
+psImageSmoothCacheData *psImageSmoothCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *modelPSF) {
+
+    psAssert (modelPSF, "psf model must be defined");
+    
+    psEllipseAxes axes;
+    bool useReff = pmModelUseReff (modelPSF->type);
+    psF32 *PAR = modelPSF->params->data.F32;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
+    
+    *sigma = NAN;
+    *kappa = NAN;
+
+    // XXX need to do this more carefully
+    if (modelPSF->type == modelType_GAUSS) {
+        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
+        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+        *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+        *sigma = 0.5 * (axes.major + axes.minor);
+        *kappa = PAR[PM_PAR_7];
+    }
+    psAssert (isfinite(*sigma), "invalid model type");
+
+    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
+    psImageSmoothCacheData *smdata = psImageSmoothCacheAlloc (flux, *sigma, nsigma);
+
+    if (modelPSF->type == modelType_GAUSS) {
+        psImageSmoothCacheKernel_Gauss (smdata, *sigma);
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+        psImageSmoothCacheKernel_PS1_V1 (smdata, *sigma, *kappa);
+    }
+
+    return smdata;
+}
+
+psImageSmooth2dCacheData *psImageSmooth2dCacheSetKernel (float *sigma, float *kappa, float nsigma, psImage *flux, pmModel *modelPSF) {
+
+    psAssert (modelPSF, "psf model must be defined");
+    
+    psEllipseAxes axes;
+    bool useReff = pmModelUseReff (modelPSF->type);
+    psF32 *PAR = modelPSF->params->data.F32;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
+    
+    *sigma = NAN;
+    *kappa = NAN;
+
+    // XXX need to do this more carefully
+    if (modelPSF->type == modelType_GAUSS) {
+        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
+        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+        *sigma = 0.50 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+        *sigma = 0.5 * (axes.major + axes.minor);
+        *kappa = PAR[PM_PAR_7];
+    }
+    psAssert (isfinite(*sigma), "invalid model type");
+
+    // psImageSmoothCacheAlloc generates a structure but does not assign the smoothing vector
+    psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (nsigma);
+
+    if (modelPSF->type == modelType_GAUSS) {
+        psImageSmooth2dCacheKernel_Gauss (smdata, *sigma);
+    }
+    if (modelPSF->type == modelType_PS1_V1) {
+        psImageSmooth2dCacheKernel_PS1_V1 (smdata, *sigma, *kappa);
+    }
+
+    return smdata;
+}
+
 pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize) {
 
-    // make sure we save a cached copy of the psf flux
-    pmSourceCachePSF (source, maskVal);
-
-    // convert the cached cached psf model for this source to a psKernel
-    psKernel *psf = pmPCMkernelFromPSF (source, psfSize);
-    if (!psf) {
-        // NOTE: this only happens if the source is too close to an edge
-        model->flags |= PM_MODEL_STATUS_BADARGS;
-        return NULL;
-    }
-
-# if (USE_DELTA_PSF)
-    psImageInit (psf->image, 0.0);
-    psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0;
-# endif
+    modelType_GAUSS = pmModelClassGetType ("PS_MODEL_GAUSS");
+    modelType_PS1_V1 = pmModelClassGetType ("PS_MODEL_PS1_V1");
 
     // count the number of unmasked pixels:
@@ -219 +395 @@
     constraint->checkLimits = model->modelLimits;
 
-    // set parameter mask based on fitting mode
-    int nParams = 0;
-    switch (fitOptions->mode) {
-      case PM_SOURCE_FIT_NORM:
-        // NORM-only model fits only source normalization (Io)
-        nParams = 1;
-        psVectorInit (constraint->paramMask, 1);
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        break;
-      case PM_SOURCE_FIT_PSF:
-        // PSF model only fits x,y,Io
-        nParams = 3;
-        psVectorInit (constraint->paramMask, 1);
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
-        break;
-      case PM_SOURCE_FIT_EXT:
-        // EXT model fits all params (except sky)
-        nParams = params->n - 1;
-        psVectorInit (constraint->paramMask, 0);
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
-        break;
-      case PM_SOURCE_FIT_INDEX:
-        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
-        psVectorInit (constraint->paramMask, 1);
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        if (params->n == 7) {
-            nParams = 1;
-        } else {
-            nParams = 2;
-            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
-        }
-        break;
-      case PM_SOURCE_FIT_NO_INDEX:
-        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
-        psVectorInit (constraint->paramMask, 0);
-        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
-        if (params->n == 7) {
-            nParams = params->n - 1;
-        } else {
-            nParams = params->n - 2;
-            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
-        }
-        break;
-      default:
-        psAbort("invalid fitting mode");
-    }
+    int nParams = pmPCMsetParams (constraint, fitOptions->mode);
 
     if (nPix <  nParams + 1) {
         psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
-        psFree (psf);
         psFree (constraint);
         model->flags |= PM_MODEL_STATUS_BADARGS;
@@ -278 +406 @@
     // generate PCM data storage structure
     pmPCMdata *pcm = pmPCMdataAlloc (params, constraint->paramMask, source);
-
-    pcm->psf = psf;
     pcm->modelConv = psMemIncrRefCounter(model);
     pcm->constraint = constraint;
+
+    pcm->poissonErrors = fitOptions->poissonErrors;
+    pcm->nsigma = fitOptions->nsigma;
 
     pcm->nPix = nPix;
@@ -288 +417 @@
 
 # if (USE_1D_GAUSS)
-    pmModel *modelPSF = source->modelPSF;
-    psAssert (modelPSF, "psf model must be defined");
-    
-    psEllipseAxes axes;
-    bool useReff = pmModelUseReff (modelPSF->type);
-    psF32 *PAR = modelPSF->params->data.F32;
-    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
-    
-    float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
-    float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
 
     pcm->use1Dgauss = true;
-    pcm->sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
-    pcm->nsigma = 2.0;
-
-    pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma);
+    if (USE_1D_CACHE) {
+        pcm->smdata = psImageSmoothCacheSetKernel (&pcm->sigma, &pcm->kappa, pcm->nsigma, source->pixels, source->modelPSF);
+    } else {
+        pcm->smdata2d = psImageSmooth2dCacheSetKernel (&pcm->sigma, &pcm->kappa, pcm->nsigma, source->pixels, source->modelPSF);
+    }
+
 # else
+    // make sure we save a cached copy of the psf flux
+    pmSourceCachePSF (source, maskVal);
+
+    // convert the cached cached psf model for this source to a psKernel
+    psKernel *psf = pmPCMkernelFromPSF (source, psfSize);
+    if (!psf) {
+        // NOTE: this only happens if the source is too close to an edge
+        model->flags |= PM_MODEL_STATUS_BADARGS;
+        return NULL;
+    }
+
+# if (USE_DELTA_PSF)
+    psImageInit (psf->image, 0.0);
+    psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0;
+# endif
+    pcm->psf = psf;
     pcm->smdata = NULL;
 # endif
@@ -341 +478 @@
     }
 
-    // if we changed the fit mode, we need to update nDOF
-    int nParams = 0;
-    // set parameter mask based on fitting mode
-    switch (fitOptions->mode) {
-      case PM_SOURCE_FIT_NORM:
-        // NORM-only model fits only source normalization (Io)
-        nParams = 1;
-        psVectorInit (pcm->constraint->paramMask, 1);
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        break;
-      case PM_SOURCE_FIT_PSF:
-        // PSF model only fits x,y,Io
-        nParams = 3;
-        psVectorInit (pcm->constraint->paramMask, 1);
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
-        break;
-      case PM_SOURCE_FIT_EXT:
-        // EXT model fits all params (except sky)
-        nParams = model->params->n - 1;
-        psVectorInit (pcm->constraint->paramMask, 0);
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
-        break;
-      case PM_SOURCE_FIT_INDEX:
-        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
-        psVectorInit (pcm->constraint->paramMask, 1);
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
-        if (model->params->n == 7) {
-            nParams = 1;
-        } else {
-            nParams = 2;
-            pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
-        }
-        break;
-      case PM_SOURCE_FIT_NO_INDEX:
-        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
-        psVectorInit (pcm->constraint->paramMask, 0);
-        pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
-        if (model->params->n == 7) {
-            nParams = model->params->n - 1;
-        } else {
-            nParams = model->params->n - 2;
-            pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
-        }
-        break;
-      default:
-        psAbort("invalid fitting mode");
-    }
+    int nParams = pmPCMsetParams (pcm->constraint, fitOptions->mode);
 
     if (pcm->nPix <  nParams + 1) {
@@ -415 +504 @@
             pcm->dmodelsConvFlux->data[n] = psImageCopy (pcm->dmodelsConvFlux->data[n], source->pixels, PS_TYPE_F32);
         }
-        psFree(pcm->smdata);
-        pcm->smdata = psImageSmooth_PreAlloc_DataAlloc (source->pixels, pcm->sigma, pcm->nsigma);
+
+        // If we have changed the window, we need to redefine the smoothing target vectors (but pcm->sigma,kappa,nsigma remain)
+        if (USE_1D_CACHE) {
+            psFree(pcm->smdata);
+            pcm->smdata = psImageSmoothCacheAlloc (source->pixels, pcm->sigma, pcm->nsigma);
+
+            pmModel *modelPSF = source->modelPSF;
+            if (modelPSF->type == modelType_GAUSS) {
+                psImageSmoothCacheKernel_Gauss (pcm->smdata, pcm->sigma);
+            }
+            if (modelPSF->type == modelType_PS1_V1) {
+                psImageSmoothCacheKernel_PS1_V1 (pcm->smdata, pcm->sigma, pcm->kappa);
+            }
+        } else {
+            psFree(pcm->smdata2d);
+            pcm->smdata2d = psImageSmooth2dCacheAlloc (pcm->nsigma);
+
+            pmModel *modelPSF = source->modelPSF;
+            if (modelPSF->type == modelType_GAUSS) {
+                // psImageSmooth2dCacheKernel_Gauss (pcm->smdata2d, pcm->sigma);
+            }
+            if (modelPSF->type == modelType_PS1_V1) {
+                psImageSmooth2dCacheKernel_PS1_V1 (pcm->smdata2d, pcm->sigma, pcm->kappa);
+            }
+        }
     }
 
@@ -423 +535 @@
 
 // construct a realization of the source model
-bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize) {
+bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize, float nsigma) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -440 +552 @@
     // convolve the model image with the PSF
     if (USE_1D_GAUSS) {
-        // do not use the threaded, mask-aware version of this code (psImageSmoothMaskPixelsThread):
-        // * the model flux is not masked
-        // * threading takes place above this level
         
-        // define the Gauss parameters from the psf
-        pmModel *modelPSF = source->modelPSF;
-        psAssert (modelPSF, "psf model must be defined");
-    
-        psEllipseAxes axes;
-        bool useReff = pmModelUseReff (modelPSF->type);
-        psF32 *PAR = modelPSF->params->data.F32;
-        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
-    
-        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
-        float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
-
-        float sigma = 0.5 * (FWHM_MAJOR + FWHM_MINOR) / 2.35;
-        float nsigma = 2.0;
-
-        psImageSmooth (source->modelFlux, sigma, nsigma);
+        float sigma = NAN;
+        float kappa = NAN;
+
+        if (USE_1D_CACHE) {
+            psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, nsigma, source->modelFlux, source->modelPSF);
+            psImageSmoothCache_F32 (source->modelFlux, smdata);
+            psFree (smdata);
+        } else {
+            psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, nsigma, source->modelFlux, source->modelPSF);
+            psImageSmooth2dCache_F32 (source->modelFlux, smdata);
+            psFree (smdata);
+        }
+        // old call: psImageSmooth (source->modelFlux, sigma, nsigma);
     } else {
         // make sure we save a cached copy of the psf flux
@@ -478 +584 @@
 }
 
+// construct a realization of the source model
+bool pmPCMMakeModel (pmSource *source, pmModel *model, float Nsigma, psImageMaskType maskVal, int psfSize) {
+
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    // if we already have a cached image, re-use that memory
+    source->modelFlux = psImageCopy (source->modelFlux, source->pixels, PS_TYPE_F32);
+    psImageInit (source->modelFlux, 0.0);
+
+    // modelFlux always has unity normalization (I0 = 1.0)
+    // pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+    pmModelAdd (source->modelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_SKY | PM_MODEL_OP_NORM, maskVal);
+
+    // convolve the model image with the PSF
+    if (USE_1D_GAUSS) {
+
+        float sigma = NAN;
+        float kappa = NAN;
+
+        if (USE_1D_CACHE) {
+            psImageSmoothCacheData *smdata = psImageSmoothCacheSetKernel (&sigma, &kappa, Nsigma, source->modelFlux, source->modelPSF);
+            psImageSmoothCache_F32 (source->modelFlux, smdata);
+            psFree (smdata);
+        } else {
+            psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheSetKernel (&sigma, &kappa, Nsigma, source->modelFlux, source->modelPSF);
+            psImageSmooth2dCache_F32 (source->modelFlux, smdata);
+            psFree (smdata);
+        }
+        // old call: psImageSmooth (source->modelFlux, sigma, nsigma);
+    } else {
+        // make sure we save a cached copy of the psf flux
+        pmSourceCachePSF (source, maskVal);
+
+        // convert the cached cached psf model for this source to a psKernel
+        psKernel *psf = pmPCMkernelFromPSF (source, psfSize);
+        if (!psf) {
+            // NOTE: this only happens if the source is too close to an edge
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+            return NULL;
+        }
+
+        // XXX not sure if I can place the output on top of the input
+        psImageConvolveFFT (source->modelFlux, source->modelFlux, NULL, 0, psf);
+    }
+    return true;
+}

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmPCMdata.h

@@ -14 +14 @@
 /// @addtogroup Objects Object Detection / Analysis Functions
 /// @{
+
+// XXX this is basically for testing -- when I am happy with the convolution process, I'll strip this out
+# define USE_1D_CACHE 0
+# define USE_1D_GAUSS 1
 
 /** pmPCMdata : PSF Convolved Model data storage structure
@@ -36 +40 @@
     int nDOF;
 
+    bool poissonErrors;
+
     bool use1Dgauss;
+    float kappa;
     float sigma;
     float nsigma;
 
-    psImageSmooth_PreAlloc_Data *smdata;
+    // psArray *smdata;
+    psImageSmoothCacheData *smdata;
+    psImageSmooth2dCacheData *smdata2d;
 } pmPCMdata;
 
@@ -96 +105 @@
 bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
 
-bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize);
+bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize, float nsigma);
+
+bool pmPCMMakeModel (pmSource *source, pmModel *model, float Nsigma, psImageMaskType maskVal, int psfSize);
 
 /// @}

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -54 +54 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -36 +36 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -36 +36 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceUtils.h"

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

@@ -34 +34 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -75 +76 @@
             psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
             psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
-            return false;
+            continue;
         }
 
@@ -118 +119 @@
     psfTry->psf->nPSFstars = Npsf;
 
-    // DEBUG code: save the PSF model fit data in detail
-# ifdef DEBUG
-
-    char filename[64];
-    snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy);
-    FILE *f = fopen (filename, "w");
-    psAssert (f, "failed open");
-
-    for (int i = 0; i < psfTry->sources->n; i++) {
-
-        // skip masked sources
-        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
-
-        pmSource *source = psfTry->sources->data[i];
-
-        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
-                 source->peak->xf, source->peak->yf, 
-                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], 
-                 source->psfMag, source->apMag, source->psfMagErr,
-                 source->modelPSF->params->data.F32[PM_PAR_I0], 
-                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], 
-                 source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_7]);
-    }
-    fclose (f);
-# endif
-
     pmSourceVisualShowModelFits (psfTry->psf, psfTry->sources, maskVal);
 
-    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, psfTry->sources->n);
+    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources (%d x %d model)\n", psTimerMark ("psf.fit"), Npsf, psfTry->sources->n, psfTry->psf->trendNx, psfTry->psf->trendNy);
     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, psfTry->sources->n);
 
     if (Npsf == 0) {
-        psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built.");
-        return false;
+#if 0
+        // DEBUG code: save the PSF model fit data in detail
+
+        char hostname[256];
+        gethostname (hostname, 256);
+
+        int pid = getpid();
+
+        char filename[64];
+        snprintf (filename, 64, "psffit.%s.%d.%dx%d.dat", hostname, pid, psfTry->psf->trendNx, psfTry->psf->trendNy);
+        FILE *f = fopen (filename, "w");
+        psAssert (f, "failed open");
+
+        for (int i = 0; i < psfTry->sources->n; i++) {
+
+            // skip masked sources
+          // if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
+
+            pmSource *source = psfTry->sources->data[i];
+
+            if (!source->modelPSF) continue;
+
+            float par7 = (source->modelPSF->params->n == 7) ? -100 : source->modelPSF->params->data.F32[PM_PAR_7];
+            fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f %d\n",
+                     source->peak->xf, source->peak->yf, 
+                     source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], 
+                     source->psfMag, source->apMag, source->psfMagErr,
+                     source->modelPSF->params->data.F32[PM_PAR_I0], 
+                     source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], 
+                     source->modelPSF->params->data.F32[PM_PAR_SYY], par7,
+                     psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
+        }
+        fclose (f);
+#endif
+        psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built.");
+        return false;
     }
 

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

@@ -35 +35 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -35 +35 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -36 +36 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -111 +112 @@
 
     // set the max order (0 = constant) which the number of psf stars can support:
+    int MaxOrderForStars = 0;
+
+    // we require only 3 stars for n = 0, increase stars / cell for higher order
+    if (sources->n >=  16) MaxOrderForStars = 1; //  4 cells, 4 per cell
+    if (sources->n >=  54) MaxOrderForStars = 2; //  9 cells, 6 per cell
+    if (sources->n >= 128) MaxOrderForStars = 3; // 16 cells, 8 per cell
+    if (sources->n >= 300) MaxOrderForStars = 4; // 25 cells, 12 per cell
+    if (sources->n >  576) MaxOrderForStars = 5; // 36 cells, 16 per cell
+
     // rule of thumb: require 3 stars per 'cell' (order+1)^2
-    int MaxOrderForStars = 0;
-    if (sources->n >= 12) MaxOrderForStars = 1; // 4 cells
-    if (sources->n >= 27) MaxOrderForStars = 2; // 9 cells
-    if (sources->n >= 48) MaxOrderForStars = 3; // 16 cells
-    if (sources->n >  75) MaxOrderForStars = 4; // 25 cells
+    // if (sources->n >= 12) MaxOrderForStars = 1; // 4 cells
+    // if (sources->n >= 27) MaxOrderForStars = 2; // 9 cells
+    // if (sources->n >= 48) MaxOrderForStars = 3; // 16 cells
+    // if (sources->n >= 75) MaxOrderForStars = 4; // 25 cells
+    // if (sources->n > 108) MaxOrderForStars = 5; // 36 cells
 
     int orderMax = PS_MAX (options->psfTrendNx, options->psfTrendNy);
@@ -236 +246 @@
     // linear clipped fit of chisq trend vs flux
     if (options->chiFluxTrend) {
+
+        if (0) {
+            FILE *f = fopen ("test.psf.dat", "w");
+            int fd = fileno (f);
+            p_psVectorPrint (fd, flux, "flux");
+            p_psVectorPrint (fd, chisq, "chisq");
+            p_psVectorPrint (fd, mask, "mask");
+            fclose (f);
+        }
+
         bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats,
                                                   mask, 0xff, chisq, NULL, flux);

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

@@ -36 +36 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 

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

@@ -41 +41 @@
 #include "pmSourcePhotometry.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 
@@ -66 +67 @@
     psFree(tmp->extpars);
     psFree(tmp->diffStats);
+    psFree(tmp->galaxyFits);
     psFree(tmp->radialAper);
+    psFree(tmp->lensingOBJ);
+    psFree(tmp->lensingPSF);
     psTrace("psModules.objects", 10, "---- end ----\n");
 }
@@ -164 +168 @@
     source->extpars = NULL;
     source->diffStats = NULL;
+    source->galaxyFits = NULL;
+    source->lensingOBJ = NULL;
+    source->lensingPSF = NULL;
     source->radialAper = NULL;
     source->parent = NULL;
@@ -247 +254 @@
 
     source->region           = in->region;
+
+    // XXX I am not copying the pointers to things like the blends, satstar profile, galaxyFits, etc
 
     return(source);
@@ -690 +699 @@
             // why do we recalculate moments here?
             // we already attempt to do this in psphotSourceStats
-            // pmSourceMoments (source, INNER_RADIUS);
             Nsatstar ++;
             continue;
@@ -804 +812 @@
     return true;
 }
-
-/******************************************************************************
-pmSourceMoments(source, radius): this function takes a subImage defined in the
-pmSource data structure, along with the peak location, and determines the
-various moments associated with that peak.
-
-Requires the following to have been created:
-    pmSource
-    pmSource->peak
-    pmSource->pixels
-    pmSource->variance
-    pmSource->mask
-
-XXX: The peak calculations are done in image coords, not subImage coords.
-
-XXX EAM : this version clips input pixels on S/N
-XXX EAM : this version returns false for several reasons
-*****************************************************************************/
-# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
-
-/*** this been moved to pmSourceMoments.c ***/
-# if (0)
-bool pmSourceMoments(pmSource *source,
-                     psF32 radius)
-{
-    psTrace("psModules.objects", 10, "---- begin ----\n");
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
-
-    //
-    // XXX: Verify the setting for sky if source->moments == NULL.
-    //
-    psF32 sky = 0.0;
-    if (source->moments == NULL) {
-        source->moments = pmMomentsAlloc();
-    } else {
-        sky = source->moments->Sky;
-    }
-
-    //
-    // Sum = SUM (z - sky)
-    // X1  = SUM (x - xc)*(z - sky)
-    // X2  = SUM (x - xc)^2 * (z - sky)
-    // XY  = SUM (x - xc)*(y - yc)*(z - sky)
-    //
-    psF32 peakPixel = -PS_MAX_F32;
-    psS32 numPixels = 0;
-    psF32 Sum = 0.0;
-    psF32 Var = 0.0;
-    psF32 X1 = 0.0;
-    psF32 Y1 = 0.0;
-    psF32 X2 = 0.0;
-    psF32 Y2 = 0.0;
-    psF32 XY = 0.0;
-    psF32 x  = 0;
-    psF32 y  = 0;
-    psF32 R2 = PS_SQR(radius);
-
-    psF32 xPeak = source->peak->x;
-    psF32 yPeak = source->peak->y;
-    psF32 xOff = source->pixels->col0 - source->peak->x;
-    psF32 yOff = source->pixels->row0 - source->peak->y;
-
-    // XXX why do I get different results for these two methods of finding Sx?
-    // XXX Sx, Sy would be better measured if we clip pixels close to sky
-    // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
-    // We loop through all pixels in this subimage (source->pixels), and for each
-    // pixel that is not masked, AND within the radius of the peak pixel, we
-    // proceed with the moments calculation.  need to do two loops for a
-    // numerically stable result.  first loop: get the sums.
-    // XXX EAM : mask == 0 is valid
-
-    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
-        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
-
-        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
-            if (vMsk) {
-                if (*vMsk) {
-                    vMsk++;
-                    psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to mask: %d\n",
-                            col, row, (int)*vMsk);
-                    continue;
-                }
-                vMsk++;
-            }
-            if (isnan(*vPix)) continue;
-
-            psF32 xDiff = col + xOff;
-            psF32 yDiff = row + yOff;
-
-            // radius is just a function of (xDiff, yDiff)
-            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
-#if 1
-                psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to position: %f %f\n",
-                        col, row, xDiff, yDiff);
-#endif
-                continue;
-            }
-
-            psF32 pDiff = *vPix - sky;
-            psF32 wDiff = *vWgt;
-
-            // XXX EAM : check for valid S/N in pixel
-            // XXX EAM : should this limit be user-defined?
-#if 1
-            if (PS_SQR(pDiff) < wDiff) {
-                psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to insignificance: %f, %f\n",
-                        col, row, pDiff, wDiff);
-                continue;
-            }
-#endif
-
-            Var += wDiff;
-            Sum += pDiff;
-
-            psF32 xWght = xDiff * pDiff;
-            psF32 yWght = yDiff * pDiff;
-
-            X1  += xWght;
-            Y1  += yWght;
-
-            XY  += xDiff * yWght;
-            X2  += xDiff * xWght;
-            Y2  += yDiff * yWght;
-
-            peakPixel = PS_MAX (*vPix, peakPixel);
-            numPixels++;
-        }
-    }
-
-    // if we have less than (1/4) of the possible pixels, force a retry
-    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
-    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
-        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n",
-                 numPixels, (int)(0.75*R2), Sum);
-        psTrace("psModules.objects", 10, "---- end (false) ----\n");
-        return (false);
-    }
-
-    psTrace ("psModules.objects", 4, "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
-             sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
-
-    //
-    // first moment X  = X1/Sum + xc
-    // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
-    // Sxy             = XY / Sum
-    //
-    x = X1/Sum;
-    y = Y1/Sum;
-    if ((fabs(x) > radius) || (fabs(y) > radius)) {
-        psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n",
-                 source->peak->x, source->peak->y);
-        psTrace("psModules.objects", 10, "---- end(false)  ----\n");
-        return (false);
-    }
-
-    source->moments->Mx = x + xPeak;
-    source->moments->My = y + yPeak;
-
-    // XXX EAM : Sxy needs to have x*y subtracted
-    source->moments->Mxy = XY/Sum - x*y;
-    source->moments->Sum = Sum;
-    source->moments->SN  = Sum / sqrt(Var);
-    source->moments->Peak = peakPixel;
-    source->moments->nPixels = numPixels;
-
-    // XXX EAM : these values can be negative, so we need to limit the range
-    // XXX EAM : make the use of this consistent: should this be the second moment or sqrt?
-    // source->moments->Mxx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
-    // source->moments->Myy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
-    source->moments->Mxx = PS_MAX(X2/Sum - PS_SQR(x), 0);
-    source->moments->Myy = PS_MAX(Y2/Sum - PS_SQR(y), 0);
-
-    psTrace ("psModules.objects", 4,
-             "sky: %f  Sum: %f  Mx: %f  My: %f  Mxx: %f  Myy: %f  Mxy: %f\n",
-             sky, Sum, source->moments->Mx, source->moments->My,
-             source->moments->Mxx, source->moments->Myy, source->moments->Mxy);
-
-    psTrace("psModules.objects", 10, "---- end ----\n");
-    return(true);
-}
-# endif
 
 // construct a realization of the source model

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSource.h

@@ -40 +40 @@
     PM_SOURCE_TMPF_PETRO_KEEP        = 0x0100,
     PM_SOURCE_TMPF_PETRO_SKIP        = 0x0200,
+    PM_SOURCE_TMPF_EXT_FIT           = 0x0400,  // not just galaxies (trails as well)
+    PM_SOURCE_TMPF_PETRO             = 0x0800,
 } pmSourceTmpF;
 
@@ -117 +119 @@
     pmSourceExtendedPars *extpars;      ///< extended source parameters
     pmSourceDiffStats *diffStats;       ///< extra parameters for difference detections
+    pmSourceGalaxyFits *galaxyFits;     ///< fits to galaxy models (psphotFullForce only)
+    pmSourceLensing *lensingOBJ;        ///< lensing moments parameters (per object)
+    pmSourceLensing *lensingPSF;        ///< lensing moments parameters (psf, interpolated)
     psArray *radialAper;                ///< radial flux in circular apertures
     pmSource *parent;                   ///< reference to the master source from which this is derived

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

@@ -40 +40 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 

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

@@ -286 +286 @@
     return pars;
 }
+
+// *** pmSourceExtFitPars describes extra metadata related to an extended fit
+static void pmSourceGalaxyFitsFree (pmSourceGalaxyFits *tmp) {
+  
+    psFree (tmp->Flux);
+    psFree (tmp->dFlux);
+    psFree (tmp->chisq);
+
+    return;
+}
+
+pmSourceGalaxyFits *pmSourceGalaxyFitsAlloc (void) {
+
+    pmSourceGalaxyFits *tmp = (pmSourceGalaxyFits *) psAlloc(sizeof(pmSourceGalaxyFits));
+    psMemSetDeallocator(tmp, (psFreeFunc) pmSourceGalaxyFitsFree);
+
+    tmp->Flux  = psVectorAllocEmpty (25, PS_TYPE_F32);
+    tmp->dFlux = psVectorAllocEmpty (25, PS_TYPE_F32);
+    tmp->chisq = psVectorAllocEmpty (25, PS_TYPE_F32);
+    tmp->nPix = 0;
+
+    return tmp;
+}

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceExtendedPars.h

@@ -82 +82 @@
 } pmSourceExtFitPars;
 
+typedef struct {
+  psVector *Flux;
+  psVector *dFlux;
+  psVector *chisq;
+  int nPix;
+} pmSourceGalaxyFits;
+
 pmSourceRadialFlux *pmSourceRadialFluxAlloc();
 bool psMemCheckSourceRadialFlux(psPtr ptr);
@@ -109 +116 @@
 pmSourceExtFitPars *pmSourceExtFitParsAlloc (void);
 
+pmSourceGalaxyFits *pmSourceGalaxyFitsAlloc (void);
+
 /// @}
 # endif /* PM_SOURCE_H */

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

@@ -40 +40 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourcePhotometry.h"
@@ -59 +60 @@
     opt->maxTol = 1.00;
     opt->weight = 1.00;
+    opt->nsigma = 5.00;
     opt->maxChisqDOF = NAN;
     opt->poissonErrors = true;
@@ -66 +68 @@
     opt->gainFactorMode = 0;
     opt->chisqConvergence = true;
+    opt->isInteractive = false;
 
     return opt;
@@ -247 +250 @@
     myMin->gainFactorMode = options->gainFactorMode;
     myMin->chisqConvergence = options->chisqConvergence;
+    myMin->isInteractive = options->isInteractive;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
@@ -279 +283 @@
     // set the model success or failure status
     model->flags |= PM_MODEL_STATUS_FITTED;
-    if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
+    if (!fitStatus) {
+        if (isnan(myMin->value)) {
+          model->flags |= PM_MODEL_STATUS_NAN_CHISQ;
+        } else {
+          model->flags |= PM_MODEL_STATUS_NONCONVERGE;
+        }
+    }
 
     if (myMin->chisqConvergence) {

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceFitModel.h

@@ -21 +21 @@
     PM_SOURCE_FIT_EXT_AND_SKY,
     PM_SOURCE_FIT_INDEX,
+    PM_SOURCE_FIT_SHAPE,
     PM_SOURCE_FIT_NO_INDEX,
     PM_SOURCE_FIT_TRAIL,
@@ -33 +34 @@
     float weight;                       ///< use this weight for constant-weight fits
     float covarFactor;                  ///< covariance factor for calculating the chisq
+    float nsigma;                       ///< how far out to convolve
     bool poissonErrors;                 ///< use poisson errors for fits?
     bool saveCovariance;
     int gainFactorMode;
     bool chisqConvergence; 
+    bool isInteractive;
 } pmSourceFitOptions;
 

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

@@ -38 +38 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourcePhotometry.h"
@@ -51 +52 @@
 # define TIMING 0
 
+bool pmSourceChisqModelFlux (pmSource *source, pmModel *model, psImageMaskType maskVal);
+
 bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
     
@@ -68 +71 @@
     myMin->chisqConvergence = fitOptions->chisqConvergence;
     myMin->gainFactorMode = fitOptions->gainFactorMode;
+    myMin->isInteractive = fitOptions->isInteractive;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
@@ -112 +116 @@
     } else {
         // xxx this is wrong because it does not convolve with the psf
-        pmSourceChisqUnsubtracted (source, pcm->modelConv, maskVal);
+        pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
+        pmSourceChisqModelFlux (source, pcm->modelConv, maskVal);
     }
     if (TIMING) { t4 = psTimerMark ("pmSourceFitPCM"); }
@@ -118 +123 @@
     // set the model success or failure status
     pcm->modelConv->flags |= PM_MODEL_STATUS_FITTED;
-    if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE;
+
+    if (!fitStatus) {
+        if (isnan(myMin->value)) {
+            pcm->modelConv->flags |= PM_MODEL_STATUS_NAN_CHISQ;
+        } else {
+            pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE;
+        }
+    }
 
     if (myMin->chisqConvergence) {

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

@@ -39 +39 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourcePhotometry.h"
@@ -352 +353 @@
         // set the model success or failure status
         model->flags |= PM_MODEL_STATUS_FITTED;
-        if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
+        if (!fitStatus) {
+          if (isnan(myMin->value)) {
+            model->flags |= PM_MODEL_STATUS_NAN_CHISQ;
+          } else {
+            model->flags |= PM_MODEL_STATUS_NONCONVERGE;
+          }
+        }
 
         // models can go insane: reject these
@@ -570 +577 @@
     myMin->gainFactorMode = options->gainFactorMode;
     myMin->chisqConvergence = options->chisqConvergence;
+    myMin->isInteractive = options->isInteractive;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);

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

@@ -25 +25 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -47 +47 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -61 +62 @@
 static bool pmReadoutReadXFIT(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex);
 static bool pmReadoutReadXRAD(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex);
+static bool pmReadoutReadXGAL(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex);
 
 // lookup the EXTNAME values used for table data and image header segments
@@ -69 +71 @@
                         psString *xfitname,    // Extension name for extended fitted measurements
                         psString *xradname,    // Extension name for radial apertures
+                        psString *xgalname,    // Extension name for galaxy shapes
                         const pmFPAfile *file, // File of interest
                         const pmFPAview *view  // View to level of interest
@@ -140 +143 @@
         }
         *xradname = pmFPAfileNameFromRule (rule, file, view);
+    }
+
+    // EXTNAME for radial apertures
+    if (xgalname) {
+        const char *rule = psMetadataLookupStr(&status, menu, "CMF.XGAL");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, true, "missing entry for CMF.XGAL in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        *xgalname = pmFPAfileNameFromRule (rule, file, view);
     }
 
@@ -362 +375 @@
             status &= pmSourcesWrite_##TYPE##_XRAD (file->fits, readout, sources, file->header, xradname, recipe); \
         }                                                               \
+        if (xgalname) {                                                 \
+            status &= pmSourcesWrite_##TYPE##_XGAL (file->fits, readout, sources, xgalname, recipe); \
+        }                                                               \
     }
 
@@ -464 +480 @@
         }
 
-        // if this is not TRUE, the output files only contain the psf measurements.
-        bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS");
+        // if none of these are TRUE, the output files only contain the psf measurements.
+        bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
+        bool doAnnuli    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
+        bool XSRC_OUTPUT = doPetrosian || doAnnuli;
         bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS");
         bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES");
+        bool XGAL_OUTPUT = psMetadataLookupBool(&status, recipe, "GALAXY_SHAPES");
 
         // define the EXTNAME values for the different data segments:
@@ -476 +495 @@
         psString xfitname = NULL;
         psString xradname = NULL;
+        psString xgalname = NULL;
         if (!pmSourceIOextnames(&headname, &dataname, &deteffname, 
                                 XSRC_OUTPUT ? &xsrcname : NULL,
                                 XFIT_OUTPUT ? &xfitname : NULL, 
                                 XRAD_OUTPUT ? &xradname : NULL, 
+                                XGAL_OUTPUT ? &xgalname : NULL, 
                                 file, view)) {
             return false;
@@ -563 +584 @@
                 psMetadataAddStr (outhead, PS_LIST_TAIL, "XRADNAME", PS_META_REPLACE, "name of XRAD table extension", xradname);
             }
+            if (xgalname) {
+                psMetadataAddStr (outhead, PS_LIST_TAIL, "XGALNAME", PS_META_REPLACE, "name of XGAL table extension", xgalname);
+            }
 
             // these are case-sensitive since the EXTYPE is case-sensitive
@@ -574 +598 @@
             PM_SOURCES_WRITE("PS1_V3",    CMF_PS1_V3);
             PM_SOURCES_WRITE("PS1_V4",    CMF_PS1_V4);
+            PM_SOURCES_WRITE("PS1_V5",    CMF_PS1_V5);
             PM_SOURCES_WRITE("PS1_SV1",   CMF_PS1_SV1);
             PM_SOURCES_WRITE("PS1_SV2",   CMF_PS1_SV2);
@@ -609 +634 @@
         psFree (xfitname);
         psFree (xradname);
+        psFree (xgalname);
         psFree (deteffname);
 
@@ -620 +646 @@
         psFree (xfitname);
         psFree (xradname);
+        psFree (xgalname);
         psFree (deteffname);
         return false;
 
+      case PM_FPA_FILE_CFF: {
+        // determine the output table format
+        psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT");
+        if (!status) {
+            psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data");
+            return false;
+        }
+
+        hdu = pmFPAviewThisHDU (view, fpa);
+        pmConfigConformHeader(hdu->header, file->format);
+        psFitsWriteBlank (file->fits, hdu->header, NULL);
+        file->header = hdu->header;
+        file->wrote_phu = true;
+        if (!pmSourcesWrite_CFF(readout, file->fits, sources, hdu->header, recipe)) {
+            psError(PS_ERR_UNKNOWN, false, "failed to write CFF");
+            return false;
+        }
+        break;
+      }
+        
       default:
         fprintf (stderr, "warning: type mismatch\n");
@@ -914 +961 @@
     psArray *sources = NULL;
     pmHDU *hdu;
+
+    // define the EXTNAME values for the different data segments:
+    psString headname = NULL;
+    psString dataname = NULL;
+    psString deteffname = NULL;
+    psString xsrcname = NULL;
+    psString xfitname = NULL;
+    psString xradname = NULL;
+    psString xgalname = NULL;
+
+    psMetadata *tableHeader = NULL;
+    char *xtension = NULL;
 
     switch (file->type) {
@@ -963 +1022 @@
         hdu = pmFPAviewThisHDU (view, file->fpa);
 
-        // define the EXTNAME values for the different data segments:
-        psString headname = NULL;
-        psString dataname = NULL;
-        psString deteffname = NULL;
-        psString xsrcname = NULL;
-        psString xfitname = NULL;
-        psString xradname = NULL;
-
         // determine the output table format. Assume if we need to output extendend source
         // parameters that they may exist in the input. 
@@ -981 +1032 @@
         }
 
-        // if this is not TRUE, the output files only contain the psf measurements.
-        bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS");
+        // if none of these are TRUE, we only read the psf measurements
+        // XXX: shouldn't we look for these extensions and read the regardless of the recipe values?
+        bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
+        bool doAnnuli    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
+        bool XSRC_OUTPUT = doPetrosian || doAnnuli;
         bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS");
         bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES");
+        bool XGAL_OUTPUT = psMetadataLookupBool(&status, recipe, "GALAXY_SHAPES");
 
         if (!pmSourceIOextnames(&headname, &dataname, &deteffname, 
@@ -990 +1045 @@
                 XFIT_OUTPUT ? &xfitname : NULL, 
                 XRAD_OUTPUT ? &xradname : NULL,
+                XGAL_OUTPUT ? &xgalname : NULL,
                 file, view)) {
             return false;
@@ -1033 +1089 @@
         }
 
-        psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header
+        tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header
         if (!tableHeader) psAbort("cannot read table header");
 
-        char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
+        xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
         if (!xtension) psAbort("cannot read table type");
         if (strcmp (xtension, "BINTABLE")) {
@@ -1062 +1118 @@
             PM_SOURCES_READ_PSF("PS1_V3",    CMF_PS1_V3);
             PM_SOURCES_READ_PSF("PS1_V4",    CMF_PS1_V4);
+            PM_SOURCES_READ_PSF("PS1_V5",    CMF_PS1_V5);
             PM_SOURCES_READ_PSF("PS1_SV1",   CMF_PS1_SV1);
             PM_SOURCES_READ_PSF("PS1_SV2",   CMF_PS1_SV2);
@@ -1069 +1126 @@
 
             long *sourceIndex = NULL;
-            if (XSRC_OUTPUT || XFIT_OUTPUT || XRAD_OUTPUT) {
+            if (XSRC_OUTPUT || XFIT_OUTPUT || XRAD_OUTPUT || XGAL_OUTPUT) {
+                // Build sourceIndex. Lookup table from source->seq to index in sources array.
+                // Consists of an array of length max(source->seq) + 1.
+
+                // find maximum sequence number
                 long seq_max = -1;
                 for (long i = sources->n -1; i >= 0; i--) {
@@ -1082 +1143 @@
                     }
                 }
+                // allocate and initialize the index
                 sourceIndex = psAlloc((seq_max + 1) * sizeof(long));
                 for (long i = 0; i < seq_max; i++) {
                     sourceIndex[i] = -1;
                 }
+                // populate the index
                 for (long i = 0; i < sources->n; i++) {
                     pmSource *source = sources->data[i];
@@ -1112 +1175 @@
                 psFree(xradname);
             }
+            if (XGAL_OUTPUT && xgalname) {
+                // a cmf file may have an XGAL extension, but it is not required
+                if (!pmReadoutReadXGAL(file, readout, exttype, hdu->header, xgalname, sources, sourceIndex)) {
+                    // do anything?
+                }
+                psFree(xgalname);
+            }
             psFree(sourceIndex);
 
@@ -1132 +1202 @@
         break;
 
+      case PM_FPA_FILE_CFF:
+        // read in header, if not yet loaded
+        hdu = pmFPAviewThisHDU (view, file->fpa);
+
+        // look these up in the camera config?
+        // headrule = {CHIP.NAME}.hdr
+        // datarule = {CHIP.NAME}.cff
+
+        // define the EXTNAME values for the different data segments:
+        headname = pmFPAfileNameFromRule("{CHIP.NAME}.hdr", file, view);
+        dataname = pmFPAfileNameFromRule("{CHIP.NAME}.cff", file, view);
+
+        // advance to the IMAGE HEADER extension
+        if (hdu->header == NULL) {
+            // if the IMAGE header does not exist, we have no data for this view
+            if (!psFitsMoveExtNameClean (file->fits, headname)) {
+                readout->data_exists = false;
+                psFree (headname);
+                psFree (dataname);
+                return true;
+            }
+            hdu->header = psFitsReadHeader (NULL, file->fits);
+        }
+
+        // advance to the table data extension
+        // since we have read the IMAGE header, the TABLE header should exist
+        if (!psFitsMoveExtName (file->fits, dataname)) {
+            psAbort("cannot find data extension %s in %s", dataname, file->filename);
+        }
+
+        tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header
+        if (!tableHeader) psAbort("cannot read table header");
+
+        // verify this is a binary table
+        char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
+        if (!xtension) psAbort("cannot read table type");
+        if (strcmp (xtension, "BINTABLE")) {
+            psWarning ("no binary table in extension %s, skipping\n", dataname);
+            psFree(tableHeader);
+            return false;
+        }
+
+        sources = pmSourcesRead_CFF(file->fits, hdu->header);
+
+        psTrace("psModules.objects", 6, "read CMF table from %s : %s : %s", file->filename, headname, dataname);
+        psFree (headname);
+        psFree (dataname);
+        psFree (tableHeader);
+        break;
+
       default:
         fprintf (stderr, "warning: type mismatch\n");
@@ -1272 +1392 @@
         PM_SOURCES_READ_XSRC("PS1_V3",    CMF_PS1_V3);
         PM_SOURCES_READ_XSRC("PS1_V4",    CMF_PS1_V4);
+        PM_SOURCES_READ_XSRC("PS1_V5",    CMF_PS1_V5);
         PM_SOURCES_READ_XSRC("PS1_SV1",   CMF_PS1_SV1);
         PM_SOURCES_READ_XSRC("PS1_SV2",   CMF_PS1_SV2);
@@ -1311 +1432 @@
         PM_SOURCES_READ_XFIT("PS1_V3",    CMF_PS1_V3);
         PM_SOURCES_READ_XFIT("PS1_V4",    CMF_PS1_V4);
+        PM_SOURCES_READ_XFIT("PS1_V5",    CMF_PS1_V5);
         PM_SOURCES_READ_XFIT("PS1_SV1",   CMF_PS1_SV1);
         PM_SOURCES_READ_XFIT("PS1_SV2",   CMF_PS1_SV2);
@@ -1349 +1471 @@
         PM_SOURCES_READ_XRAD("PS1_V3",    CMF_PS1_V3);
         PM_SOURCES_READ_XRAD("PS1_V4",    CMF_PS1_V4);
+        PM_SOURCES_READ_XRAD("PS1_V5",    CMF_PS1_V5);
         PM_SOURCES_READ_XRAD("PS1_SV1",   CMF_PS1_SV1);
         PM_SOURCES_READ_XRAD("PS1_SV2",   CMF_PS1_SV2);
@@ -1358 +1481 @@
     return status;
 }
+static bool pmReadoutReadXGAL(pmFPAfile *file, pmReadout *readout, char *exttype, psMetadata *hduHeader, psString xgalname, psArray *sources, long *sourceIndex) 
+{
+    if (!psFitsMoveExtNameClean (file->fits, xgalname)) {
+        psTrace ("pmFPAfile", 1, "cannot find xgal extension %s in %s, skipping", xgalname, file->filename);
+        return false;
+    }
+
+    psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header
+    if (!tableHeader) psAbort("cannot read table header");
+
+    char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
+    if (!xtension) psAbort("cannot read table type");
+    if (strcmp (xtension, "BINTABLE")) {
+        psFree(tableHeader);
+        psWarning ("no binary table in extension %s, skipping\n", xgalname);
+        return false;
+    }
+
+# define PM_SOURCES_READ_XGAL(NAME,TYPE)                                \
+    if (!strcmp (exttype, NAME)) {                                      \
+        status = pmSourcesRead_##TYPE##_XGAL(file->fits, readout, hduHeader, tableHeader, sources, sourceIndex); \
+    }                                                                   
+
+    bool status = false;
+    if (file->type == PM_FPA_FILE_CMF) {
+        PM_SOURCES_READ_XGAL("PS1_V1",    CMF_PS1_V1);
+        PM_SOURCES_READ_XGAL("PS1_V2",    CMF_PS1_V2);
+        PM_SOURCES_READ_XGAL("PS1_V3",    CMF_PS1_V3);
+        PM_SOURCES_READ_XGAL("PS1_V4",    CMF_PS1_V4);
+        PM_SOURCES_READ_XGAL("PS1_V5",    CMF_PS1_V5);
+        PM_SOURCES_READ_XGAL("PS1_SV1",   CMF_PS1_SV1);
+        PM_SOURCES_READ_XGAL("PS1_SV2",   CMF_PS1_SV2);
+        PM_SOURCES_READ_XGAL("PS1_DV1",   CMF_PS1_DV1);
+        PM_SOURCES_READ_XGAL("PS1_DV2",   CMF_PS1_DV2);
+        PM_SOURCES_READ_XGAL("PS1_DV3",   CMF_PS1_DV3);
+    }
+    psFree(tableHeader);
+    return status;
+}

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO.h

@@ -21 +21 @@
   bool pmSourcesWrite_##TYPE##_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname); \
   bool pmSourcesWrite_##TYPE##_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe); \
+  bool pmSourcesWrite_##TYPE##_XGAL(psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe); \
   psArray *pmSourcesRead_##TYPE (psFits *fits, psMetadata *header); \
   bool pmSourcesRead_##TYPE##_XSRC (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index); \
   bool pmSourcesRead_##TYPE##_XFIT (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index); \
   bool pmSourcesRead_##TYPE##_XRAD (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index);\
+  bool pmSourcesRead_##TYPE##_XGAL (psFits *fits, pmReadout *readout, psMetadata *header, psMetadata *tableHeader, psArray *sources, long *index);\
   
 // All of these functions need to use the same API, even if not all elements are used in a specific case
@@ -35 +37 @@
 MK_PROTO(CMF_PS1_V3);
 MK_PROTO(CMF_PS1_V4);
+MK_PROTO(CMF_PS1_V5);
 MK_PROTO(CMF_PS1_SV1);
 MK_PROTO(CMF_PS1_SV2);
@@ -52 +55 @@
 
 psArray *pmSourcesReadCMP (char *filename, psMetadata *header);
+psArray *pmSourcesRead_CFF (psFits *fits, psMetadata *header);
+bool pmSourcesWrite_CFF (pmReadout *readout, psFits *fits, psArray *sources, psMetadata *header, psMetadata *recipe);
 
 bool pmSourcesWritePSFs (psArray *sources, char *filename);

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceIO_CMF.c.in

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -178 +179 @@
         @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
         @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
+
+        // Lensing parameters:
+        if (source->lensingOBJ && source->lensingOBJ->smear) {
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X11);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X12);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X22);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e1);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e2);
+        }
+
+        if (source->lensingOBJ && source->lensingOBJ->shear) {
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X11);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X12);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X22);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e1); 
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e2); 
+        }
+
+        if (source->lensingOBJ && source->lensingPSF->smear) {
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X11);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X12);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X22);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e1); 
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e2); 
+        }
+
+        if (source->lensingOBJ && source->lensingPSF->shear) {
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X11);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X12);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X22);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e1); 
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e2); 
+        }
 
         @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
@@ -708 +742 @@
             return false;
         }
-        // Find the source with this sequence number. 
-        // XXX: I am assuming that sources is sorted in order of seq
+        // Find the source with this sequence number using the sourceIndex. 
         long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
-        pmSource *source = NULL;
-#ifndef ASSUME_SORTED
-        long j = seq < sources->n ? seq : sources->n - 1;
-        for (; j >= 0; j--) {
-            source = sources->data[j];
-            if (source->seq == seq) {
-                break;
-            }
-        }
-#else
         long j = sourceIndex[seq];
         psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
-        source = sources->data[j];
-#endif
+        pmSource *source = sources->data[j];
         if (!source) {
             psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
@@ -789 +811 @@
     char name[64];
 
+    pmModelType modelTypeTrail = pmModelClassGetType("PS_MODEL_TRAIL");
+
     // create a header to hold the output data
     psMetadata *outhead = psMetadataAlloc ();
 
+    pmModelClassWriteHeader(outhead);
+
     // write the links to the image header
     psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);
@@ -798 +824 @@
     sources = psArraySort (sources, pmSourceSortByFlux);
 
-    @>PS1_DV2@ float magOffset; 
-    @>PS1_DV2@ float zeroptErr; 
-    @>PS1_DV2@ float fwhmMajor; 
-    @>PS1_DV2@ float fwhmMinor;
-    @>PS1_DV2@ pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
+    float magOffset; 
+    float zeroptErr; 
+    float fwhmMajor; 
+    float fwhmMinor;
+    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
 
     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
@@ -823 +849 @@
     @>PS1_DV2@ pmChip *chip = readout->parent->parent;
 
+    pmModelStatus badModel = PM_MODEL_STATUS_NONE;
+    badModel |= PM_MODEL_STATUS_BADARGS;
+    badModel |= PM_MODEL_STATUS_OFFIMAGE;
+    badModel |= PM_MODEL_STATUS_NAN_CHISQ;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID;
+    badModel |= PM_MODEL_PCM_FAIL_GUESS;
+
     table = psArrayAllocEmpty (sources->n);
 
@@ -847 +881 @@
 
             // skip models which were not actually fitted
-            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
+            // XXX 
+            if (model->flags & badModel) continue;
 
             PAR = model->params->data.F32;
@@ -853 +888 @@
             xPos = PAR[PM_PAR_XPOS];
             yPos = PAR[PM_PAR_YPOS];
-            xErr = dPAR[PM_PAR_XPOS];
-            yErr = dPAR[PM_PAR_YPOS];
+
+            // for the extended source models, we do not always fit the centroid in the non-linear fitting process
+            // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM,
+            // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM)
+            // TRAIL : X,Y are fitted
+            // 
+            
+            // XXX this should be based on what happened, not on the model type
+            if (model->type == modelTypeTrail) {
+                xErr = dPAR[PM_PAR_XPOS];
+                yErr = dPAR[PM_PAR_YPOS];
+            } else {
+                // this is definitely an underestimate since it does not
+                // account for the extent of the source
+                xErr = fwhmMajor * model->magErr / 2.35;
+                yErr = fwhmMinor * model->magErr / 2.35;
+            }
 
             @>PS1_DV2@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
@@ -870 +920 @@
             row = psMetadataAlloc ();
 
-            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
             // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
             // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment"
@@ -888 +937 @@
             @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN;
             @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration",   calMag);
-            @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Magnitude using supplied calibration",   model->chisq);
-            @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Magnitude using supplied calibration",   model->nDOF);
+            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Model Chisq",   model->chisq);
+            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Model num degrees of freedom",   model->nDOF);
+            @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE",    PS_DATA_S32, "type for chosen EXT_MODEL",   source->modelEXT ? source->modelEXT->type : -1);
+
+            // EAM : adding for PV2 outputs:
+            @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", source->modelEXT ? source->modelEXT->flags : 0);
 
             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE",   0, "position angle at source (degrees)",         posAngle);
@@ -946 +999 @@
 
                 snprintf (name, 64, "EXT_PAR_%02d", k);
-
+               
                 if (k < model->params->n) {
                     psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]);
@@ -956 +1009 @@
             // optionally, write out the covariance matrix values
             // XXX do I need to pad this to match the biggest covar matrix?
-            if (model->covar) {
+            if (false && model->covar) {
                 for (int iy = 0; iy < model->covar->numCols; iy++) {
                     for (int ix = iy; ix < model->covar->numCols; ix++) {
@@ -1005 +1058 @@
         return false;
     }
+    // set up the lookup table to translate between input model types and output model types
+    // if not defined it is assumed that the tables are the same
+    pmModelClassReadHeader(tableHeader);
 
     for (long i = 0; i < numSources; i++) {
@@ -1013 +1069 @@
             return false;
         }
-        // Find the source with this sequence number. 
-        // XXX: I am assuming that sources is sorted in order of seq.
         long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
-        long j = seq < sources->n ? seq : sources->n - 1;
-        pmSource *source = NULL;
-        for (; j >= 0; j--) {
-            source = sources->data[j];
-            if (source->seq == seq) {
-                break;
-            }
-        }
+        long j = sourceIndex[seq];
+        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
+        pmSource *source = sources->data[j];
         if (!source) {
             psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
@@ -1057 +1106 @@
         model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG");
         model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG");
+
+        model->chisq = psMetadataLookupF32(&status, row, "EXT_CHISQ");
+        model->nDOF = psMetadataLookupF32(&status, row, "EXT_NDOF");
+
+        // EXT_MODEL_TYPE gives the model chosen by psphot as the best.
+        // Putting this into the XFIT table makes 3 copies of it (one for each model)
+        // but since we have many fewer XFIT rows than psf rows that is cheaper than putting it
+        // in the psf table.
+        psS32 extModelType = psMetadataLookupS32(&status, row, "EXT_MODEL_TYPE");
+        if (status) {
+            // translate between the type value in xfit and values used by this program
+            extModelType = pmModelClassGetLocalType(extModelType);
+        } else {
+            // older cmfs don't have this column
+            extModelType = -1;
+        } 
 
         psEllipseAxes axes;
@@ -1072 +1137 @@
         if (model->params->n > 7) {
             PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07");
-        }
-        // read the covariance matrix
-        int nparams = model->params->n;
-        psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
-        for (int y = 0; y < nparams; y++) {
-            for (int x = 0; x < nparams; x++) {
-                char name[64];
-                snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
-                covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
+            // XXX add an error:
+            // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_");
+        }
+
+        // XXX : make this depend on what is in the cmf
+        if (0) {
+            // read the covariance matrix
+            int nparams = model->params->n;
+            psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
+            for (int y = 0; y < nparams; y++) {
+                for (int x = 0; x < nparams; x++) {
+                    char name[64];
+                    snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
+                    covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
+                }
+            }
+            model->covar = covar;
+        }
+
+        if (modelType == extModelType) {
+            // The software that created this source picked this model as the best of the fits. 
+            // Set the extModel to point to it.
+            // This is important for programs like psastro (skycal) so that its output cmfs
+            // will have valid EXT_MODEL_TYPE
+            psFree(source->modelEXT);
+            source->modelEXT = psMemIncrRefCounter(model);
+            source->type = PM_SOURCE_TYPE_EXTENDED;
+            if (0) {
+                // since FLAGS were read we don't need to do this
+                source->mode |= PM_SOURCE_MODE_EXTMODEL;
+                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
             }
         }
-        model->covar = covar;
 
         psArrayAdd(source->modelFits, 1, model);
         psFree(model);
-
         psFree(row);
     }
@@ -1209 +1294 @@
 
         write_annuli:
-            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",       PS_DATA_VECTOR, "flux within annuli",       radFlux);
-            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_DATA_VECTOR, "flux error in annuli",     radFluxErr);
-            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation",  radFluxStdev);
-            psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",       PS_DATA_VECTOR, "fill factor of annuli",    radFill);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX",       PS_META_REPLACE, "flux within annuli",       radFlux);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_META_REPLACE, "flux error in annuli",     radFluxErr);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation",  radFluxStdev);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL",       PS_META_REPLACE, "fill factor of annuli",    radFill);
             psFree (radFlux);
             psFree (radFluxErr);
@@ -1277 +1362 @@
             return false;
         }
-        // Find the source with this sequence number. 
-        // XXX: I am assuming that sources is sorted in order of seq.
         long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
-        long j = seq < sources->n ? seq : sources->n - 1;
-        pmSource *source = NULL;
-        for (; j >= 0; j--) {
-            source = sources->data[j];
-            if (source->seq == seq) {
-                break;
-            }
-        }
+        long j = sourceIndex[seq];
+        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
+        pmSource *source = sources->data[j];
         if (!source) {
             psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
@@ -1343 +1421 @@
     return true;
 }
+
+// XXX where should I record the number of columns??
+bool pmSourcesWrite_CMF_@CMFMODE@_XGAL (psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe)
+{
+    bool status = false;
+
+    // perform full non-linear fits / extended source analysis?
+    if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) {
+        psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n");
+        return true;
+    }
+
+    // create a header to hold the output data
+    psMetadata *outhead = psMetadataAlloc ();
+
+    // write the links to the image header
+    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "galaxy table extension", extname);
+
+    psMetadataAddStr (outhead, PS_LIST_TAIL, "HI", PS_META_REPLACE, "does this get through?", "THERE");
+
+    // let's write these out in S/N order
+    sources = psArraySort (sources, pmSourceSortByFlux);
+
+    psArray *table = psArrayAllocEmpty (sources->n);
+
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *thisSource = sources->data[i];
+
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+
+        // if we did not fit the galaxy model, modelFits will be NULL
+        if (source->modelFits == NULL) continue;
+
+        // if we did not fit the galaxy model, galaxyFits will also be NULL
+        if (source->galaxyFits == NULL) continue;
+
+        pmModel *model = source->modelFits->data[0];
+        if (!model) return false;
+
+        // X,Y coordinates are stored with the model parameters
+        psF32 *PAR = model->params->data.F32;
+
+        psMetadata *row = psMetadataAlloc ();
+
+        // we write out the x,y positions so people can link to the psf either way (position or ID)
+        psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index", source->seq);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "X_FIT",            0, "model x coordinate",             PAR[PM_PAR_XPOS]);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "Y_FIT",            0, "model y coordinate",             PAR[PM_PAR_YPOS]);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "NPIX",             0, "number of pixels for fits",      source->galaxyFits->nPix);
+
+        psVector *Flux = source->galaxyFits->Flux;
+        psVector *dFlux = source->galaxyFits->dFlux;
+        psVector *chisq = source->galaxyFits->chisq;
+
+        psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX",     PS_META_REPLACE, "normalization for galaxy flux", Flux);
+        psMetadataAddVector (row, PS_LIST_TAIL, "GAL_FLUX_ERR", PS_META_REPLACE, "error on normalization", dFlux);
+        psMetadataAddVector (row, PS_LIST_TAIL, "GAL_CHISQ",    PS_META_REPLACE, "galaxy fit chisq", chisq);
+
+        psArrayAdd (table, 100, row);
+        psFree (row);
+    }
+
+    if (table->n == 0) {
+        if (!psFitsWriteBlank (fits, outhead, extname)) {
+            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
+            psFree(outhead);
+            psFree(table);
+            return false;
+        }
+        psFree (outhead);
+        psFree (table);
+        return true;
+    }
+
+    psTrace ("pmFPAfile", 5, "writing galaxy data %s\n", extname);
+    if (!psFitsWriteTable (fits, outhead, table, extname)) {
+        psError(psErrorCodeLast(), false, "writing galaxy data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
+    }
+    psFree (outhead);
+    psFree (table);
+    return true;
+}
+
+bool pmSourcesRead_CMF_@CMFMODE@_XGAL(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psMetadata *tableHeader, psArray *sources, long *sourceIndex)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(sources, false);
+
+    bool status;
+    long numSources = psFitsTableSize(fits); // Number of sources in table
+    if (numSources == 0) {
+        psError(psErrorCodeLast(), false, "XGAL Table contains no entries\n");
+        return false;
+    }
+
+    for (long i = 0; i < numSources; i++) {
+        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
+        if (!row) {
+            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
+            psFree(row);
+            return false;
+        }
+        // Find the source with this sequence number. 
+        // XXX: I am assuming that sources is sorted in order of seq
+        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
+        long j = sourceIndex[seq];
+        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
+
+        pmSource *source = sources->data[j];
+        if (!source) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
+            psFree(row);
+            return false;
+        }
+
+        psVector *Flux  = psMetadataLookupVector(&status, row, "GAL_FLUX");
+        psVector *dFlux = psMetadataLookupVector(&status, row, "GAL_FLUX_ERR");
+        psVector *chisq = psMetadataLookupVector(&status, row, "GAL_CHISQ");
+
+        if (Flux && Flux->n > 0) {
+            psFree(source->galaxyFits);
+            source->galaxyFits = pmSourceGalaxyFitsAlloc();
+            source->galaxyFits->nPix = psMetadataLookupF32(&status, row, "NPIX");
+
+            psFree(source->galaxyFits->Flux);
+            source->galaxyFits->Flux  = psMemIncrRefCounter(Flux);
+            psFree(source->galaxyFits->dFlux);
+            source->galaxyFits->dFlux = psMemIncrRefCounter(dFlux);
+            psFree(source->galaxyFits->chisq);
+            source->galaxyFits->chisq = psMemIncrRefCounter(chisq);
+        }
+
+        psFree(row);
+    }
+
+    return true;
+}

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -713 +714 @@
     return true;
 }
+
+bool pmSourcesWrite_PS1_CAL_0_XGAL (psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe)
+{
+    return true;
+}

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -255 +256 @@
     return true;
 }
+
+bool pmSourcesWrite_PS1_DEV_0_XGAL(psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe)
+{
+    return true;
+}

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -595 +596 @@
     return true;
 }
+
+bool pmSourcesWrite_PS1_DEV_1_XGAL(psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe)
+{
+    return true;
+}

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -225 +226 @@
     return true;
 } 
+
+bool pmSourcesWrite_SMPDATA_XGAL(psFits *fits, pmReadout *readout, psArray *sources, char *extname, psMetadata *recipe)
+{
+    return true;
+} 

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

@@ -44 +44 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

branches/eam_branches/ps2-tc3-20130727/psModules/src/objects/pmSourceMasks.h

@@ -56 +56 @@
     PM_SOURCE_MODE2_DIFF_SELF_MATCH       = 0x00000800, ///< positive detection match is probably this source 
     PM_SOURCE_MODE2_SATSTAR_PROFILE       = 0x00001000, ///< saturated source is modeled with a radial profile
+
+    PM_SOURCE_MODE2_ECONTOUR_FEW_PTS      = 0x00002000, ///< too few points to measure the elliptical contour
+    PM_SOURCE_MODE2_RADBIN_NAN_CENTER     = 0x00004000, ///< radial bins failed with too many NaN center bin
+    PM_SOURCE_MODE2_PETRO_NAN_CENTER      = 0x00008000, ///< petrosian radial bins failed with too many NaN center bin
+    PM_SOURCE_MODE2_PETRO_NO_PROFILE      = 0x00010000, ///< petrosian not build because radial bins missing
+
+    PM_SOURCE_MODE2_PETRO_INSIG_RATIO     = 0x00020000, ///< insignificant measurement of petrosian ratio
+    PM_SOURCE_MODE2_PETRO_RATIO_ZEROBIN   = 0x00040000, ///< petrosian ratio in the 0th bin (likely bad)
+    
+    PM_SOURCE_MODE2_EXT_FITS_RUN          = 0x00080000, ///< we attempted to run extended fits on this source
+    PM_SOURCE_MODE2_EXT_FITS_FAIL         = 0x00100000, ///< at least one of the model fits failed
+    PM_SOURCE_MODE2_EXT_FITS_RETRY        = 0x00200000, ///< one of the model fits was re-tried with new window
+    PM_SOURCE_MODE2_EXT_FITS_NONE         = 0x00400000, ///< ALL of the model fits failed
+    
+    
 } pmSourceMode2;
 

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

@@ -25 +25 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -42 +42 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 
@@ -65 +66 @@
 void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
 
+bool pmSourceMomentsHighOrder    (pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal);
+bool pmSourceMomentsRadialMoment (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal);
+bool pmSourceMomentsKronFluxes   (pmSource *source, float sigma,  float minSN, psImageMaskType maskVal);
+
 // if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
-
 bool pmSourceMoments(pmSource *source, float radius, float sigma, float minSN, float minKronRadius, psImageMaskType maskVal)
 {
@@ -74 +78 @@
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
-    // this function assumes the sky has been well-subtracted for the image
-    float sky = 0.0;
-
     if (source->moments == NULL) {
       source->moments = pmMomentsAlloc();
     }
-
-    float Sum = 0.0;
-    float Var = 0.0;
-    float SumCore = 0.0;
-    float VarCore = 0.0;
-    float R2 = PS_SQR(radius);
-    float minSN2 = PS_SQR(minSN);
-    float rsigma2 = 0.5 / PS_SQR(sigma);
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -110 +103 @@
     // of any object drops pretty quickly outside 1-2 sigmas.  (The exception is bright
     // saturated stars, for which we need to use a very large radius here)
+    // NOTE: if (mode & EXTERNAL) or (mode2 & MATCHED), do not re-calculate the centroid (use peak as centroid)
+    // (we still call this function because it sets moments->Sum,SN,Peak,nPixels
     if (!pmSourceMomentsGetCentroid (source, 1.5*sigma, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
         return false;
     }
 
+    pmSourceMomentsHighOrder (source, radius, sigma, minSN, maskVal);
+
+    // now calculate the 1st radial moment (for kron flux) using symmetrical averaging
+    pmSourceMomentsRadialMoment (source, radius, minKronRadius, maskVal);
+
+    // now calculate the kron flux values using the 1st radial moment
+    pmSourceMomentsKronFluxes (source, sigma, minSN, maskVal);
+
+    psTrace ("psModules.objects", 4, "Mrf: %f  KronFlux: %f  Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
+             source->moments->Mrf,   source->moments->KronFlux, 
+             source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
+             source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
+             source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
+
+    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  Npix: %d\n",
+             source->peak->xf, source->peak->yf, 
+             source->peak->rawFlux, sqrt(source->peak->detValue), 
+             source->moments->Mx, source->moments->My, 
+             source->moments->Sum, 
+             source->moments->Mxx, source->moments->Mxy, source->moments->Myy, 
+             source->moments->nPixels);
+
+    return(true);
+}
+
+bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
+
+    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
+    // Sum = SUM (z - sky)
+    // X1  = SUM (x - xc)*(z - sky)
+    // .. etc
+
+    float sky = 0.0;
+
+    float peakPixel = -PS_MAX_F32;
+    psS32 numPixels = 0;
+    float Sum = 0.0;
+    float Var = 0.0;
+    float X1 = 0.0;
+    float Y1 = 0.0;
+    float R2 = PS_SQR(radius);
+    float minSN2 = PS_SQR(minSN);
+    float rsigma2 = 0.5 / PS_SQR(sigma);
+
+    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
+    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
+
+    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
+    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
+    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
+    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
+    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
+
+    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
+    // not depend on the fractional pixel location of the source.  However, the aperture
+    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
+    // position of the expected centroid
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+        float yDiff = row + 0.5 - yPeak;
+        if (fabs(yDiff) > radius) continue;
+
+        float *vPix = source->pixels->data.F32[row];
+        float *vWgt = source->variance ? source->variance->data.F32[row] : source->pixels->data.F32[row];
+
+        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+        // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
+
+        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
+            if (vMsk) {
+                if (*vMsk & maskVal) {
+                    vMsk++;
+                    continue;
+                }
+                vMsk++;
+            }
+            if (isnan(*vPix)) continue;
+
+            float xDiff = col + 0.5 - xPeak;
+            if (fabs(xDiff) > radius) continue;
+
+            // radius is just a function of (xDiff, yDiff)
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            if (r2 > R2) continue;
+
+            float pDiff = *vPix - sky;
+            float wDiff = *vWgt;
+
+            // skip pixels below specified significance level.  for a PSFs, this
+            // over-weights the wings of bright stars compared to those of faint stars.
+            // for the estimator used for extended source analysis (where the window
+            // function is allowed to be arbitrarily large), we need to clip to avoid
+            // negative second moments.
+            if (PS_SQR(pDiff) < minSN2*wDiff) continue; // 
+            if ((minSN > 0.0) && (pDiff < 0)) continue; // 
+
+            // Apply a Gaussian window function.  Be careful with the window function.  S/N
+            // weighting over weights the sky for faint sources
+            if (sigma > 0.0) {
+                float z  = r2*rsigma2;
+                assert (z >= 0.0);
+                float weight  = exp(-z);
+
+                wDiff *= weight;
+                pDiff *= weight;
+            } 
+
+            Var += wDiff;
+            Sum += pDiff;
+
+            float xWght = xDiff * pDiff;
+            float yWght = yDiff * pDiff;
+
+            X1  += xWght;
+            Y1  += yWght;
+
+            peakPixel = PS_MAX (*vPix, peakPixel);
+            numPixels++;
+        }
+    }
+
+    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
+    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
+
+    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
+    if ((numPixels < minPixels) || (Sum <= 0)) {
+        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
+        return (false);
+    }
+
+    // calculate the first moment.
+    float Mx = X1/Sum;
+    float My = Y1/Sum;
+    if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
+        psTrace ("psModules.objects", 3, "extreme centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
+        return (false);
+    }
+    if ((fabs(Mx) > 2.0) || (fabs(My) > 2.0)) {
+        psTrace ("psModules.objects", 3, " big centroid swing; ok peak? %d, %d\n", source->peak->x, source->peak->y);
+    }
+
+    psTrace ("psModules.objects", 5, "id: %d, sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", source->id, sky, Mx, My, Sum, X1, Y1, numPixels);
+
+    // add back offset of peak in primary image
+    // also offset from pixel index to pixel coordinate
+    // (the calculation above uses pixel index instead of coordinate)
+    // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
+
+    // we only update the centroid if the position is not supplied from elsewhere
+    bool skipCentroid = false;
+    skipCentroid |= (source->mode  & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions
+    skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions
+
+    if (skipCentroid) {
+        source->moments->Mx = source->peak->xf;
+        source->moments->My = source->peak->yf;
+    } else {
+        source->moments->Mx = Mx + xGuess;
+        source->moments->My = My + yGuess;
+    }
+
+    source->moments->Sum = Sum;
+    source->moments->SN  = Sum / sqrt(Var);
+    source->moments->Peak = peakPixel;
+    source->moments->nPixels = numPixels;
+
+    return true;
+}
+
+float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) {
+
+    psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *src = sources->data[i]; // Source of interest
+        if (!src || !src->moments) {
+            continue;
+        }
+
+        if (src->mode & PM_SOURCE_MODE_BLEND) {
+            continue;
+        }
+
+        if (!src->moments->nPixels) continue;
+
+        if (src->moments->SN < PSF_SN_LIM) continue;
+
+        // XXX put in Mxx,Myy cut based on clump location
+
+        psVectorAppend(radii, src->moments->Mrf);
+    }
+
+    // find the peak in this image
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    if (!psVectorStats (stats, radii, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
+        psFree(stats);
+        return NAN;
+    }
+
+    float minRadius = stats->sampleMedian;
+
+    psFree(radii);
+    psFree(stats);
+    return minRadius;
+}
+
+bool pmSourceMomentsHighOrder (pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal) {
+
+    // this function assumes the sky has been well-subtracted for the image
+    float Sum = 0.0;
+    float R2 = PS_SQR(radius);
+    float minSN2 = PS_SQR(minSN);
+    float rsigma2 = 0.5 / PS_SQR(sigma);
+
     // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
-    // Xn  = SUM (x - xc)^n * (z - sky)
+    // Xn  = SUM (x - xc)^n * (z - sky) -- note that sky is 0.0 by definition here
     float XX = 0.0;
     float XY = 0.0;
@@ -129 +341 @@
     float YYYY = 0.0;
 
-    Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
-
-    // float dX = source->moments->Mx - source->peak->xf;
-    // float dY = source->moments->My - source->peak->yf;
-    // float dR = hypot(dX, dY);
-    // float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
-    // float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
     float Xo = source->moments->Mx;
     float Yo = source->moments->My;
@@ -154 +359 @@
 
         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
-        // psImageMaskType  *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
 
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -173 +377 @@
             if (r2 > R2) continue;
 
-            float fDiff = *vPix - sky;
+            float fDiff = *vPix;
             float pDiff = fDiff;
             float wDiff = *vWgt;
@@ -181 +385 @@
             // stars.
             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
-            if ((minSN > 0.0) && (pDiff < 0)) continue; // 
+            if ((minSN > 0.0) && (pDiff < 0)) continue; 
 
             // Apply a Gaussian window function.  Be careful with the window function.  S/N
-            // weighting over weights the sky for faint sources
+            // weighting over-weights the sky for faint sources
             if (sigma > 0.0) {
                 float z = r2 * rsigma2;
@@ -230 +434 @@
             XYYY  += xyyy;
             YYYY  += yyyy;
-
-            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
-            // XXX float r = sqrt(r2);
-            // XXX float rf = r * fDiff;
-            // XXX float rh = sqrt(r) * fDiff;
-            // XXX float rs = fDiff;
-            // XXX 
-            // XXX float rfw = r * pDiff;
-            // XXX float rhw = sqrt(r) * pDiff;
-            // XXX 
-            // XXX RF  += rf;
-            // XXX RH  += rh;
-            // XXX RS  += rs;
-            // XXX 
-            // XXX RFW  += rfw;
-            // XXX RHW  += rhw;
         }
     }
+    // NOT needed : source->moments->wSum = Sum;
+
     source->moments->Mxx = XX/Sum;
     source->moments->Mxy = XY/Sum;
@@ -263 +453 @@
     source->moments->Myyyy = YYYY/Sum;
 
-    // *** now calculate the 1st radial moment (for kron flux) -- symmetrical averaging
+    return true;
+}
+
+bool pmSourceMomentsRadialMoment (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) {
+
 
     float **vPix = source->pixels->data.F32;
-    float **vWgt = source->variance ? source->variance->data.F32 : source->pixels->data.F32;
     psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
 
@@ -272 +465 @@
     float RH = 0.0;
     float RS = 0.0;
+
+    // centroid around which to calculate the moments
+    float Xo = source->moments->Mx;
+    float Yo = source->moments->My;
+
+    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
+    // xCM, yCM from pixel coords to pixel index here.
+    float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
+    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
+
+    float R2 = PS_SQR(radius);
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -304 +508 @@
             if (r2 > R2) continue;
 
-            float fDiff1 = vPix[row][col] - sky;
-            float fDiff2 = vPix[yFlip][xFlip] - sky;
+            float fDiff1 = vPix[row][col];
+            float fDiff2 = vPix[yFlip][xFlip];
             float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
 
@@ -329 +533 @@
         kronRefRadius = MIN(radius, kronRefRadius);
     }
-    source->moments->Mrf = kronRefRadius;
-
-    // *** now calculate the kron flux values using the 1st radial moment
-
-    float radKinner = 1.0*kronRefRadius;
-    float radKron   = 2.5*kronRefRadius;
-    float radKouter = 4.0*kronRefRadius;
+
+    // if source is externally supplied and it already has a finite Mrf do not change it
+    if (! ((source->mode & PM_SOURCE_MODE_EXTERNAL) && isfinite(source->moments->Mrf))) {
+        source->moments->Mrf = kronRefRadius;
+    }
+
+    return true;
+}
+
+bool pmSourceMomentsKronFluxes (pmSource *source, float sigma, float minSN, psImageMaskType maskVal) {
+
+    float **vPix = source->pixels->data.F32;
+    float **vWgt = source->variance ? source->variance->data.F32 : source->pixels->data.F32;
+    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+
+    float radKinner = 1.0*source->moments->Mrf;
+    float radKron   = 2.5*source->moments->Mrf;
+    float radKouter = 4.0*source->moments->Mrf;
 
     int nKronPix = 0;
@@ -341 +556 @@
     int nInner = 0;
     int nOuter = 0;
-    Sum = Var = 0.0;
+    
+    float Sum = 0.0;
+    float Var = 0.0;
+    float SumCore = 0.0;
+    float VarCore = 0.0;
     float SumInner = 0.0;
     float SumOuter = 0.0;
+
+    // centroid around which to calculate the moments
+    float Xo = source->moments->Mx;
+    float Yo = source->moments->My;
+
+    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
+    // xCM, yCM from pixel coords to pixel index here.
+    float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
+    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
+
+    float minSN2 = PS_SQR(minSN);
 
     // calculate the Kron flux, and related fluxes (NO symmetrical averaging)
@@ -362 +592 @@
             float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
 
-            float fDiff1 = vPix[row][col] - sky;
+            float fDiff1 = vPix[row][col];
             float pDiff = fDiff1;
             float wDiff = vWgt[row][col];
@@ -376 +606 @@
                 Var += wDiff;
                 nKronPix ++;
-                // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
             }
 
@@ -397 +626 @@
     }
     // *** should I rescale these fluxes by pi R^2 / nNpix?
-    // XXX source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
-    // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
-    // XXX source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
-    // XXX source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
-    // XXX source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
-    // XXX source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
+    // XXX source->moments->KronCore    = SumCore       * M_PI *  PS_SQR(sigma) / nCorePix;
+    // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI *  PS_SQR(sigma) / nCorePix;
+    // XXX source->moments->KronFlux    = Sum           * M_PI *  PS_SQR(radKron) / nKronPix;
+    // XXX source->moments->KronFluxErr = sqrt(Var)     * M_PI *  PS_SQR(radKron) / nKronPix;
+    // XXX source->moments->KronFinner  = SumInner      * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
+    // XXX source->moments->KronFouter  = SumOuter      * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
 
     source->moments->KronCore    = SumCore;
@@ -408 +637 @@
     source->moments->KronFlux    = Sum;
     source->moments->KronFluxErr = sqrt(Var);
-    source->moments->KronFinner = SumInner;
-    source->moments->KronFouter = SumOuter;
+    source->moments->KronFinner  = SumInner;
+    source->moments->KronFouter  = SumOuter;
 
     // XXX not sure I should save this here...
@@ -416 +645 @@
     source->moments->KronRadiusPSF  = source->moments->Mrf;
 
-    psTrace ("psModules.objects", 4, "Mrf: %f  KronFlux: %f  Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
-             source->moments->Mrf,   source->moments->KronFlux, 
-             source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
-             source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
-             source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
-
-    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  sky: %f  Npix: %d\n",
-             source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, source->moments->nPixels);
-
-    return(true);
-}
-
-bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
-
-    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
-    // Sum = SUM (z - sky)
-    // X1  = SUM (x - xc)*(z - sky)
-    // .. etc
-
-    float sky = 0.0;
-
-    float peakPixel = -PS_MAX_F32;
-    psS32 numPixels = 0;
-    float Sum = 0.0;
-    float Var = 0.0;
-    float X1 = 0.0;
-    float Y1 = 0.0;
-    float R2 = PS_SQR(radius);
-    float minSN2 = PS_SQR(minSN);
-    float rsigma2 = 0.5 / PS_SQR(sigma);
-
-    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
-    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
-
-    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
-    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
-    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
-    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
-    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
-
-    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
-    // not depend on the fractional pixel location of the source.  However, the aperture
-    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
-    // position of the expected centroid
-
-    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-
-        float yDiff = row + 0.5 - yPeak;
-        if (fabs(yDiff) > radius) continue;
-
-        float *vPix = source->pixels->data.F32[row];
-        float *vWgt = source->variance ? source->variance->data.F32[row] : source->pixels->data.F32[row];
-
-        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
-        // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
-
-        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
-            if (vMsk) {
-                if (*vMsk & maskVal) {
-                    vMsk++;
-                    continue;
-                }
-                vMsk++;
-            }
-            if (isnan(*vPix)) continue;
-
-            float xDiff = col + 0.5 - xPeak;
-            if (fabs(xDiff) > radius) continue;
-
-            // radius is just a function of (xDiff, yDiff)
-            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-            if (r2 > R2) continue;
-
-            float pDiff = *vPix - sky;
-            float wDiff = *vWgt;
-
-            // skip pixels below specified significance level.  for a PSFs, this
-            // over-weights the wings of bright stars compared to those of faint stars.
-            // for the estimator used for extended source analysis (where the window
-            // function is allowed to be arbitrarily large), we need to clip to avoid
-            // negative second moments.
-            if (PS_SQR(pDiff) < minSN2*wDiff) continue; // 
-            if ((minSN > 0.0) && (pDiff < 0)) continue; // 
-
-            // Apply a Gaussian window function.  Be careful with the window function.  S/N
-            // weighting over weights the sky for faint sources
-            if (sigma > 0.0) {
-                float z  = r2*rsigma2;
-                assert (z >= 0.0);
-                float weight  = exp(-z);
-
-                wDiff *= weight;
-                pDiff *= weight;
-            } 
-
-            Var += wDiff;
-            Sum += pDiff;
-
-            float xWght = xDiff * pDiff;
-            float yWght = yDiff * pDiff;
-
-            X1  += xWght;
-            Y1  += yWght;
-
-            peakPixel = PS_MAX (*vPix, peakPixel);
-            numPixels++;
-        }
-    }
-
-    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
-    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
-
-    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
-    if ((numPixels < minPixels) || (Sum <= 0)) {
-        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
-        return (false);
-    }
-
-    // calculate the first moment.
-    float Mx = X1/Sum;
-    float My = Y1/Sum;
-    if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
-        psTrace ("psModules.objects", 3, "extreme centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
-        return (false);
-    }
-    if ((fabs(Mx) > 2.0) || (fabs(My) > 2.0)) {
-        psTrace ("psModules.objects", 3, " big centroid swing; ok peak? %d, %d\n", source->peak->x, source->peak->y);
-    }
-
-    psTrace ("psModules.objects", 5, "id: %d, sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", source->id, sky, Mx, My, Sum, X1, Y1, numPixels);
-
-    // add back offset of peak in primary image
-    // also offset from pixel index to pixel coordinate
-    // (the calculation above uses pixel index instead of coordinate)
-    // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
-
-    // we only update the centroid if the position is not supplied from elsewhere
-    bool skipCentroid = false;
-    skipCentroid |= (source->mode  & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions
-    skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions
-
-    if (skipCentroid) {
-        source->moments->Mx = source->peak->xf;
-        source->moments->My = source->peak->yf;
-    } else {
-        source->moments->Mx = Mx + xGuess;
-        source->moments->My = My + yGuess;
-    }
-
-    source->moments->Sum = Sum;
-    source->moments->SN  = Sum / sqrt(Var);
-    source->moments->Peak = peakPixel;
-    source->moments->nPixels = numPixels;
-
     return true;
 }
-
-float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) {
-
-    psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32);
-
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *src = sources->data[i]; // Source of interest
-        if (!src || !src->moments) {
-            continue;
-        }
-
-        if (src->mode & PM_SOURCE_MODE_BLEND) {
-            continue;
-        }
-
-        if (!src->moments->nPixels) continue;
-
-        if (src->moments->SN < PSF_SN_LIM) continue;
-
-        // XXX put in Mxx,Myy cut based on clump location
-
-        psVectorAppend(radii, src->moments->Mrf);
-    }
-
-    // find the peak in this image
-    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-
-    if (!psVectorStats (stats, radii, NULL, NULL, 0)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
-        psFree(stats);
-        return NAN;
-    }
-
-    float minRadius = stats->sampleMedian;
-
-    psFree(radii);
-    psFree(stats);
-    return minRadius;
-}
-

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

@@ -33 +33 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -40 +40 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"
@@ -113 +114 @@
     source->apFluxErr = NAN;
 
+    pmModelStatus badModel = PM_MODEL_STATUS_NONE;
+    badModel |= PM_MODEL_STATUS_BADARGS;
+    badModel |= PM_MODEL_STATUS_OFFIMAGE;
+    badModel |= PM_MODEL_STATUS_NAN_CHISQ;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID;
+    badModel |= PM_MODEL_PCM_FAIL_GUESS;
+
     // XXXXXX review:
     // Select the 'best' model -- this is used for PSF_QF,_PERFECT & ???. isPSF is true if this
@@ -162 +171 @@
         for (int i = 0; i < source->modelFits->n; i++) {
             pmModel *model = source->modelFits->data[i];
-            if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
+            if (model->flags & badModel) continue;
             status = pmSourcePhotometryModel (&model->mag, NULL, model);
             if (model == source->modelEXT) foundEXT = true;
@@ -902 +911 @@
 }
 
+bool pmSourceChisqModelFlux (pmSource *source, pmModel *model, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(model, false);
+
+    float dC = 0.0;
+    int Npix = 0;
+
+    psVector *params = model->params;
+    psImage  *image = source->pixels;
+    psImage  *modelFlux = source->modelFlux;
+    psImage  *mask = source->maskObj;
+    psImage  *variance = source->variance;
+
+    float Io = params->data.F32[PM_PAR_I0];
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+        for (int ix = 0; ix < image->numCols; ix++) {
+
+            // skip pixels which are masked
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue;
+
+            if (variance->data.F32[iy][ix] <= 0) continue;
+
+            dC += PS_SQR (image->data.F32[iy][ix] - Io*modelFlux->data.F32[iy][ix]) / variance->data.F32[iy][ix];
+            Npix ++;
+        }
+    }
+    model->nPix = Npix;
+    model->nDOF = Npix - model->nPar;
+    model->chisq = dC;
+    model->chisqNorm = dC / model->nDOF;
+
+    return (true);
+}
+
 double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
 {

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

@@ -42 +42 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -45 +45 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -46 +46 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

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

@@ -41 +41 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 

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

@@ -41 +41 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 

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

@@ -23 +23 @@
 #include "pmSourceDiffStats.h"
 #include "pmSourceSatstar.h"
+#include "pmSourceLensing.h"
 #include "pmSource.h"
 #include "pmSourceFitModel.h"

@@ -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