Changeset 32347

Timestamp:

Sep 6, 2011, 1:02:53 PM (15 years ago)

Author:

eugene

Message:

• add concept of saddlePoints to peaks (not actually used in the end)
• add tmp flags to mark sources for analysis or not in psphotStack
• autocode the pmSourceIO_CMF_PS1_* functions
• use 1D gauss approx for convolution in PCM fitting
• added pmSourceExtFitPars (not actually used)
• in model guess, use 1st radial moments to define size (if it exists)
• include PSF_INST_MAG, AP_MAG, KRON_MAG in xfit output
• fix the position for extended source fits (avoid instability)
• Sersic-like models (incl. Exp and Dev) use Reff, not sigma; conversion tools need to respect this
• only use a single pass on the centroid (unwindowed, but limited to 1.5 sigma radius) - this avoids moving the centroid because of nearby neighbors
• use symmetrical averaging (geometric mean) to calculated 1st radial moment (and avoid neighbor biases), do not use symm. averaging for the flux
• fix the integration of the sersic, pgauss, and related model functions.
• fix the central pixel to have the full flux for sersic-like models (interpolated value)

Location:

trunk/psModules/src

Files:

• camera/pmReadoutFake.c (modified) (2 diffs)
• objects (modified) (1 prop)
• objects/Makefile.am (modified) (1 diff)
• objects/mksource.pl (modified) (2 diffs)
• objects/models/pmModel_DEV.c (modified) (7 diffs)
• objects/models/pmModel_EXP.c (modified) (6 diffs)
• objects/models/pmModel_PGAUSS.c (modified) (2 diffs)
• objects/models/pmModel_PS1_V1.c (modified) (3 diffs)
• objects/models/pmModel_QGAUSS.c (modified) (2 diffs)
• objects/models/pmModel_RGAUSS.c (modified) (2 diffs)
• objects/models/pmModel_SERSIC.CP.h (copied) (copied from branches/eam_branches/ipp-20110710/psModules/src/objects/models/pmModel_SERSIC.CP.h )
• objects/models/pmModel_SERSIC.c (modified) (7 diffs)
• objects/pmFootprintCullPeaks.c (modified) (1 diff)
• objects/pmModelUtils.c (modified) (1 diff)
• objects/pmMoments.c (modified) (2 diffs)
• objects/pmMoments.h (modified) (1 diff)
• objects/pmPCM_MinimizeChisq.c (modified) (6 diffs)
• objects/pmPCMdata.c (modified) (4 diffs)
• objects/pmPCMdata.h (modified) (2 diffs)
• objects/pmPSF.c (modified) (7 diffs)
• objects/pmPSF.h (modified) (1 diff)
• objects/pmPeaks.c (modified) (2 diffs)
• objects/pmPeaks.h (modified) (1 diff)
• objects/pmSource.c (modified) (3 diffs)
• objects/pmSource.h (modified) (2 diffs)
• objects/pmSourceExtendedPars.c (modified) (1 diff)
• objects/pmSourceExtendedPars.h (modified) (2 diffs)
• objects/pmSourceFitModel.c (modified) (2 diffs)
• objects/pmSourceFitPCM.c (modified) (6 diffs)
• objects/pmSourceIO_CMF.c.in (modified) (11 diffs)
• objects/pmSourceIO_CMF_PS1_DV1.c (modified) (2 diffs)
• objects/pmSourceIO_CMF_PS1_DV2.c (modified) (2 diffs)
• objects/pmSourceIO_CMF_PS1_SV1.c (modified) (2 diffs)
• objects/pmSourceIO_CMF_PS1_V1.c (deleted)
• objects/pmSourceIO_CMF_PS1_V2.c (deleted)
• objects/pmSourceIO_CMF_PS1_V3.c (deleted)
• objects/pmSourceIO_CMP.c (modified) (2 diffs)
• objects/pmSourceIO_OBJ.c (modified) (1 diff)
• objects/pmSourceIO_PS1_CAL_0.c (modified) (3 diffs)
• objects/pmSourceIO_PS1_DEV_0.c (modified) (2 diffs)
• objects/pmSourceIO_PS1_DEV_1.c (modified) (3 diffs)
• objects/pmSourceIO_SMPDATA.c (modified) (2 diffs)
• objects/pmSourceIO_SX.c (modified) (1 diff)
• objects/pmSourceMoments.c (modified) (10 diffs)
• objects/pmSourceOutputs.c (modified) (2 diffs)
• objects/pmSourceOutputs.h (modified) (1 diff)
• objects/pmSourcePhotometry.c (modified) (1 diff)
• objects/pmSourcePlotPSFModel.c (modified) (1 diff)

trunk/psModules/src/camera/pmReadoutFake.c

@@ -51 +51 @@
 
     psF32 *params = model->params->data.F32; // Model parameters
-    psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO); // Ellipse axes
+    psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO, model->type); // Ellipse axes
     // Curiously, the minor axis can be larger than the major axis, so need to check.
     if (axes.major >= axes.minor) {
@@ -58 +58 @@
         axes.major = axes.minor;
     }
-    return pmPSF_AxesToModel(params, axes);
+    return pmPSF_AxesToModel(params, axes, model->type);
 }
 

trunk/psModules/src/objects

@@ -5 +5 @@
 *.la
 *.lo
+pmSourceIO_CMF_PS1_V1.c
+pmSourceIO_CMF_PS1_V2.c
+pmSourceIO_CMF_PS1_V3.c

@@ -5 +5 @@
 *.la
 *.lo
+pmSourceIO_CMF_PS1_V1.c
+pmSourceIO_CMF_PS1_V2.c
+pmSourceIO_CMF_PS1_V3.c

trunk/psModules/src/objects/Makefile.am

@@ -130 +130 @@
 
 # pmSourceID_CMF_* functions use a common framework
-BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.v1.c pmSourceIO_CMF_PS1_V2.v1.c pmSourceIO_CMF_PS1_V3.v1.c
+BUILT_SOURCES = pmSourceIO_CMF_PS1_V1.c pmSourceIO_CMF_PS1_V2.c pmSourceIO_CMF_PS1_V3.c
 
-pmSourceIO_CMF_PS1_V1.v1.c : pmSourceIO_CMF.c.in mksource.pl
-        mksource.pl pmSourceIO_CMF.c.in PS1_V1 pmSourceIO_CMF_PS1_V1.v1.c
+pmSourceIO_CMF_PS1_V1.c : pmSourceIO_CMF.c.in mksource.pl
+        mksource.pl pmSourceIO_CMF.c.in PS1_V1 pmSourceIO_CMF_PS1_V1.c
 
-pmSourceIO_CMF_PS1_V2.v1.c : pmSourceIO_CMF.c.in mksource.pl
-        mksource.pl pmSourceIO_CMF.c.in PS1_V2 pmSourceIO_CMF_PS1_V2.v1.c
+pmSourceIO_CMF_PS1_V2.c : pmSourceIO_CMF.c.in mksource.pl
+        mksource.pl pmSourceIO_CMF.c.in PS1_V2 pmSourceIO_CMF_PS1_V2.c
 
-pmSourceIO_CMF_PS1_V3.v1.c : pmSourceIO_CMF.c.in mksource.pl
-        mksource.pl pmSourceIO_CMF.c.in PS1_V2 pmSourceIO_CMF_PS1_V3.v1.c
+pmSourceIO_CMF_PS1_V3.c : pmSourceIO_CMF.c.in mksource.pl
+        mksource.pl pmSourceIO_CMF.c.in PS1_V3 pmSourceIO_CMF_PS1_V3.c
 
 # EXTRA_DIST = pmErrorCodes.h.in pmErrorCodes.dat pmErrorCodes.c.in

trunk/psModules/src/objects/mksource.pl

@@ -14 +14 @@
 
 # see if we can add in PS1_DV* and PS1_SV* as well...
-@cmfmodes = ("PS1_V1", 1,
+%cmfmodes = ("PS1_V1", 1,
              "PS1_V2", 2,
              "PS1_V3", 3);
+
+print "1: $cmfmodes{1}\n";
+print "PS1_V1: $cmfmodes{'PS1_V1'}\n";
 
 open (FILE, "$template") || die "failed to open template $template\n";
@@ -50 +53 @@
 
     if ($gtMode) {
+        # print "gtMode : $line\n";
         $thisLevel = $cmfmodes{$gtMode};
         $myLevel = $cmfmodes{$cmfmode};
-        if ($thisLevel <= $myLevel) { next; }
-        $line =~ s|\@<\S*\@\s*||;
+        print "gtMode : $gtMode vs $cmfmode, $thisLevel, $myLevel\n";
+        if ($myLevel <= $thisLevel) { next; }
+        $line =~ s|\@>\S*\@\s*||;
     }
 
     if ($ltMode) {
+        # print "ltMode : $line\n";
         $thisLevel = $cmfmodes{$ltMode};
         $myLevel = $cmfmodes{$cmfmode};
-        if ($thisLevel >= $myLevel) { next; }
-        $line =~ s|\@>\S*\@\s*||;
+        print "ltMode : $ltMode vs $cmfmode, $thisLevel, $myLevel\n";
+        if ($myLevel >= $thisLevel) { next; }
+        $line =~ s|\@<\S*\@\s*||;
     }
 

trunk/psModules/src/objects/models/pmModel_DEV.c

@@ -89 +89 @@
 static bool limitsApply = true;         // Apply limits?
 
+# include "pmModel_SERSIC.CP.h"
+
 psF32 PM_MODEL_FUNC (psVector *deriv,
                      const psVector *params,
@@ -94 +96 @@
 {
     psF32 *PAR = params->data.F32;
-
-    float index = 0.5 / ALPHA;
-    float bn = 1.9992*index - 0.3271;
-    float Io = exp(bn);
 
     psF32 X  = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS];
@@ -105 +103 @@
     psF32 z  = (PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y);
 
-    assert (z >= 0);
+    // If the elliptical contour is defined in a valid way, we should never trigger this
+    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
+    // NaN values for negative values of z
+    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)
+        psEllipseShape shape;
+
+        shape.sx  = PAR[PM_PAR_SXX];
+        shape.sy  = PAR[PM_PAR_SYY];
+        shape.sxy = PAR[PM_PAR_SXY];
+
+        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
+        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+
+        // 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;
@@ -120 +199 @@
         psF32 *dPAR = deriv->data.F32;
 
+        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);
-
-        dPAR[PM_PAR_SKY]  = +1.0;
-        dPAR[PM_PAR_I0]   = +2.0*f1;
 
         assert (isfinite(z1));
@@ -223 +302 @@
 
     // set the shape parameters
+    // XXX adjust this?
     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
       return false;
@@ -246 +326 @@
 }
 
+// A DeVaucouleur model is equivalent to a Sersic with index = 4.0
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    float z, norm;
     psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
+    // for a non-circular DeVaucouleur, the flux of the Rmajor equivalent is scaled by the AspectRatio
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-
-    # define DZ 0.25
-
-    float f0 = 1.0;
-    float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = exp(-pow(z,ALPHA));
-        z += DZ;
-        f2 = exp(-pow(z,ALPHA));
-        norm += f0 + 4*f1 + f2;
-        f0 = f2;
-    }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    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);
@@ -297 +369 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 

trunk/psModules/src/objects/models/pmModel_EXP.c

@@ -81 +81 @@
 static bool limitsApply = true;         // Apply limits?
 
+# include "pmModel_SERSIC.CP.h"
+
 psF32 PM_MODEL_FUNC (psVector *deriv,
                      const psVector *params,
@@ -93 +95 @@
     psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
 
-    // XXX if the elliptical contour is defined in valid way, this step should not be required.
-    // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
-    // for negative values of z
-    // XXX use an assert here to force the elliptical parameters to be correctly determined
-    // if (z < 0) z = 0;
-    assert (z >= 0);
-
-    psF32 f2 = sqrt(z);
-    psF32 f1 = exp(-f2);
+    // If the elliptical contour is defined in a valid way, we should never trigger this
+    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
+    // NaN values for negative values of z
+    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);
+
+    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)
+        psEllipseShape shape;
+
+        shape.sx  = PAR[PM_PAR_SXX];
+        shape.sy  = PAR[PM_PAR_SYY];
+        shape.sxy = PAR[PM_PAR_SXY];
+
+        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
+        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+
+        // 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;
@@ -118 +196 @@
         // 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*z0/sqrt(0.01) : 0.5*z0/sqrt(z);
+        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?
@@ -216 +294 @@
 
     // set the shape parameters
+    // XXX adjust this?
     if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
       return false;
@@ -233 +312 @@
 }
 
+// An exponential model is equivalent to a Sersic with index = 1.0
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    float z, norm;
     psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
+    // for a non-circular Exponential, the flux of the Rmajor equivalent is scaled by the AspectRatio
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z) = exp(-sqrt(z)) [0 to infinity]
-    norm = 0.0;
-
-    # define DZ 0.25
-
-    float f0 = 1.0;
-    float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = exp(-sqrt(z));
-        z += DZ;
-        f2 = exp(-sqrt(z));
-        norm += f0 + 4*f1 + f2;
-        f0 = f2;
-    }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    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);
@@ -284 +355 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 

trunk/psModules/src/objects/models/pmModel_PGAUSS.c

@@ -217 +217 @@
 }
 
-psF64 PM_MODEL_FLUX(const psVector *params)
-{
-    float norm, z;
+// integrate the model to get the full flux
+psF64 PM_MODEL_FLUX (const psVector *params)
+{
+    float z, norm;
     psEllipseShape shape;
 
@@ -228 +229 @@
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
+    float AspectRatio = axes.minor / axes.major;
+
+    // flux = 2 \pi \int f(r) r dr
     norm = 0.0;
 
-    # define DZ 0.25
-
-    float f0 = 1.0;
+    # define DR 0.25
+
+    // f = f(r) * r
+    float f0 = 0.0;
     float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
-        z += DZ;
-        f2 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
+    for (float r = DR; r < 150; r += DR) {
+        z = 0.5 * PS_SQR(r / axes.major);
+        f1 = r / (1 + z + z*z/2.0 + z*z*z/6.0);
+        r += DR;
+        z = 0.5 * PS_SQR(r / axes.major);
+        f2 = r / (1 + z + z*z/2.0 + z*z*z/6.0);
         norm += f0 + 4*f1 + f2;
         f0 = f2;
     }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    norm *= DR / 3.0;
+
+    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
 
     return(Flux);

trunk/psModules/src/objects/models/pmModel_PS1_V1.c

@@ -14 +14 @@
    * PM_PAR_XPOS 2  - X center of object
    * PM_PAR_YPOS 3  - Y center of object
-   * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
-   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
+   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
+   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    * PM_PAR_SXY 6   - X*Y term of elliptical contour
    * PM_PAR_7   7   - amplitude of the linear component (k)
@@ -239 +239 @@
 }
 
+// integrate the model to get the full flux
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
@@ -250 +251 @@
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
+    float AspectRatio = axes.minor / axes.major;
+
+    // flux = 2 \pi \int f(r) r dr
     norm = 0.0;
 
-    # define DZ 0.25
-
-    float f0 = 1.0;
+    # define DR 0.25
+
+    // f = f(r) * r
+    float f0 = 0.0;
     float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
-        z += DZ;
-        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+    for (float r = DR; r < 150; r += DR) {
+        z = 0.5 * PS_SQR(r / axes.major);
+        f1 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+        r += DR;
+        z = 0.5 * PS_SQR(r / axes.major);
+        f2 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         norm += f0 + 4*f1 + f2;
         f0 = f2;
     }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    norm *= DR / 3.0;
+
+    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
 
     return(Flux);

trunk/psModules/src/objects/models/pmModel_QGAUSS.c

@@ -240 +240 @@
 }
 
+// integrate the model to get the full flux
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
@@ -251 +252 @@
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
+    float AspectRatio = axes.minor / axes.major;
+
+    // flux = 2 \pi \int f(r) r dr
     norm = 0.0;
 
-    # define DZ 0.25
-
-    float f0 = 1.0;
+    # define DR 0.25
+
+    // f = f(r) * r
+    float f0 = 0.0;
     float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
-        z += DZ;
-        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+    for (float r = DR; r < 150; r += DR) {
+        z = 0.5 * PS_SQR(r / axes.major);
+        f1 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+        r += DR;
+        z = 0.5 * PS_SQR(r / axes.major);
+        f2 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         norm += f0 + 4*f1 + f2;
         f0 = f2;
     }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    norm *= DR / 3.0;
+
+    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
 
     return(Flux);

trunk/psModules/src/objects/models/pmModel_RGAUSS.c

@@ -229 +229 @@
 }
 
+// integrate the model to get the full flux
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    float norm, z;
+    float z, norm;
     psEllipseShape shape;
 
@@ -240 +241 @@
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
+    float AspectRatio = axes.minor / axes.major;
+
+    // flux = 2 \pi \int f(r) r dr
     norm = 0.0;
 
-    # define DZ 0.25
-
-    float f0 = 1.0;
+    # define DR 0.25
+
+    // f = f(r) * r
+    float f0 = 0.0;
     float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        f1 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
-        z += DZ;
-        f2 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
+    for (float r = DR; r < 150; r += DR) {
+        z = 0.5 * PS_SQR(r / axes.major);
+        f1 = r / (1 + z + pow(z, PAR[PM_PAR_7]));
+        r += DR;
+        z = 0.5 * PS_SQR(r / axes.major);
+        f2 = r / (1 + z + pow(z, PAR[PM_PAR_7]));
         norm += f0 + 4*f1 + f2;
         f0 = f2;
     }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    norm *= DR / 3.0;
+
+    psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio;
 
     return(Flux);

trunk/psModules/src/objects/models/pmModel_SERSIC.c

@@ -13 +13 @@
    * PM_PAR_XPOS 2  - X center of object
    * PM_PAR_YPOS 3  - Y center of object
-   * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
-   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
+   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
+   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    * PM_PAR_SXY 6   - X*Y term of elliptical contour
    * PM_PAR_7   7   - normalized sersic parameter
+
+   * 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:
+     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)
 
    note that a standard sersic model uses exp(-K*(z^(1/2n) - 1).  the additional elements (K,
@@ -85 +91 @@
 static bool limitsApply = true;         // Apply limits?
 
+# include "pmModel_SERSIC.CP.h"
+
 psF32 PM_MODEL_FUNC (psVector *deriv,
                      const psVector *params,
@@ -97 +105 @@
     psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
 
-    // XXX if the elliptical contour is defined in valid way, this step should not be required.
-    // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values
-    // for negative values of z
-    // XXX use an assert here to force the elliptical parameters to be correctly determined
-    // if (z < 0) z = 0;
-    assert (z >= 0);
+    // If the elliptical contour is defined in a valid way, we should never trigger this
+    // assert.  Other models (like PGAUSS) don't use fractional powers, and thus do not have
+    // NaN values for negative values of z
+    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,PAR[PM_PAR_7]);
+    psF32 f2 = bn*pow(z,par7);
     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)
+        psEllipseShape shape;
+
+        shape.sx  = PAR[PM_PAR_SXX];
+        shape.sy  = PAR[PM_PAR_SYY];
+        shape.sxy = PAR[PM_PAR_SXY];
+
+        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
+        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+
+        // 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;
@@ -121 +201 @@
         psF32 *dPAR = deriv->data.F32;
 
-        // gradient is infinite for z = 0; saturate at z = 0.01
-        psF32 z1 = (z < 0.01) ? z0*bn*PAR[PM_PAR_7]*pow(0.01,PAR[PM_PAR_7] - 1.0) : z0*bn*PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0);
-
         dPAR[PM_PAR_SKY]  = +1.0;
         dPAR[PM_PAR_I0]   = +f1;
-        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z);
+
+        // 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;
 
@@ -269 +350 @@
     float Io = exp(0.5*bn);
 
-    // XXX do we need this factor of sqrt(2)?
-    // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx);
-    // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy);
-
     float Sxx = PS_MAX(0.5, shape.sx);
     float Syy = PS_MAX(0.5, shape.sy);
@@ -294 +371 @@
 }
 
+// A sersic model has an integral flux which can be analytically represented
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    float z, norm;
     psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 
-    // Area is equivalent to 2 pi sigma^2
+    // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-
-    # define DZ 0.25
-
-    float f0 = 1.0;
-    float f1, f2;
-    for (z = DZ; z < 150; z += DZ) {
-        // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
-        f1 = exp(-pow(z,PAR[PM_PAR_7]));
-        z += DZ;
-        f2 = exp(-pow(z,PAR[PM_PAR_7]));
-        norm += f0 + 4*f1 + f2;
-        f0 = f2;
-    }
-    norm *= DZ / 3.0;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
+    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);
@@ -346 +414 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
+    shape.sx  = PAR[PM_PAR_SXX];
+    shape.sy  = PAR[PM_PAR_SYY];
     shape.sxy = PAR[PM_PAR_SXY];
 

trunk/psModules/src/objects/pmFootprintCullPeaks.c

@@ -181 +181 @@
             if (!myFP->n) {
                 psWarning ("empty footprint?");
+                psFree (myFP);
                 continue;
             }

trunk/psModules/src/objects/pmModelUtils.c

@@ -122 +122 @@
     if (!isfinite(axes.theta)) return false;
 
+    // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width 
+    float scale = (isfinite(moments->Mrf) && (moments->Mrf > 0.0)) ? moments->Mrf / axes.major : 1.0;
+    axes.major *= scale;
+    axes.minor *= scale;
+
     psEllipseShape shape = psEllipseAxesToShape (axes);
 

trunk/psModules/src/objects/pmMoments.c

@@ -32 +32 @@
     tmp->Mrf = 0.0;
     tmp->Mrh = 0.0;
+
+    tmp->KronCore = 0.0;
+    tmp->KronCoreErr = 0.0;
+
     tmp->KronFlux = 0.0;
     tmp->KronFluxErr = 0.0;
@@ -37 +41 @@
     tmp->KronFinner = 0.0;
     tmp->KronFouter = 0.0;
+
+    tmp->KronFluxPSF = 0.0;
+    tmp->KronFluxPSFErr = 0.0;
+    tmp->KronRadiusPSF = 0.0;
 
     tmp->Mx = 0.0;

trunk/psModules/src/objects/pmMoments.h

@@ -51 +51 @@
     int nPixels;  ///< Number of pixels used.
 
+    float KronFluxPSF; ///< Kron Flux using PSF-optimized window
+    float KronFluxPSFErr; ///< Kron Flux Error using PSF-optimized window
+    float KronRadiusPSF; ///< Kron Radius using PSF-optimized window (Flux in 2.5 Radius)
+
     float KronCore;    ///< flux in r < 1.0*Mrf
     float KronCoreErr;    ///< error on flux in r < 1.0*Mrf

trunk/psModules/src/objects/pmPCM_MinimizeChisq.c

@@ -45 +45 @@
 # define USE_FFT 1
 # define PRE_CONVOLVE 1
+# define TESTCOPY 0
 
 bool pmPCM_MinimizeChisq (
@@ -93 +94 @@
         psFree (pcm->psfFFT);
     }
-    pcm->psfFFT = psImageConvolveKernelInit(pcm->modelFlux, pcm->psf);
+# if (!TESTCOPY)
+    if (!pcm->use1Dgauss) {
+        pcm->psfFFT = psImageConvolveKernelInit(pcm->modelFlux, pcm->psf);
+    }
+# endif
 # endif    
 
@@ -285 +290 @@
 
             // Convert i/j to image space:
-            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
-            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+            coord->data.F32[0] = (psF32) (j + 0.5 + source->pixels->col0);
+            coord->data.F32[1] = (psF32) (i + 0.5 + source->pixels->row0);
 
             pcm->modelFlux->data.F32[i][j] = pcm->modelConv->modelFunc (deriv, params, coord);
@@ -309 +314 @@
 # if (PRE_CONVOLVE)
     // convolve model image and derivative images with pre-convolved kernel
-    psImageConvolveKernel (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psfFFT);
+
+// XXX for a test, just copy, rather than convolve
+# if (TESTCOPY)
+    psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
+# 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 (pcm->modelConvFlux, pcm->sigma, pcm->nsigma);
+    } else {
+        psImageConvolveKernel (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psfFFT);
+    }
+# endif
+
     for (int n = 0; n < pcm->dmodelsFlux->n; n++) {
         if (pcm->dmodelsFlux->data[n] == NULL) continue;
@@ -315 +335 @@
         psImage *dmodel = pcm->dmodelsFlux->data[n];
         psImage *dmodelConv = pcm->dmodelsConvFlux->data[n];
-        psImageConvolveKernel (dmodelConv, dmodel, NULL, 0, pcm->psfFFT);
+# if (TESTCOPY)
+        psImageCopy (dmodelConv, dmodel, dmodel->type.type);
+# 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 (dmodelConv, pcm->sigma, pcm->nsigma);
+        } else {
+            psImageConvolveKernel (dmodelConv, dmodel, NULL, 0, pcm->psfFFT);
+        }
+# endif
     }
 # else
@@ -325 +357 @@
         psImage *dmodel = pcm->dmodelsFlux->data[n];
         psImage *dmodelConv = pcm->dmodelsConvFlux->data[n];
-        psImageConvolveFFT (dmodelConv, dmodel, NULL, 0, pcm->psf);
+
+        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 (dmodelConv, pcm->sigma, pcm->nsigma);
+        } else {
+            psImageConvolveFFT (dmodelConv, dmodel, NULL, 0, pcm->psf);
+        }
     }
 # endif // PRE-CONVOLVE

trunk/psModules/src/objects/pmPCMdata.c

@@ -41 +41 @@
 #include "pmPCMdata.h"
 
+# define USE_DELTA_PSF 0
+# define USE_1D_GAUSS 1
+
 static void pmPCMdataFree (pmPCMdata *pcm) {
 
@@ -88 +91 @@
     pcm->constraint = NULL;
     pcm->nDOF = 0;
+
+    // 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->sigma = 1.0; // this should be set to something sensible when the psf is known
 
     return pcm;
@@ -257 +266 @@
     pcm->nDOF = nPix - nParams;
 
+# if (USE_1D_GAUSS)
+    pmModel *modelPSF = source->modelPSF;
+    psAssert (modelPSF, "psf model must be defined");
+    
+    psEllipseShape shape;
+    psEllipseAxes axes;
+
+    shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
+    shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
+    shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
+    axes = psEllipseShapeToAxes (shape, 20.0);
+    
+    float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[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;
+# endif
+
     return pcm;
 }
@@ -368 +397 @@
     return true;
 }
+
+// construct a realization of the source model
+bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize) {
+
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    // select appropriate model
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+
+    // 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);
+
+    // 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");
+    
+        psEllipseShape shape;
+        psEllipseAxes axes;
+
+        shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
+        shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
+        shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
+        axes = psEllipseShapeToAxes (shape, 20.0);
+    
+        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[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);
+    } 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;
+}
+

trunk/psModules/src/objects/pmPCMdata.h

@@ -35 +35 @@
     int nPar;
     int nDOF;
+
+    bool use1Dgauss;
+    float sigma;
+    float nsigma;
 } pmPCMdata;
 
@@ -90 +94 @@
 bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
 
+bool pmPCMCacheModel (pmSource *source, psImageMaskType maskVal, int psfSize);
 
 /// @}

trunk/psModules/src/objects/pmPSF.c

@@ -281 +281 @@
 // convert the parameters used in the fitted source model
 // to the parameters used in the 2D PSF model
+// XXX this function may be invalid for SERSIC, DEV, EXP models (SQRT2 not used?)
 bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis)
 {
@@ -307 +308 @@
 // convert the PSF parameters used in the 2D PSF model fit into the
 // parameters used in the source model
+// XXX this function may be invalid for SERSIC, DEV, EXP models (SQRT2 not used?)
 psEllipsePol pmPSF_ModelToFit (psF32 *modelPar)
 {
@@ -329 +331 @@
 // convert the parameters used in the fitted source model to the psEllipseAxes representation
 // (major,minor,theta)
-psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR)
+psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR, pmModelType type)
 {
     psEllipseShape shape;
@@ -336 +338 @@
     axes.minor = NAN;
     axes.theta = NAN;
-//   XXX: must assert non-NULL input parameter
+    //   XXX: must assert non-NULL input parameter
     PS_ASSERT_PTR_NON_NULL(modelPar, axes);
 
-    shape.sx  = modelPar[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = modelPar[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = modelPar[PM_PAR_SXY];
+    bool useReff = true;
+    useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC"));
+    useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV"));
+    useReff |= (type == pmModelClassGetType ("PS_MODEL_EXP"));
+
+    if (useReff) {
+        shape.sx  = modelPar[PM_PAR_SXX];
+        shape.sy  = modelPar[PM_PAR_SYY];
+        shape.sxy = modelPar[PM_PAR_SXY];
+    } else {
+        shape.sx  = modelPar[PM_PAR_SXX] / M_SQRT2;
+        shape.sy  = modelPar[PM_PAR_SYY] / M_SQRT2;
+        shape.sxy = modelPar[PM_PAR_SXY];
+    }
 
     if ((shape.sx == 0) || (shape.sy == 0)) {
@@ -357 +370 @@
 // convert the psEllipseAxes representation (major,minor,theta) to the parameters used in the
 // fitted source model
-bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes)
+bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, pmModelType type)
 {
     PS_ASSERT_PTR_NON_NULL(modelPar, false);
@@ -370 +383 @@
     psEllipseShape shape = psEllipseAxesToShape (axes);
 
-    modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2;
-    modelPar[PM_PAR_SYY] = shape.sy * M_SQRT2;
-    modelPar[PM_PAR_SXY] = shape.sxy;
-
+    bool useReff = true;
+    useReff |= pmModelClassGetType ("PS_MODEL_SERSIC");
+    useReff |= pmModelClassGetType ("PS_MODEL_DEV");
+    useReff |= pmModelClassGetType ("PS_MODEL_EXP");
+
+    if (useReff) {
+        modelPar[PM_PAR_SXX] = shape.sx;
+        modelPar[PM_PAR_SYY] = shape.sy;
+        modelPar[PM_PAR_SXY] = shape.sxy;
+    } else {
+        modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2;
+        modelPar[PM_PAR_SYY] = shape.sy * M_SQRT2;
+        modelPar[PM_PAR_SXY] = shape.sxy;
+    }
     return true;
 }
@@ -438 +461 @@
 # if (0)
     psF32 *params = model->params->data.F32; // Model parameters
-    psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO); // Ellipse axes
+    psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO, model->type); // Ellipse axes
 
     // Curiously, the minor axis can be larger than the major axis, so need to check.

trunk/psModules/src/objects/pmPSF.h

@@ -106 +106 @@
 pmPSF *pmPSFBuildSimple (char *typeName, float sxx, float syy, float sxy, ...);
 
-bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes);
+bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, pmModelType type);
 bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis);
 
 psEllipsePol pmPSF_ModelToFit (psF32 *modelPar);
-psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR);
+psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR, pmModelType type);
 
 /// Calculate FWHM value from a PSF

trunk/psModules/src/objects/pmPeaks.c

@@ -137 +137 @@
 *****************************************************************************/
 static void peakFree(pmPeak *tmp)
-{} //
+{
+    if (!tmp) return;
+    psFree (tmp->saddlePoints);
+    return;
+}
 
 pmPeak *pmPeakAlloc(psS32 x,
@@ -162 +166 @@
     tmp->type = type;
     tmp->footprint = NULL;
+    tmp->saddlePoints = NULL;
 
     psMemSetDeallocator(tmp, (psFreeFunc) peakFree);

trunk/psModules/src/objects/pmPeaks.h

@@ -72 +72 @@
     pmPeakType type;                    ///< Description of peak.
     pmFootprint *footprint;             ///< reference to containing footprint
+    psArray *saddlePoints;              ///< set of saddle points between this peak and near neighbors
 }
 pmPeak;

trunk/psModules/src/objects/pmSource.c

@@ -58 +58 @@
     psFree(tmp->modelEXT);
     psFree(tmp->modelFits);
+    psFree(tmp->extFitPars);
     psFree(tmp->extpars);
     psFree(tmp->moments);
@@ -119 +120 @@
     source->modelEXT = NULL;
     source->modelFits = NULL;
+    source->extFitPars = NULL;
     source->type = PM_SOURCE_TYPE_UNKNOWN;
     source->mode = PM_SOURCE_MODE_DEFAULT;
@@ -196 +198 @@
     // NOTE : because of the const id element, we cannot just assign *source = *in
 
+    source->imageID          = in->imageID;
     source->type             = in->type;
     source->mode             = in->mode;

trunk/psModules/src/objects/pmSource.h

@@ -36 +36 @@
     PM_SOURCE_TMPF_MOMENTS_MEASURED  = 0x0010,
     PM_SOURCE_TMPF_CANDIDATE_PSFSTAR = 0x0020,
+    PM_SOURCE_TMPF_STACK_KEEP        = 0x0040,
+    PM_SOURCE_TMPF_STACK_SKIP        = 0x0080,
 } pmSourceTmpF;
 
@@ -76 +78 @@
     pmModel *modelEXT;                  ///< EXT Model fit used for subtraction (parameters and type)
     psArray *modelFits;                 ///< collection of extended source models (best == modelEXT)
+    psArray *extFitPars;                ///< extra extended fit parameters
     pmSourceType type;                  ///< Best identification of object.
     pmSourceMode mode;                  ///< analysis flags set for object.

trunk/psModules/src/objects/pmSourceExtendedPars.c

@@ -258 +258 @@
     return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceExtendedFluxFree);
 }
+
+// *** pmSourceExtFitPars describes extra metadata related to an extended fit
+static void pmSourceExtFitParsFree (pmSourceExtFitPars *pars) {
+    return;
+}
+
+pmSourceExtFitPars *pmSourceExtFitParsAlloc (void) {
+
+    pmSourceExtFitPars *pars = (pmSourceExtFitPars *) psAlloc(sizeof(pmSourceExtFitPars));
+    psMemSetDeallocator(pars, (psFreeFunc) pmSourceExtFitParsFree);
+
+    pars->Mxx = NAN;
+    pars->Mxy = NAN;
+    pars->Myy = NAN;
+
+    pars->Mrf    = NAN;
+    pars->Mrh    = NAN;
+
+    pars->apMag  = NAN;
+    pars->krMag  = NAN;
+    pars->psfMag = NAN;
+    pars->peakMag = NAN;
+
+    return pars;
+}

trunk/psModules/src/objects/pmSourceExtendedPars.h

@@ -67 +67 @@
 } pmSourceExtendedPars;
 
+// additional measurements related to the model fits
+typedef struct {
+    float Mxx;
+    float Mxy;
+    float Myy;
+    
+    float Mrf;
+    float Mrh;
+
+    float apMag;
+    float krMag;
+    float psfMag;
+    float peakMag;
+} pmSourceExtFitPars;
+
 pmSourceRadialFlux *pmSourceRadialFluxAlloc();
 bool psMemCheckSourceRadialFlux(psPtr ptr);
@@ -92 +107 @@
 bool pmSourceRadialProfileSortPair(psVector *index, psVector *extra);
 
-
+pmSourceExtFitPars *pmSourceExtFitParsAlloc (void);
 
 /// @}

trunk/psModules/src/objects/pmSourceFitModel.c

@@ -176 +176 @@
         break;
       case PM_SOURCE_FIT_EXT:
-        // EXT model fits all params (except sky)
-        nParams = params->n - 1;
+        // EXT model fits all shape params and Io (not Xo, Yo, sky)
+        nParams = params->n - 3;
         psVectorInit (constraint->paramMask, 0);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 1;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 1;
         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
         break;
@@ -193 +195 @@
         break;
       case PM_SOURCE_FIT_NO_INDEX:
-        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
+        // PSF model only fits Io, Sxx, Sxy, Syy
         psVectorInit (constraint->paramMask, 0);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 1;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 1;
         constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
         if (params->n == 7) {
-            nParams = params->n - 1;
+            nParams = params->n - 3;
         } else {
-            nParams = params->n - 2;
+            nParams = params->n - 4;
             constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
         }

trunk/psModules/src/objects/pmSourceFitPCM.c

@@ -48 +48 @@
 // convolved model image.
 
+# define TIMING 0
+
 bool pmSourceFitPCM (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
     
+    if (TIMING) { psTimerStart ("pmSourceFitPCM"); }
+
     psVector *params  = pcm->modelConv->params;
     psVector *dparams  = pcm->modelConv->dparams;
@@ -64 +68 @@
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
 
+    float t1, t2, t3, t4, t5;
+    if (TIMING) { t1 = psTimerMark ("pmSourceFitPCM"); }
+
     bool fitStatus = pmPCM_MinimizeChisq (myMin, covar, params, source, pcm);
+    if (TIMING) { t2 = psTimerMark ("pmSourceFitPCM"); }
+
     for (int i = 0; i < dparams->n; i++) {
         if ((pcm->constraint->paramMask != NULL) && pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
@@ -76 +85 @@
     }
     psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
+    if (TIMING) { t3 = psTimerMark ("pmSourceFitPCM"); }
 
     // renormalize output model image (generated by fitting process)
@@ -97 +107 @@
         pmSourceChisqUnsubtracted (source, pcm->modelConv, maskVal);
     }
+    if (TIMING) { t4 = psTimerMark ("pmSourceFitPCM"); }
 
     // set the model success or failure status
@@ -114 +125 @@
 
     source->mode |= PM_SOURCE_MODE_FITTED; // XXX is this needed?
+    if (TIMING) { t5 = psTimerMark ("pmSourceFitPCM"); }
+
+     if (TIMING) {
+        fprintf (stderr, "nIter: %2d, npix: %5d, t1: %6.4f, t2: %6.4f, t3: %6.4f, t4: %6.4f, t5: %6.4f\n", myMin->iter, pcm->nPix, t1, t2, t3, t4, t5);
+     }
 
     psFree(myMin);
@@ -121 +137 @@
 }
 
+// XXX deprecate this function or merge with the empirical model
 bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) {
 

trunk/psModules/src/objects/pmSourceIO_CMF.c.in

@@ -55 +55 @@
 // followed by a zero-size matrix, followed by the table data
 
-// # define MODE @CMFMODE@
 bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
 {
@@ -171 +170 @@
         @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
 
+        // XXX do not keep this long term, just a TEST:
+        // @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_PSF",    PS_DATA_F32, "Kron Flux",                                  moments.KronPSF);
+        // @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_PSF_SIG",PS_DATA_F32, "Kron Flux",                                  moments.KronPSFErr);
         // Do NOT write these : not consistent with the definition of PS1_V3 in Ohana/src/libautocode/dev/cmf-ps1-v3.d
         // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX",   PS_DATA_F32, "Kron Flux (in 1.0 R1)",                      moments.KronCore);
         // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR",  PS_DATA_F32, "Kron Error (in 1.0 R1)",                     moments.KronCoreErr);
+
         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
         @=PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
@@ -222 +225 @@
 
 // read in a readout from the fits file
-psArray *pmSourcesRead_CMF_PS1_V3 (psFits *fits, psMetadata *header)
+psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header)
 {
     PS_ASSERT_PTR_NON_NULL(fits, false);
@@ -281 +284 @@
         // XXX use these to determine PAR[PM_PAR_I0]?
         @ALL@     source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        @ALL@     source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        @ALL@     source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         @ALL@     source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
+        @=PS1_V3@ source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
+
+        // XXX use these to determine PAR[PM_PAR_I0] if they exist?
+        @=PS1_V3@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
+        @=PS1_V3@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
@@ -288 +296 @@
         @ALL@     dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
-        pmPSF_AxesToModel (PAR, axes);
+        pmPSF_AxesToModel (PAR, axes, modelType);
 
         @ALL@     float peakMag     = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
@@ -353 +361 @@
 }
 
-bool pmSourcesWrite_CMF_PS1_V3_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+bool pmSourcesWrite_CMF_@CMFMODE@_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
 {
     bool status;
@@ -541 +549 @@
 
 // XXX this layout is still the same as PS1_DEV_1
-bool pmSourcesWrite_CMF_PS1_V3_XFIT (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
+bool pmSourcesWrite_CMF_@CMFMODE@_XFIT (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
 {
 
@@ -590 +598 @@
             assert (model);
 
+            // pmSourceExtFitPars *extPars = source->extFitPars->data[j];
+            // assert (extPars);
+
             // skip models which were not actually fitted
             if (model->flags & PM_MODEL_STATUS_BADARGS) continue;
@@ -600 +611 @@
             yErr = dPAR[PM_PAR_YPOS];
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
+            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+
+            float kronFlux = source->moments ? source->moments->KronFlux : NAN;
+            float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN;
 
             row = psMetadataAlloc ();
@@ -612 +626 @@
             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
+
+            // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XX",       0, "second moment in x",                      extPars->Mxx);
+            // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XY",       0, "second moment in x,y",                    extPars->Mxy);
+            // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_YY",       0, "second moment in y",                      extPars->Myy);
+            // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_R1",       0, "first radial moment",                     extPars->Mrf);
+            // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_RH",       0, "half radial moment",                      extPars->Mrh);
+
+            psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_INST_MAG",     0, "PSF fit instrumental magnitude",             source->psfMag);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "AP_MAG",           0, "PSF-sized aperture magnitude",               source->apMag);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "KRON_MAG",         0, "Kron Mag",                                   kronMag);
 
             psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
@@ -673 +697 @@
 }
 
-bool pmSourcesWrite_CMF_PS1_V3_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
+bool pmSourcesWrite_CMF_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
 {
     return true;

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c

@@ -263 +263 @@
         dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
-        pmPSF_AxesToModel (PAR, axes);
+        pmPSF_AxesToModel (PAR, axes, modelType);
 
         float peakMag     = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
@@ -547 +547 @@
             yErr = dPAR[PM_PAR_YPOS];
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
+            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
 
             row = psMetadataAlloc ();

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c

@@ -281 +281 @@
         dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
-        pmPSF_AxesToModel (PAR, axes);
+        pmPSF_AxesToModel (PAR, axes, modelType);
 
         float peakMag     = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
@@ -572 +572 @@
             yErr = dPAR[PM_PAR_YPOS];
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
+            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
 
             row = psMetadataAlloc ();

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

@@ -280 +280 @@
         dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
-        pmPSF_AxesToModel (PAR, axes);
+        pmPSF_AxesToModel (PAR, axes, modelType);
 
         float peakMag     = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
@@ -607 +607 @@
             yErr = dPAR[PM_PAR_YPOS];
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
+            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
 
             row = psMetadataAlloc ();

trunk/psModules/src/objects/pmSourceIO_CMP.c

@@ -135 +135 @@
         lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky);
 
-        axes = pmPSF_ModelToAxes (PAR, 20.0);
+        axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
 
         float psfMagErr = isfinite(source->psfMagErr) ? source->psfMagErr : 999;
@@ -293 +293 @@
                 goto skip_source;
 
-            pmPSF_AxesToModel (PAR, axes);
+            pmPSF_AxesToModel (PAR, axes, modelType);
 
             psArrayAdd (sources, 100, source);

trunk/psModules/src/objects/pmSourceIO_OBJ.c

@@ -91 +91 @@
         }
 
-        axes = pmPSF_ModelToAxes (PAR, 20.0);
+        axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
 
         psLineInit (line);

trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c

@@ -113 +113 @@
             yErr = dPAR[PM_PAR_YPOS];
             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0);
+                axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
             } else {
                 axes.major = NAN;
@@ -288 +288 @@
         dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
-        pmPSF_AxesToModel (PAR, axes);
+        pmPSF_AxesToModel (PAR, axes, modelType);
 
         float peakMag     = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
@@ -618 +618 @@
             yErr = dPAR[PM_PAR_YPOS];
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
+            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
 
             // generate RA,DEC

trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c

@@ -89 +89 @@
             yErr = dPAR[PM_PAR_YPOS];
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
+            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
         } else {
             // XXX: This code seg faults if source->peak is NULL.
@@ -214 +214 @@
         source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
 
-        pmPSF_AxesToModel (PAR, axes);
+        pmPSF_AxesToModel (PAR, axes, modelType);
 
         float peakMag = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");

trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

@@ -95 +95 @@
             yErr = dPAR[PM_PAR_YPOS];
             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0);
+                axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
             } else {
                 axes.major = NAN;
@@ -257 +257 @@
         dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
-        pmPSF_AxesToModel (PAR, axes);
+        pmPSF_AxesToModel (PAR, axes, modelType);
 
         float peakMag     = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
@@ -522 +522 @@
             yErr = dPAR[PM_PAR_YPOS];
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
+            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
 
             row = psMetadataAlloc ();

trunk/psModules/src/objects/pmSourceIO_SMPDATA.c

@@ -92 +92 @@
             lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky);
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
+            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
 
         } else {
@@ -190 +190 @@
         axes.theta       = psMetadataLookupF32 (&status, row, "THETA");
 
-        pmPSF_AxesToModel (PAR, axes);
+        pmPSF_AxesToModel (PAR, axes, modelType);
 
         source->psfMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT;

trunk/psModules/src/objects/pmSourceIO_SX.c

@@ -81 +81 @@
         // pmSourceSextractType (source, &type, &flags);
 
-        axes = pmPSF_ModelToAxes (PAR, 20.0);
+        axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
 
         psLineInit (line);

trunk/psModules/src/objects/pmSourceMoments.c

@@ -96 +96 @@
     // (int) so they can be used in the image index below.
 
-    // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to 
-    // get an unbiased (but probably noisy) centroid
-    if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
-        return false;
-    }
-    // second pass applies the Gaussian window and uses the centroid from the first pass
-    if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) {
+    // XXX // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to 
+    // XXX // get an unbiased (but probably noisy) centroid
+    // XXX if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
+    // XXX      return false;
+    // XXX }
+    // XXX // second pass applies the Gaussian window and uses the centroid from the first pass
+    // XXX if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) {
+    // XXX      return false;
+    // XXX }
+
+    // If we use a large radius for the centroid, it will be biased by any neighbors.  The flux
+    // 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)
+    if (!pmSourceMomentsGetCentroid (source, 1.5*sigma, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
         return false;
     }
@@ -108 +115 @@
     // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
     // Xn  = SUM (x - xc)^n * (z - sky)
-
-    float RFW = 0.0;
-    float RHW = 0.0;
-
-    float RF = 0.0;
-    float RH = 0.0;
-    float RS = 0.0;
     float XX = 0.0;
     float XY = 0.0;
@@ -143 +143 @@
     float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
 
+    // calculate the higher-order moments using Xo,Yo
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
@@ -194 +195 @@
             Sum += pDiff;
 
-            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
             float r = sqrt(r2);
-            float rf = r * fDiff;
-            float rh = sqrt(r) * fDiff;
-            float rs = fDiff;
-
-            float rfw = r * pDiff;
-            float rhw = sqrt(r) * pDiff;
 
             float x = xDiff * pDiff;
@@ -221 +215 @@
             float yyyy = yDiff * yyy / r2;
 
-            RF  += rf;
-            RH  += rh;
-            RS  += rs;
-
-            RFW  += rfw;
-            RHW  += rhw;
-
             XX  += xx;
             XY  += xy;
@@ -242 +229 @@
             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;
         }
     }
-
-    source->moments->Mrf = RF/RS;
-    source->moments->Mrh = RH/RS;
-
     source->moments->Mxx = XX/Sum;
     source->moments->Mxy = XY/Sum;
@@ -262 +261 @@
     source->moments->Mxyyy = XYYY/Sum;
     source->moments->Myyyy = YYYY/Sum;
+
+    // *** now calculate the 1st radial moment (for kron flux) -- symmetrical averaging
+
+    float **vPix = source->pixels->data.F32;
+    float **vWgt = source->variance->data.F32;
+    psImageMaskType  **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+
+    float RF = 0.0;
+    float RH = 0.0;
+    float RS = 0.0;
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+
+        float yDiff = row - yCM;
+        if (fabs(yDiff) > radius) continue;
+
+        // coordinate of mirror pixel
+        int yFlip = yCM - yDiff;
+        if (yFlip < 0) continue;
+        if (yFlip >= source->pixels->numRows) continue;
+
+        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+            // check mask and value for this pixel
+            if (vMsk && (vMsk[row][col] & maskVal)) continue;
+            if (isnan(vPix[row][col])) continue;
+
+            float xDiff = col - xCM;
+            if (fabs(xDiff) > radius) continue;
+
+            // coordinate of mirror pixel
+            int xFlip = xCM - xDiff;
+            if (xFlip < 0) continue;
+            if (xFlip >= source->pixels->numCols) continue;
+
+            // check mask and value for mirror pixel
+            if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
+            if (isnan(vPix[yFlip][xFlip])) continue;
+
+            // radius is just a function of (xDiff, yDiff)
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            if (r2 > R2) continue;
+
+            float fDiff1 = vPix[row][col] - sky;
+            float fDiff2 = vPix[yFlip][xFlip] - sky;
+            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
+
+            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
+            float r = sqrt(r2);
+            float rf = r * pDiff;
+            float rh = sqrt(r) * pDiff;
+            float rs = 0.5 * (fDiff1 + fDiff2);
+
+            RF  += rf;
+            RH  += rh;
+            RS  += rs;
+        }
+    }
+
+    source->moments->Mrf = RF/RS;
+    source->moments->Mrh = RH/RS;
 
     // if Mrf (first radial moment) is very small, we are getting into low-significance
@@ -270 +329 @@
         kronRefRadius = MIN(radius, kronRefRadius);
     }
+
+    // *** now calculate the kron flux values using the 1st radial moment
 
     float radKinner = 1.0*kronRefRadius;
@@ -283 +344 @@
     float SumOuter = 0.0;
 
+    // calculate the Kron flux, and related fluxes (NO symmetrical averaging)
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
-
+        
         float yDiff = row - yCM;
         if (fabs(yDiff) > radKouter) continue;
+        
+        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+            // check mask and value for this pixel
+            if (vMsk && (vMsk[row][col] & maskVal)) continue;
+            if (isnan(vPix[row][col])) continue;
+            
+            float xDiff = col - xCM;
+            if (fabs(xDiff) > radKouter) continue;
+            
+            // radKron is just a function of (xDiff, yDiff)
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+
+            float fDiff1 = vPix[row][col] - sky;
+            float pDiff = fDiff1;
+            float wDiff = vWgt[row][col];
+                                    
+            // skip pixels below specified significance level.  this is allowed, but should be
+            // avoided -- the over-weights the wings of bright stars compared to those of faint
+            // stars.
+            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
+            
+            float r  = sqrt(r2);
+            if (r < radKron) {
+                Sum += pDiff;
+                Var += wDiff;
+                nKronPix ++;
+                // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
+            }
+
+            // use sigma (fixed by psf) not a radKron based value
+            if (r < sigma) {
+                SumCore += pDiff;
+                VarCore += wDiff;
+                nCorePix ++;
+            }
+
+            if ((r > radKinner) && (r < radKron)) {
+                SumInner += pDiff;
+                nInner ++;
+            }
+            if ((r > radKron)  && (r < radKouter)) {
+                SumOuter += pDiff;
+                nOuter ++;
+            }
+        }
+    }
+    // *** 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;
+
+    source->moments->KronCore    = SumCore;
+    source->moments->KronCoreErr = sqrt(VarCore);
+    source->moments->KronFlux    = Sum;
+    source->moments->KronFluxErr = sqrt(Var);
+    source->moments->KronFinner = SumInner;
+    source->moments->KronFouter = SumOuter;
+
+    // XXX not sure I should save this here...
+    source->moments->KronFluxPSF    = source->moments->KronFlux;
+    source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
+    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->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];
+        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++) {
@@ -304 +481 @@
             if (isnan(*vPix)) continue;
 
-            float xDiff = col - xCM;
-            if (fabs(xDiff) > radKouter) continue;
-
-            // radKron is just a function of (xDiff, yDiff)
-            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-
-            float pDiff = *vPix - sky;
-            float wDiff = *vWgt;
-
-            // skip pixels below specified significance level.  this is allowed, but should be
-            // avoided -- the over-weights the wings of bright stars compared to those of faint
-            // stars.
-            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
-
-            float r  = sqrt(r2);
-            if (r < radKron) {
-                Sum += pDiff;
-                Var += wDiff;
-                nKronPix ++;
-                // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
-            }
-
-            // use sigma (fixed by psf) not a radKron based value
-            if (r < sigma) {
-                SumCore += pDiff;
-                VarCore += wDiff;
-                nCorePix ++;
-            }
-
-            if ((r > radKinner) && (r < radKron)) {
-                SumInner += pDiff;
-                nInner ++;
-            }
-            if ((r > radKron)  && (r < radKouter)) {
-                SumOuter += pDiff;
-                nOuter ++;
-            }
-        }
-    }
-    // *** should I rescale these fluxes by pi R^2 / nNpix?
-    source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
-    source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
-    source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
-    source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
-    source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
-    source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
-
-    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->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;

trunk/psModules/src/objects/pmSourceOutputs.c

@@ -107 +107 @@
         }
         if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXY]) && isfinite(PAR[PM_PAR_SYY])) {
-            axes = pmPSF_ModelToAxes (PAR, 20.0);
+            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
             outputs->psfMajor = axes.major;
             outputs->psfMinor = axes.minor;
@@ -178 +178 @@
     moments->KronCore    = source->moments ? source->moments->KronCore : NAN;
     moments->KronCoreErr = source->moments ? source->moments->KronCoreErr : NAN;
-
-    return true;
-}
+    moments->KronPSF    = source->moments ? source->moments->KronFluxPSF : NAN;
+    moments->KronPSFErr = source->moments ? source->moments->KronFluxPSFErr : NAN;
+
+    return true;
+}
+
+bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos) {
+
+    pmFPA *fpa = chip->parent;
+
+    if (!chip->toFPA) goto escape;
+    if (!fpa->toTPA) goto escape;
+    if (!fpa->toSky) goto escape;
+
+    // generate RA,DEC
+    psPlane ptCH, ptFP, ptTP_o, ptTP_x, ptTP_y;
+
+    // calculate the astrometry for the coordinate of interest
+    ptCH.x = xPos;
+    ptCH.y = yPos;
+    psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
+    psPlaneTransformApply (&ptTP_o, fpa->toTPA, &ptFP);
+    psDeproject (ptSky, &ptTP_o, fpa->toSky);
+
+    // calculate the astrometry for the coordinate + 1pix in X
+    ptCH.x = xPos + 1.0;
+    ptCH.y = yPos;
+    psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
+    psPlaneTransformApply (&ptTP_x, fpa->toTPA, &ptFP);
+
+    // calculate the astrometry for the coordinate + 1pix in Y
+    ptCH.x = xPos;
+    ptCH.y = yPos + 1.0;
+    psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
+    psPlaneTransformApply (&ptTP_y, fpa->toTPA, &ptFP);
+
+    // the resulting Tangent Plane coordinates are in TP pixels; convert to local Tangent Plane
+    // degrees
+
+    float dTPx_dCHx = fpa->toSky->Xs * (ptTP_x.x - ptTP_o.x);
+    float dTPy_dCHx = fpa->toSky->Ys * (ptTP_x.y - ptTP_o.y);
+
+    float dTPx_dCHy = fpa->toSky->Xs * (ptTP_y.x - ptTP_o.x);
+    float dTPy_dCHy = fpa->toSky->Ys * (ptTP_y.y - ptTP_o.y);
+
+    float pltScale_x = hypot(dTPx_dCHx, dTPy_dCHx);
+    float pltScale_y = hypot(dTPx_dCHy, dTPy_dCHy);
+    *pltScale = 0.5*(pltScale_x + pltScale_y);
+
+    float posAngle_x = atan2 (+dTPy_dCHx, +dTPx_dCHx);
+    float posAngle_y = atan2 (-dTPy_dCHy, +dTPx_dCHy);
+    *posAngle = 0.5*(posAngle_x + posAngle_y);
+
+    return true;
+
+escape:
+    // no astrometry calibration, give up
+    ptSky->r = NAN;
+    ptSky->d = NAN;
+    *posAngle = NAN;
+    *pltScale = NAN;
+
+    return false;
+}
+

trunk/psModules/src/objects/pmSourceOutputs.h

@@ -52 +52 @@
     float KronCore;
     float KronCoreErr;
+    float KronPSF;
+    float KronPSFErr;
 } pmSourceOutputsMoments;
 

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -94 +94 @@
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(psf, false);
+    // PS_ASSERT_PTR_NON_NULL(psf, false);
 
     int status = false;

trunk/psModules/src/objects/pmSourcePlotPSFModel.c

@@ -145 +145 @@
         // force the axis ratio to be < 20.0
         psEllipseAxes axes_mnt = psEllipseMomentsToAxes (moments, 20.0);
-        psEllipseAxes axes_psf = pmPSF_ModelToAxes (PAR, 20.0);
+        psEllipseAxes axes_psf = pmPSF_ModelToAxes (PAR, 20.0, model->type);
 
         // moments major axis

@@ -5 +5 @@
 *.la
 *.lo
+pmSourceIO_CMF_PS1_V1.c
+pmSourceIO_CMF_PS1_V2.c
+pmSourceIO_CMF_PS1_V3.c