Changeset 31153 for trunk/psModules/src/objects

Timestamp:

Apr 4, 2011, 1:04:41 PM (15 years ago)

Author:

eugene

Message:

updates to pmPeak to better distinguish peak flux versions; updates to visualization; add bits for substantial suspect masking; consolidate assignment of source position and flux based on peak, moments, etc; improve footprint culling process; fix PSF_QF and PSF_QF_PERFECT calculations; fix source model chisq values

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/objects/Makefile.am (modified) (2 diffs)
• src/objects/models/pmModel_DEV.c (modified) (2 diffs)
• src/objects/models/pmModel_EXP.c (modified) (2 diffs)
• src/objects/models/pmModel_GAUSS.c (modified) (2 diffs)
• src/objects/models/pmModel_PGAUSS.c (modified) (2 diffs)
• src/objects/models/pmModel_PS1_V1.c (modified) (2 diffs)
• src/objects/models/pmModel_QGAUSS.c (modified) (2 diffs)
• src/objects/models/pmModel_RGAUSS.c (modified) (2 diffs)
• src/objects/models/pmModel_SERSIC.c (modified) (4 diffs)
• src/objects/pmFootprint.c (modified) (1 diff)
• src/objects/pmFootprint.h (modified) (2 diffs)
• src/objects/pmFootprintAssignPeaks.c (modified) (3 diffs)
• src/objects/pmFootprintCullPeaks.c (modified) (6 diffs)
• src/objects/pmFootprintFindAtPoint.c (modified) (6 diffs)
• src/objects/pmFootprintIDs.c (modified) (1 diff)
• src/objects/pmModel.c (modified) (3 diffs)
• src/objects/pmModel.h (modified) (2 diffs)
• src/objects/pmModelUtils.c (modified) (1 diff)
• src/objects/pmModelUtils.h (modified) (1 diff)
• src/objects/pmMoments.h (modified) (1 diff)
• src/objects/pmPCMdata.c (modified) (2 diffs)
• src/objects/pmPCMdata.h (modified) (1 diff)
• src/objects/pmPSF.c (modified) (2 diffs)
• src/objects/pmPSFtryFitPSF.c (modified) (1 diff)
• src/objects/pmPeaks.c (modified) (7 diffs)
• src/objects/pmPeaks.h (modified) (4 diffs)
• src/objects/pmPhotObj.c (modified) (3 diffs)
• src/objects/pmPhotObj.h (modified) (2 diffs)
• src/objects/pmSource.c (modified) (9 diffs)
• src/objects/pmSource.h (modified) (2 diffs)
• src/objects/pmSourceExtendedPars.c (modified) (2 diffs)
• src/objects/pmSourceExtendedPars.h (modified) (2 diffs)
• src/objects/pmSourceFitModel.c (modified) (3 diffs)
• src/objects/pmSourceFitModel.h (modified) (1 diff)
• src/objects/pmSourceFitPCM.c (modified) (3 diffs)
• src/objects/pmSourceFitSet.c (modified) (11 diffs)
• src/objects/pmSourceFitSet.h (modified) (2 diffs)
• src/objects/pmSourceIO_CMF_PS1_DV1.c (modified) (10 diffs)
• src/objects/pmSourceIO_CMF_PS1_DV2.c (modified) (11 diffs)
• src/objects/pmSourceIO_CMF_PS1_SV1.c (modified) (20 diffs)
• src/objects/pmSourceIO_CMF_PS1_V1.c (modified) (9 diffs)
• src/objects/pmSourceIO_CMF_PS1_V2.c (modified) (9 diffs)
• src/objects/pmSourceIO_CMF_PS1_V3.c (modified) (10 diffs)
• src/objects/pmSourceIO_PS1_CAL_0.c (modified) (5 diffs)
• src/objects/pmSourceIO_PS1_DEV_0.c (modified) (2 diffs)
• src/objects/pmSourceIO_PS1_DEV_1.c (modified) (5 diffs)
• src/objects/pmSourceIO_RAW.c (modified) (5 diffs)
• src/objects/pmSourceIO_SMPDATA.c (modified) (1 diff)
• src/objects/pmSourceMasks.h (modified) (1 diff)
• src/objects/pmSourceMoments.c (modified) (6 diffs)
• src/objects/pmSourceOutputs.c (copied) (copied from branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.c )
• src/objects/pmSourceOutputs.h (copied) (copied from branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.h )
• src/objects/pmSourcePhotometry.c (modified) (11 diffs)
• src/objects/pmSourcePhotometry.h (modified) (1 diff)
• src/objects/pmSourceVisual.c (modified) (5 diffs)

trunk/psModules

@@ -28 +28 @@
 ChangeLog
 psmodules-*.tar.*
+a.out.dSYM

@@ -28 +28 @@
 ChangeLog
 psmodules-*.tar.*
+a.out.dSYM

trunk/psModules/src/objects/Makefile.am

@@ -32 +32 @@
         pmSourceFitSet.c \
         pmSourcePhotometry.c \
+        pmSourceOutputs.c \
         pmSourceIO.c \
         pmSourceIO_RAW.c \
@@ -102 +103 @@
         pmSourceFitSet.h \
         pmSourcePhotometry.h \
+        pmSourceOutputs.h \
         pmSourceIO.h \
         pmSourcePlots.h \

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

@@ -217 +217 @@
 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
 {
-    pmMoments *moments = source->moments;
-    pmPeak    *peak    = source->peak;
-    psF32     *PAR  = model->params->data.F32;
-
-    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;
-
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
-    if (!isfinite(shape.sx))  return false;
-    if (!isfinite(shape.sy))  return false;
-    if (!isfinite(shape.sxy)) return false;
-
-    // the other parameters depend on the guess for PAR_7
+    psF32 *PAR  = model->params->data.F32;
+
+    // sky is set to 0.0
+    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)) {
+      return false;
+    }
+
+    // the normalization is modified by the slope
     float index = 0.5 / ALPHA;
     float bn = 1.9992*index - 0.3271;
-    // float fR = 1.0 / (sqrt(2.0) * pow (bn, index));
     float Io = exp(0.5*bn);
 
-    float Sxx = PS_MAX(0.5, shape.sx);
-    float Syy = PS_MAX(0.5, shape.sy);
-
-    PAR[PM_PAR_SKY]  = 0.0;
-    PAR[PM_PAR_I0]   = peak->flux / Io;
-    PAR[PM_PAR_XPOS] = peak->xf;
-    PAR[PM_PAR_YPOS] = peak->yf;
-    // PAR[PM_PAR_SXX]  = Sxx * fR;
-    // PAR[PM_PAR_SYY]  = Syy * fR;
-    PAR[PM_PAR_SXX]  = Sxx;
-    PAR[PM_PAR_SYY]  = Syy;
-    PAR[PM_PAR_SXY]  = shape.sxy;
+    // set the model normalization
+    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
+      return false;
+    }
+    PAR[PM_PAR_I0] /= Io;
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
 
     return(true);
@@ -322 +307 @@
     psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA);
 
-    psAssert (isfinite(radius), "fix this code: z should not be nan");
+    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
+              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
+
     return (radius);
 }

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

@@ -210 +210 @@
 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
 {
-    pmMoments *moments = source->moments;
-    pmPeak    *peak    = source->peak;
-    psF32     *PAR  = model->params->data.F32;
-
-    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;
-
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
-    if (!isfinite(shape.sx))  return false;
-    if (!isfinite(shape.sy))  return false;
-    if (!isfinite(shape.sxy)) return false;
-
+    psF32 *PAR  = model->params->data.F32;
+
+    // sky is set to 0.0
     PAR[PM_PAR_SKY]  = 0.0;
-    PAR[PM_PAR_I0]   = peak->flux;
-    PAR[PM_PAR_XPOS] = peak->xf;
-    PAR[PM_PAR_YPOS] = peak->yf;
-    PAR[PM_PAR_SXX]  = PS_MAX(0.5, M_SQRT2*shape.sx);
-    PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
-    PAR[PM_PAR_SXY]  = shape.sxy;
+
+    // set the shape parameters
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+      return false;
+    }
+
+    // set the model normalization
+    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
+      return false;
+    }
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
+
     return(true);
 }
@@ -303 +294 @@
     psF64 radius = axes.major * sqrt (2.0) * zn;
 
-    psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);
+    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
+              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
     return (radius);
 }

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

@@ -193 +193 @@
 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
 {
-    pmMoments *moments = source->moments;
-    pmPeak    *peak    = source->peak;
-    psF32     *PAR  = model->params->data.F32;
-
-    psEllipseMoments emoments;
-    emoments.x2 = moments->Mxx;
-    emoments.y2 = moments->Myy;
-    emoments.xy = moments->Mxy;
-
-    // force the axis ratio to be < 20.0
-    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
+    psF32 *PAR  = model->params->data.F32;
+
+    // sky is set to 0.0
     PAR[PM_PAR_SKY]  = 0.0;
-    PAR[PM_PAR_I0]   = peak->flux;
-    PAR[PM_PAR_XPOS] = peak->xf;
-    PAR[PM_PAR_YPOS] = peak->yf;
-    PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx);
-    PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy);
-    PAR[PM_PAR_SXY] = shape.sxy;
+
+    // set the shape parameters
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+      return false;
+    }
+
+    // set the model normalization
+    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
+      return false;
+    }
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
+
     return(true);
 }
@@ -258 +258 @@
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
-    psAssert (isfinite(radius), "fix this code: radius should not be nan for %f", PAR[PM_PAR_I0]);
+    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
+              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
 
     return (radius);

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

@@ -194 +194 @@
 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
 {
-    pmMoments *moments = source->moments;
-    pmPeak    *peak    = source->peak;
-    psF32     *PAR     = model->params->data.F32;
-
-    psEllipseMoments emoments;
-    emoments.x2 = moments->Mxx;
-    emoments.xy = moments->Mxy;
-    emoments.y2 = moments->Myy;
-
-    psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
+    psF32 *PAR  = model->params->data.F32;
+
+    // sky is set to 0.0
     PAR[PM_PAR_SKY]  = 0.0;
-    PAR[PM_PAR_I0]   = peak->flux;
-    PAR[PM_PAR_XPOS] = peak->xf;
-    PAR[PM_PAR_YPOS] = peak->yf;
-    PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx);
-    PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy);
-    PAR[PM_PAR_SXY] = shape.sxy;
+
+    // set the shape parameters
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+      return false;
+    }
+
+    // set the model normalization
+    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
+      return false;
+    }
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
+
     return(true);
 }
@@ -285 +286 @@
     // choose a z value guaranteed to be beyond our limit
     float z0 = pow((1.0 / limit), (1.0 / 3.0));
-    psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_I0]);
+    psAssert (isfinite(z0), "fix this code: z0 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
+              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
+
     float z1 = (1.0 / limit);
-    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_I0]);
+    psAssert (isfinite(z1), "fix this code: z1 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
+              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
+
     z1 = PS_MAX (z0, z1);
     z0 = 0.0;

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

@@ -213 +213 @@
 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
 {
-    pmMoments *moments = source->moments;
-    pmPeak    *peak    = source->peak;
-    psF32     *PAR  = model->params->data.F32;
-
-    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;
-
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
-    if (!isfinite(shape.sx))  return false;
-    if (!isfinite(shape.sy))  return false;
-    if (!isfinite(shape.sxy)) return false;
-
+    psF32 *PAR  = model->params->data.F32;
+
+    // sky is set to 0.0
     PAR[PM_PAR_SKY]  = 0.0;
-    PAR[PM_PAR_I0]   = peak->flux;
-    PAR[PM_PAR_XPOS] = peak->xf;
-    PAR[PM_PAR_YPOS] = peak->yf;
-    PAR[PM_PAR_SXX]  = PS_MAX(0.5, M_SQRT2*shape.sx);
-    PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
-    PAR[PM_PAR_SXY]  = shape.sxy;
+
+    // set the shape parameters
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+      return false;
+    }
+
+    // set the model normalization
+    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
+      return false;
+    }
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
+
+    // extra parameter
     PAR[PM_PAR_7]    = 0.5;
 
@@ -314 +306 @@
     float z0 = 0.0;
     float z1 = pow((1.0 / limit), (1.0 / ALPHA));
-    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
+    psAssert (isfinite(z1), "fix this code: z1 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
+              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
     if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
 

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

@@ -214 +214 @@
 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
 {
-    pmMoments *moments = source->moments;
-    pmPeak    *peak    = source->peak;
-    psF32     *PAR  = model->params->data.F32;
-
-    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;
-
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
-    if (!isfinite(shape.sx))  return false;
-    if (!isfinite(shape.sy))  return false;
-    if (!isfinite(shape.sxy)) return false;
-
+    psF32 *PAR  = model->params->data.F32;
+
+    // sky is set to 0.0
     PAR[PM_PAR_SKY]  = 0.0;
-    PAR[PM_PAR_I0]   = peak->flux;
-    PAR[PM_PAR_XPOS] = peak->xf;
-    PAR[PM_PAR_YPOS] = peak->yf;
-    PAR[PM_PAR_SXX]  = PS_MAX(0.5, M_SQRT2*shape.sx);
-    PAR[PM_PAR_SYY]  = PS_MAX(0.5, M_SQRT2*shape.sy);
-    PAR[PM_PAR_SXY]  = shape.sxy;
+
+    // set the shape parameters
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+      return false;
+    }
+
+    // set the model normalization
+    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
+      return false;
+    }
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
+
+    // extra parameter
     PAR[PM_PAR_7]    = 1.0;
 
@@ -315 +307 @@
     float z0 = 0.0;
     float z1 = pow((1.0 / limit), (1.0 / ALPHA));
-    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
+    psAssert (isfinite(z1), "fix this code: z1 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
+              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
     if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
 

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

@@ -203 +203 @@
 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
 {
-    pmMoments *moments = source->moments;
-    pmPeak    *peak    = source->peak;
-    psF32     *PAR  = model->params->data.F32;
-
-    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;
-
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
-    if (!isfinite(shape.sx))  return false;
-    if (!isfinite(shape.sy))  return false;
-    if (!isfinite(shape.sxy)) return false;
-
+    psF32 *PAR  = model->params->data.F32;
+
+    // sky is set to 0.0
     PAR[PM_PAR_SKY]  = 0.0;
-    PAR[PM_PAR_I0]   = peak->flux;
-    PAR[PM_PAR_XPOS] = peak->xf;
-    PAR[PM_PAR_YPOS] = peak->yf;
-    PAR[PM_PAR_SXX]  = PS_MAX(0.5, shape.sx);
-    PAR[PM_PAR_SYY]  = PS_MAX(0.5, shape.sy);
-    PAR[PM_PAR_SXY]  = shape.sxy;
+
+    // set the shape parameters
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+      return false;
+    }
+
+    // set the model normalization
+    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
+      return false;
+    }
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
+
+    // extra parameter
     PAR[PM_PAR_7]    = 1.5;
 
@@ -305 +297 @@
     // choose a z value guaranteed to be beyond our limit
     float z0 = pow((1.0 / limit), (1.0 / PAR[PM_PAR_7]));
-    psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_7]);
+    psAssert (isfinite(z0), "fix this code: z0 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f", 
+              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], PAR[PM_PAR_7]);
+
     float z1 = (1.0 / limit);
-    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
+    psAssert (isfinite(z1), "fix this code: z1 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f", 
+              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], PAR[PM_PAR_7]);
+
     z1 = PS_MAX (z0, z1);
     z0 = 0.0;

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

@@ -222 +222 @@
 {
     pmMoments *moments = source->moments;
-    pmPeak    *peak    = source->peak;
     psF32     *PAR  = model->params->data.F32;
+
+    // sky is set to 0.0
+    PAR[PM_PAR_SKY]  = 0.0;
 
     // the other parameters depend on the guess for PAR_7
@@ -236 +238 @@
     float Zero  = 1.16 - 0.615 * PAR[PM_PAR_7];
 
+    // Sersic shape is a bit special
     psEllipseMoments emoments;
     emoments.x2 = moments->Mxx;
@@ -273 +276 @@
     float Syy = PS_MAX(0.5, shape.sy);
 
-    PAR[PM_PAR_SKY]  = 0.0;
-    PAR[PM_PAR_I0]   = peak->flux / Io;
-    PAR[PM_PAR_XPOS] = peak->xf;
-    PAR[PM_PAR_YPOS] = peak->yf;
     PAR[PM_PAR_SXX]  = Sxx;
     PAR[PM_PAR_SYY]  = Syy;
     PAR[PM_PAR_SXY]  = shape.sxy;
+
+    // set the model normalization (adjust for Sersic best guess)
+    if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
+      return false;
+    }
+    PAR[PM_PAR_I0] /= Io;
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
 
     return(true);
@@ -346 +356 @@
     psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]);
 
-    // fprintf (stderr, "sersic model %f %f, n %f, radius: %f, zn: %f, f/Io: %f, major: %f\n", PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], PAR[PM_PAR_7], radius, zn, flux/PAR[PM_PAR_I0], axes.major);
-
-    psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);
+    psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f), par 7 = %f", 
+              PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], PAR[PM_PAR_7]);
     return (radius);
 }

trunk/psModules/src/objects/pmFootprint.c

@@ -98 +98 @@
 }
 
+// XXX not actually used anywhere
 pmFootprint *pmFootprintNormalize(pmFootprint *fp) {
     if (fp != NULL && !fp->normalized) {
-        fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN);
+        if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE) {
+            fp->peaks = psArraySort(fp->peaks, pmPeaksSortBySmoothFluxDescend);
+        } else {
+            fp->peaks = psArraySort(fp->peaks, pmPeaksSortByRawFluxDescend);
+        }
         fp->normalized = true;
     }

trunk/psModules/src/objects/pmFootprint.h

@@ -10 +10 @@
 #ifndef PM_FOOTPRINT_H
 #define PM_FOOTPRINT_H
+
+// We need to choose up front if the culling algorithm uses the raw or smoothed image.
+// depending on which we choose, we should produce sorted peaks based on peak->rawFlux or
+// peak->smoothFlux
+
+# define PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE 1
 
 typedef struct {
@@ -84 +90 @@
                                  const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
                                  const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
-                                 const float min_threshold // minimum permitted coll height
+                                 const float min_threshold, // minimum permitted coll height
+                                 const float max_threshold,// maximum permitted coll height
+                                 const bool isWeightVar // the input weight may be variance (sigma^2) or S/N (1/sigma)
     );
 

trunk/psModules/src/objects/pmFootprintAssignPeaks.c

@@ -46 +46 @@
      */
     psImage *ids = pmSetFootprintArrayIDs(footprints, true);
-    assert (ids != NULL);
-    assert (ids->type.type == PS_TYPE_S32);
-    const int row0 = ids->row0;
-    const int col0 = ids->col0;
-    const int numRows = ids->numRows;
-    const int numCols = ids->numCols;
+    if (ids) { assert (ids->type.type == PS_TYPE_S32); }
+    
+    const int row0 = ids ? ids->row0 : 0;
+    const int col0 = ids ? ids->col0 : 0;
+    const int numRows = ids ? ids->numRows : -1;
+    const int numCols = ids ? ids->numCols : -1;
 
     for (int i = 0; i < peaks->n; i++) {
@@ -58 +58 @@
         const int y = peak->y - row0;
         
-        assert (x >= 0 && x < numCols && y >= 0 && y < numRows);
-        int id = ids->data.S32[y][x - col0];
+        if (ids) { assert (x >= 0 && x < numCols && y >= 0 && y < numRows);}
+        int id = ids ? ids->data.S32[y][x - col0] : 0;
+        // XXX I think the '[x - col0]' above is just wrong (should be [x], but never gets triggerd.
 
         if (id == 0) {                  // peak isn't in a footprint, so make one for it
@@ -86 +87 @@
         if (fp->peaks->n == 1) continue;
 
-        fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN);
+        // make sure the peaks are sorted in a way consistent with our cull process
+        if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE) {
+            fp->peaks = psArraySort(fp->peaks, pmPeaksSortBySmoothFluxDescend);
+        } else {
+            fp->peaks = psArraySort(fp->peaks, pmPeaksSortByRawFluxDescend);
+        }
 
         // XXX check for an assert on duplicates (I don't think they can happen, but

trunk/psModules/src/objects/pmFootprintCullPeaks.c

@@ -38 +38 @@
                                  const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
                                  const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
-                                 const float min_threshold) { // minimum permitted coll height
+                                 const float min_threshold, // minimum permitted coll height
+                                 const float max_threshold,
+                                 const bool useSmoothedImage
+    ) { // maximum permitted coll height
     assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
     assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
@@ -44 +47 @@
     assert (fp != NULL);
 
-    if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
+    if (fp->peaks == NULL || fp->peaks->n < 2) { // nothing to do
         return PS_ERR_NONE;
     }
@@ -53 +56 @@
 
     psImage *subImg = psImageSubset((psImage *)img, subRegion);
-    psImage *subWt = psImageSubset((psImage *)weight, subRegion);
-    assert (subImg != NULL && subWt != NULL);
+    psAssert (subImg, "trouble making local subimage");
 
     psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32);
@@ -68 +70 @@
     psArrayAdd (brightPeaks, 128, fp->peaks->data[0]);
 
+    // XXX test point
+    // pmPeak *myPeak = fp->peaks->data[0];
+    // if ((fabs(myPeak->x - 205) < 100) && (fabs(myPeak->y - 806) < 100)) {
+    //  fprintf (stderr, "test peak\n");
+    // }
+
     // allocate the peakFootprint and peakFPSpans containers -- these are re-used by pmFootprintsFindAtPoint to minimize allocs in this function
     pmFootprint *peakFootprint = pmFootprintAlloc(fp->nspans, subImg);
@@ -77 +85 @@
     psImage *subMask = psImageCopy(NULL, subImg, PS_TYPE_IMAGE_MASK);
 
-    // The brightest peak is always safe; go through other peaks trying to cull them
-    for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
-        const pmPeak *peak = fp->peaks->data[i];
-        int x = peak->x - subImg->col0;
-        int y = peak->y - subImg->row0;
-        //
-        // Find the level nsigma below the peak that must separate the peak
-        // from any of its friends
-        //
-        assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
-
-        // const float stdev = sqrt(subWt->data.F32[y][x]);
-        // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
-
-        const float stdev = sqrt(subWt->data.F32[y][x]);
-        const float flux = subImg->data.F32[y][x];
-        const float fStdev = fabs(stdev/flux);
-        const float stdev_pad = fabs(flux * hypot(fStdev, fPadding));
-        // if flux is negative, careful with fStdev
-
-        float threshold = flux - nsigma_delta*stdev_pad;
-
-        if (isnan(threshold)) {
-            // min_threshold is assumed to be below the detection threshold,
-            // so all the peaks are pmFootprint, and this isn't the brightest
-            continue;
-        }
-
-        // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
-        if (threshold > subImg->data.F32[y][x]) {
-            threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
-        }
-
-        if (threshold < min_threshold) {
-            threshold = min_threshold;
-        }
-
-        // init peakFootprint here?
-        pmFootprintsFindAtPoint(peakFootprint, peakFPSpans, subImg, subMask, threshold, brightPeaks, peak->y, peak->x);
-        if (peakFPSpans->nStartSpans > 2000) {
-            // dumpfootprints(peakFootprint, peakFPSpans);
-            // fprintf (stderr, "big footprint %d : %d\n", peakFootprint->nspans, peakFPSpans->nStartSpans);
-            fprintf (stderr, "big test footprint: %f %f to %f %f (%d pix)\n", peakFootprint->bbox.x0, peakFootprint->bbox.y0, peakFootprint->bbox.x1, peakFootprint->bbox.y1, peakFootprint->npix);
-        }
-
-        // at this point brightPeaks only has the peaks brighter than the current
-
-        // we set the IDs to either 1 (in peak) or 0 (not in peak)
-        pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
-
-        // If this peak has not already been assigned to a source, then we can look for any
-        // brighter peaks within its footprint. Check if any of the previous (brighter) peaks
-        // are within the footprint of this peak If so, the current peak is bogus; drop it.
-        bool keep = true;
-        for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) {
-            const pmPeak *peak2 = fp->peaks->data[j];
-            int x2 = peak2->x - subImg->col0;
-            int y2 = peak2->y - subImg->row0;
-            if (idImg->data.S32[y2][x2] == IN_PEAK) {
-                // There's a brighter peak within the footprint above threshold; so cull our initial peak
-                keep = false;
-            }
-        }
-        if (!keep) {
-            psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak");  // we should not drop any already assigned peaks.
-            continue;
-        }
-
-        psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
+    // if we have many peaks in a large footprint, we waste a lot of time generating nearly identical footprints
+    // here we attempt to cull peaks drawing a single footprint for all peaks in some threshold range
+    // fprintf (stderr, "footprint: %d x %d : %d pix, %d peaks\n", subImg->numCols, subImg->numRows, fp->npix, (int) fp->peaks->n);
+    if ((fp->npix > 30000) && (fp->peaks->n > 10)) {
+
+        // max flux is above threshold for brightest peak
+        pmPeak *maxPeak = fp->peaks->data[0];
+        float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux;
+
+        // we have a relationship between the bin and the threshold of:
+        // threshold = 0.25 beta^2 bin^2 + minThreshold
+        // thus, the max bin is: sqrt(4.0*(maxThreshold - minThreshold)/ALPHA^2)
+
+# define ALPHA 0.1
+        
+        float beta = nsigma_delta * ALPHA;
+        float beta2 = PS_SQR(beta);
+        int nBins = sqrt(4.0*(maxFlux - min_threshold)/beta2) + 10; // let's be extra generous
+
+        // create a vector to store the threshold bins used for each peak
+        psVector *threshbins = psVectorAlloc(fp->peaks->n, PS_TYPE_S32);
+        threshbins->data.S32[0] = -1; // we skip the first peak
+
+        /// create a vector to track if a peak has been tried or not:
+        psVector *peaktried = psVectorAlloc(fp->peaks->n, PS_TYPE_U8);
+        psVectorInit(peaktried, 0);
+        peaktried->data.U8[0] = true; // we skip the first peak
+
+        psVector *threshbounds = psVectorAlloc(nBins, PS_TYPE_F32);
+        for (int i = 0; i < nBins; i++) {
+            threshbounds->data.F32[i] = 0.25*beta2*PS_SQR(i) + min_threshold;       
+        }
+        psAssert(threshbounds->data.F32[threshbounds->n-1] > maxFlux, "upper limit does not include max flux");
+
+        psHistogram *threshist = psHistogramAllocGeneric(threshbounds);
+
+        // assign the peaks to the histogram bins based on their nominal thresholsd
+        for (int i = 1; i < fp->peaks->n; i++) {
+            const pmPeak *peak = fp->peaks->data[i];
+            float flux = useSmoothedImage ? peak->smoothFlux : peak->rawFlux;
+            float stdev = useSmoothedImage ? peak->smoothFluxStdev : peak->rawFluxStdev;
+
+            // if flux is negative, careful with fStdev
+            const float fStdev = fabs(stdev/flux);
+            const float stdev_pad = fabs(flux * hypot(fStdev, fPadding));
+
+            float threshold = flux - nsigma_delta*stdev_pad;
+            psAssert(!isnan(threshold), "impossible");
+
+            if (threshold <= min_threshold) {
+                threshist->nums->data.F32[0] += 1.0;
+                threshbins->data.S32[i] = 0;
+                continue;
+            }
+            int bin = sqrt(4.0*(threshold - min_threshold)/beta2);
+            psAssert(bin >= 0, "impossible bin");
+
+            bin = PS_MIN(bin, threshist->nums->n - 1);
+            threshist->nums->data.F32[bin] += 1.0;
+            
+            // record the bin selected for this peak
+            threshbins->data.S32[i] = bin;
+        }
+
+        // XXX TEST: did we assign correctly
+# if (0)
+        int nPeaks = 1; // we don't put the brightest in the histogram
+        for (int i = 0; i < threshist->nums->n; i++) {
+            if (threshist->nums->data.F32[i] > 0) {
+                fprintf (stderr, "%f : %f : %d\n", threshist->bounds->data.F32[i], threshist->bounds->data.F32[i+1], (int)threshist->nums->data.F32[i]);
+                nPeaks += threshist->nums->data.F32[i];
+            }
+        }
+        fprintf (stderr, "%d peaks vs %d in histogram\n", (int) fp->peaks->n, nPeaks);
+# endif
+
+        // XXX for the moment, we will use the alternate cull for all peaks -- it might be
+        // faster to use the standard process for the peaks for which the threshold bin
+        // contains < N sources (N ~ 5-10?)
+
+        // loop over the threshold bins from brightest to faintest
+        for (int i = threshist->nums->n - 1; i >= 0; i--) {
+            if (threshist->nums->data.F32[i] == 0) continue;
+
+            // we are going to examine the footprints at this threshold
+            float threshold = 0.5*(threshist->bounds->data.F32[i] + threshist->bounds->data.F32[i+1]);
+                    
+            // generate all footprints corresponding to this threshold
+            psArray *myFP = pmFootprintsFind(subImg, threshold, 5);
+            if (!myFP) {
+                psWarning ("missing footprint?");
+                continue;
+            }
+            if (!myFP->n) {
+                psWarning ("empty footprint?");
+                continue;
+            }
+
+            // an array to track if the footprint has a brighter peak or not
+            psVector *found = psVectorAlloc(myFP->n + 1, PS_TYPE_U8);
+            psVectorInit(found, 0);
+            int nFound = 0;
+
+            // set IDs to distinguish the footprints
+            psImageInit(idImg, 0);
+            pmSetFootprintArrayIDsForImage(idImg, myFP, true);
+        
+            // check which footprints contain already-accepted (brighter) peaks 
+            // (we can give up if/when we found a peak for all footprints)
+            for (int j = 0; (j < brightPeaks->n) && (nFound < found->n); j++) {
+                const pmPeak *peak = brightPeaks->data[j];
+                int x = peak->x - subImg->col0;
+                int y = peak->y - subImg->row0;
+                int myID = idImg->data.S32[y][x];
+                psAssert(myID >= 0, "impossible");
+                psAssert(myID < found->n, "impossible");
+
+                if (myID == 0) continue; // bright peak is not in a footprint
+                if (found->data.U8[myID]) continue; // we already know this footprint contains a peak
+                found->data.U8[myID] = true;
+                nFound ++;
+            }
+        
+            // check the peaks from this threshold bin: if they land in a footprint which has
+            // been found, we should drop that peak.  otherwise, keep it
+            for (int j = 0; j < fp->peaks->n; j++) {
+                pmPeak *peak = fp->peaks->data[j];
+                if (peak->assigned) continue; // peak is already claimed by a source -- don't cull
+
+                // skip peaks if we've already tried them
+                if (peaktried->data.U8[j]) continue;
+
+                // is this peak in the threshold bin of interest?
+                if (threshbins->data.S32[j] != i) continue;
+                
+                // do not try this peak again
+                peaktried->data.U8[j] = true;
+
+                int x = peak->x - subImg->col0;
+                int y = peak->y - subImg->row0;
+                int myID = idImg->data.S32[y][x];
+                psAssert(myID < found->n, "impossible");
+
+                // a peak in this threshold bin without a valid footprint comes from a region
+                // with only a handful of pixels (1 or more from the peak itself).  It probably
+                // cannot be joined to a neighbor
+                if (myID == 0) {
+                    psArrayAdd (brightPeaks, 128, peak);
+                    continue;
+                }
+
+                // keep this peak if found->data.U8[myID] == false (no brighter peak in the footprint)
+                if (!found->data.U8[myID]) {
+                    // fprintf (stderr, "keeping %d: %d,%d\n", j, peak->x, peak->y);
+                    psArrayAdd (brightPeaks, 128, peak);
+                    continue;
+                }
+                // fprintf (stderr, "skipping %d: %d,%d\n", j, peak->x, peak->y);
+            }
+            psFree (myFP);
+            psFree (found);
+        }
+        psFree (threshist);
+        psFree (threshbounds);
+        psFree (threshbins);
+        psFree (peaktried);
+
+    } else {
+
+        // The brightest peak is always safe; go through other peaks trying to cull them
+        for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
+            const pmPeak *peak = fp->peaks->data[i];
+
+            float flux  = useSmoothedImage ? peak->smoothFlux      : peak->rawFlux;
+            float stdev = useSmoothedImage ? peak->smoothFluxStdev : peak->rawFluxStdev;
+
+            // if flux is negative, careful with fStdev
+            const float fStdev = fabs(stdev/flux);
+            const float stdev_pad = fabs(flux * hypot(fStdev, fPadding));
+
+            float threshold = flux - nsigma_delta*stdev_pad;
+
+            if (isnan(threshold)) {
+                // min_threshold is assumed to be below the detection threshold,
+                // so all the peaks are pmFootprint, and this isn't the brightest
+                continue;
+            }
+
+            // just in case, force the threshold below the peak source flux
+            if (threshold > flux) {
+                threshold = flux - 10*FLT_EPSILON;
+            }
+
+            if (threshold < min_threshold) {
+                threshold = min_threshold;
+            }
+            if (threshold > max_threshold) {
+                threshold = max_threshold;
+            }
+
+            pmFootprintsFindAtPoint(peakFootprint, peakFPSpans, subImg, subMask, threshold, brightPeaks, peak->y, peak->x);
+            // if (peakFPSpans->nStartSpans > 2000) {
+            //  // dumpfootprints(peakFootprint, peakFPSpans);
+            //  // fprintf (stderr, "big footprint %d : %d\n", peakFootprint->nspans, peakFPSpans->nStartSpans);
+            //  // fprintf (stderr, "big test footprint: %f %f to %f %f (%d pix)\n", peakFootprint->bbox.x0, peakFootprint->bbox.y0, peakFootprint->bbox.x1, peakFootprint->bbox.y1, peakFootprint->npix);
+            // }
+
+            // at this point brightPeaks only has the peaks brighter than the current
+
+            // we set the IDs to either 1 (in peak) or 0 (not in peak)
+            pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
+
+            // If this peak has not already been assigned to a source, then we can look for any
+            // brighter peaks within its footprint. Check if any of the previous (brighter) peaks
+            // are within the footprint of this peak If so, the current peak is bogus; drop it.
+            bool keep = true;
+            for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) {
+                // const pmPeak *peak2 = fp->peaks->data[j]; XXX isn't this an error?  we only care about the kept brighter peak, right?
+                const pmPeak *peak2 = brightPeaks->data[j];
+                int x2 = peak2->x - subImg->col0;
+                int y2 = peak2->y - subImg->row0;
+                if (idImg->data.S32[y2][x2] == IN_PEAK) {
+                    // There's a brighter peak within the footprint above threshold; so cull our initial peak
+                    keep = false;
+                }
+            }
+            if (!keep) {
+                psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak");  // we should not drop any already assigned peaks.
+                continue;
+            }
+            psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
+        }
     }
 
@@ -153 +324 @@
     psFree(idImg);
     psFree(subImg);
-    psFree(subWt);
     psFree(subMask);
     psFree(peakFootprint);

trunk/psModules/src/objects/pmFootprintFindAtPoint.c

@@ -31 +31 @@
  *
  * This is the guts of pmFootprintsFindAtPoint
+ * 
+ * This function is/was ill-defined if pixel values are NAN.  we should either treat NAN as >
+ * threshold or < threshold, but the current (r29004) code is ambiguous.
+ * EAM : change code so NAN is always > threshold
  */
 bool pmFootprintSpansBuild(pmFootprint *fp, // the footprint that we're building
@@ -103 +107 @@
         for (int j = x0 - 1; j >= -1; j--) {
             double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
-            if ((maskRow[j] & PM_STARTSPAN_DETECTED) || pixVal < threshold) {
+            bool belowThreshold = (pixVal < threshold) && isfinite(pixVal);
+            if ((maskRow[j] & PM_STARTSPAN_DETECTED) || belowThreshold)  {
                 if (j < x0 - 1) {       // we found some pixels above threshold
                     nx0 = j + 1;
@@ -119 +124 @@
             for (int j = nx0 + 1; j <= numCols; j++) {
                 double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
-                if ((maskRow[j] & PM_STARTSPAN_DETECTED) || pixVal < threshold) {
+                bool belowThreshold = (pixVal < threshold) && isfinite(pixVal);
+                if ((maskRow[j] & PM_STARTSPAN_DETECTED) || belowThreshold) {
                     nx1 = j - 1;
                     break;
@@ -146 +152 @@
         for (int j = nx1 + 1; j <= x1 + 1; j++) {
             double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
-            if (!(maskRow[j] & PM_STARTSPAN_DETECTED) && pixVal >= threshold) {
+            bool aboveThreshold = (pixVal >= threshold) || !isfinite(pixVal);
+            if (!(maskRow[j] & PM_STARTSPAN_DETECTED) && aboveThreshold) {
                 int sx0 = j++;          // span that we're working on is sx0:sx1
                 int sx1 = -1;           // We know that if we got here, we'll also set sx1
                 for (; j <= numCols; j++) {
                     double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
-                    if ((maskRow[j] & PM_STARTSPAN_DETECTED) || pixVal < threshold) { // end of span
+                    bool belowThreshold = (pixVal < threshold) && isfinite(pixVal);
+                    if ((maskRow[j] & PM_STARTSPAN_DETECTED) || belowThreshold) { // end of span
                         sx1 = j;
                         break;
@@ -306 +314 @@
         for (i = col; i >= 0; i--) {
             pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
-            if ((maskRow[i] & PM_STARTSPAN_DETECTED) || pixVal < threshold) {
+            bool belowThreshold = (pixVal < threshold) && isfinite(pixVal);
+            if ((maskRow[i] & PM_STARTSPAN_DETECTED) || belowThreshold) {
                 break;
             }
@@ -313 +322 @@
         for (i = col; i < numCols; i++) {
             pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
-            if ((maskRow[i] & PM_STARTSPAN_DETECTED) || pixVal < threshold) {
+            bool belowThreshold = (pixVal < threshold) && isfinite(pixVal);
+            if ((maskRow[i] & PM_STARTSPAN_DETECTED) || belowThreshold) {
                 break;
             }

trunk/psModules/src/objects/pmFootprintIDs.c

@@ -66 +66 @@
 
    if (footprints->n == 0) {
-       psError(PS_ERR_BAD_PARAMETER_SIZE, true, "You didn't provide any footprints");
+       // XXX this was an error, but is an allowed condition -- NULL returned image means no footprints for any peaks
+       // psError(PS_ERR_BAD_PARAMETER_SIZE, true, "You didn't provide any footprints");
        return NULL;
    }

trunk/psModules/src/objects/pmModel.c

@@ -41 +41 @@
     psFree(tmp->params);
     psFree(tmp->dparams);
+    psFree(tmp->covar);
     psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__);
 }
@@ -66 +67 @@
     tmp->chisqNorm = NAN;
     tmp->nDOF  = 0;
+    tmp->nPar  = 0;
     tmp->nPix  = 0;
     tmp->nIter = 0;
@@ -71 +73 @@
     tmp->flags = PM_MODEL_STATUS_NONE;
     tmp->residuals = NULL;              // do not free: the model does not own this memory
+    tmp->covar = NULL;
     tmp->isPCM = false;
 

trunk/psModules/src/objects/pmModel.h

@@ -33 +33 @@
     psVector *params;                   ///< Paramater values.
     psVector *dparams;                  ///< Parameter errors.
+    psImage *covar;                     ///< Optional covariance matrix
     float chisq;                        ///< Fit chi-squared.
     float chisqNorm;                    ///< re-normalized fit chi-squared.
@@ -38 +39 @@
     float magErr;                       ///< integrated model magnitude error
     int nPix;                           ///< number of pixels used for fit
-    int nDOF;                           ///< number of degrees of freedom
+    int nPar;                           ///< number of parameters in fit
+    int nDOF;                           ///< number of degrees of freedom (nDOF = nPix - nPar)
     int nIter;                          ///< number of iterations to reach min
     pmModelStatus flags;                ///< model status flags

trunk/psModules/src/objects/pmModelUtils.c

@@ -107 +107 @@
     return true;
 }
+
+bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments) {
+
+    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;
+
+    psEllipseShape shape = psEllipseAxesToShape (axes);
+
+    if (!isfinite(shape.sx))  return false;
+    if (!isfinite(shape.sy))  return false;
+    if (!isfinite(shape.sxy)) return false;
+
+    // set the shape parameters
+    *Sxx  = PS_MAX(0.5, M_SQRT2*shape.sx);
+    *Syy  = PS_MAX(0.5, M_SQRT2*shape.sy);
+    *Sxy  = shape.sxy;
+
+    return true;
+}
+
+bool pmModelSetNorm (float *Io, pmSource *source) {
+
+    *Io = source->peak->rawFlux;
+    if (!isfinite(*Io) && !source->moments) return false;
+
+    *Io = source->moments->Peak;
+    if (!isfinite(*Io)) return false;
+
+    return true;
+}
+
+bool pmModelSetPosition (float *Xo, float *Yo, pmSource *source) {
+
+    bool useMoments = pmSourcePositionUseMoments(source);
+    
+    if (useMoments) {
+        *Xo = source->moments->Mx;
+        *Yo = source->moments->My;
+    } else {
+        *Xo = source->peak->xf;
+        *Yo = source->peak->yf;
+    }
+
+    return true;
+}

trunk/psModules/src/objects/pmModelUtils.h

@@ -42 +42 @@
     );
 
-
+bool pmModelSetPosition (float *Xo, float *Yo, pmSource *source);
+bool pmModelSetNorm (float *Io, pmSource *source);
+bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments);
 
 /// @}

trunk/psModules/src/objects/pmMoments.h

@@ -51 +51 @@
     int nPixels;  ///< Number of pixels used.
 
-    float KronFlux;    ///< Kron flux (flux in 2.5*Mrf)
+    float KronCore;    ///< flux in r < 1.0*Mrf
+    float KronCoreErr;    ///< error on flux in r < 1.0*Mrf
+
+    float KronFlux;    ///< Kron flux (flux in r < 2.5*Mrf)
     float KronFluxErr; ///< Kron flux error
 
-    float KronFinner;    ///< Kron flux (flux in 2.5*Mrf)
-    float KronFouter;    ///< Kron flux (flux in 2.5*Mrf)
-
+    float KronFinner;    ///< Kron flux (flux in 1.0*Mrf < r < 2.5*Mrf)
+    float KronFouter;    ///< Kron flux (flux in 2.5*Mrf < r < 4.0*Mrf)
 }
 pmMoments;

trunk/psModules/src/objects/pmPCMdata.c

@@ -254 +254 @@
 
     pcm->nPix = nPix;
-    pcm->nDOF = nPix - nParams - 1;
+    pcm->nPar = nParams;
+    pcm->nDOF = nPix - nParams;
 
     return pcm;
@@ -341 +342 @@
         return false;
     }
-    pcm->nDOF = pcm->nPix - nParams - 1;
+    pcm->nPar = nParams;
+    pcm->nDOF = pcm->nPix - nParams;
 
     // has the source pixel window changed?

trunk/psModules/src/objects/pmPCMdata.h

@@ -33 +33 @@
     psMinConstraint *constraint;
     int nPix;
+    int nPar;
     int nDOF;
 } pmPCMdata;

trunk/psModules/src/objects/pmPSF.c

@@ -430 +430 @@
         return NAN;
     }
+
+    // get the model full-width at half-max
+    float fwhmMajor = 2*model->modelRadius (model->params, 0.5);
+
+# if (0)
     psF32 *params = model->params->data.F32; // Model parameters
     psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO); // Ellipse axes
@@ -439 +444 @@
 
     return fwhm;
-}
+# else
+
+    psFree(model);
+
+    return fwhmMajor;
+# endif
+}

trunk/psModules/src/objects/pmPSFtryFitPSF.c

@@ -109 +109 @@
 
         // This function calculates the psf and aperture magnitudes
-        status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal); // raw PSF mag, AP mag
+        status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag
         if (!status || isnan(source->apMag)) {
             psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask

trunk/psModules/src/objects/pmPeaks.c

@@ -115 +115 @@
 XXX: Macro this.
 *****************************************************************************/
+# if (0)
 static bool isItInThisRegion(const psRegion valid,
                              psS32 x,
@@ -130 +131 @@
     return(false);
 }
+# endif
 
 /******************************************************************************
@@ -148 +150 @@
     tmp->x = x;
     tmp->y = y;
-    tmp->value = value;
-    tmp->flux = value;
-    tmp->SN = 0;
+    tmp->detValue        = value;
+    tmp->rawFlux         = value; // set this by default: it is up to the user to supply a better value
+    tmp->rawFluxStdev    = NAN;
+    tmp->smoothFlux      = value; // set this by default: it is up to the user to supply a better value
+    tmp->smoothFluxStdev = NAN;
+    // tmp->SN = 0;
     tmp->xf = x;
     tmp->yf = y;
@@ -170 +175 @@
 
 
-// psSort comparison function for peaks
+// psSort comparison functions for peaks
 // XXX: Add error-checking for NULL args
-int pmPeaksCompareAscend (const void **a, const void **b)
-{
-    psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__);
+int pmPeaksSortByDetValueAscend (const void **a, const void **b)
+{
     pmPeak *A = *(pmPeak **)a;
     pmPeak *B = *(pmPeak **)b;
@@ -180 +184 @@
     psF32 diff;
 
-    diff = A->value - B->value;
+    diff = A->detValue - B->detValue;
     if (diff < FLT_EPSILON) {
-        psTrace("psModules.objects", 10, "---- %s(-1) end ----\n", __func__);
         return (-1);
     } else if (diff > FLT_EPSILON) {
-        psTrace("psModules.objects", 10, "---- %s(+1) end ----\n", __func__);
         return (+1);
     }
-    psTrace("psModules.objects", 10, "---- %s(0) end ----\n", __func__);
     return (0);
 }
-
-// psSort comparison function for peaks
-// XXX: Add error-checking for NULL args
-int pmPeaksCompareDescend (const void **a, const void **b)
-{
-    psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__);
+int pmPeaksSortByDetValueDescend (const void **a, const void **b)
+{
     pmPeak *A = *(pmPeak **)a;
     pmPeak *B = *(pmPeak **)b;
@@ -202 +199 @@
     psF32 diff;
 
-    diff = A->value - B->value;
+    diff = A->detValue - B->detValue;
     if (diff < FLT_EPSILON) {
-        psTrace("psModules.objects", 10, "---- %s(+1) end ----\n", __func__);
         return (+1);
     } else if (diff > FLT_EPSILON) {
-        psTrace("psModules.objects", 10, "---- %s(-1) end ----\n", __func__);
         return (-1);
     }
-    psTrace("psModules.objects", 10, "---- %s(0) end ----\n", __func__);
     return (0);
 }
-
-// sort by SN (descending)
-int pmPeakSortBySN (const void **a, const void **b)
+int pmPeaksSortByRawFluxAscend (const void **a, const void **b)
 {
     pmPeak *A = *(pmPeak **)a;
     pmPeak *B = *(pmPeak **)b;
 
-    psF32 fA = A->SN;
-    psF32 fB = B->SN;
-    if (isnan (fA)) fA = 0;
-    if (isnan (fB)) fB = 0;
-
-    psF32 diff = fA - fB;
-    if (diff > FLT_EPSILON) return (-1);
-    if (diff < FLT_EPSILON) return (+1);
+    psF32 diff;
+
+    diff = A->rawFlux - B->rawFlux;
+    if (diff < FLT_EPSILON) {
+        return (-1);
+    } else if (diff > FLT_EPSILON) {
+        return (+1);
+    }
     return (0);
 }
-
-// sort by Y (ascending)
-int pmPeakSortByY (const void **a, const void **b)
+int pmPeaksSortByRawFluxDescend (const void **a, const void **b)
 {
     pmPeak *A = *(pmPeak **)a;
     pmPeak *B = *(pmPeak **)b;
 
-    psF32 fA = A->y;
-    psF32 fB = B->y;
-
-    psF32 diff = fA - fB;
-    if (diff > FLT_EPSILON) return (+1);
-    if (diff < FLT_EPSILON) return (-1);
+    psF32 diff;
+
+    diff = A->rawFlux - B->rawFlux;
+    if (diff < FLT_EPSILON) {
+        return (+1);
+    } else if (diff > FLT_EPSILON) {
+        return (-1);
+    }
     return (0);
 }
+int pmPeaksSortBySmoothFluxAscend (const void **a, const void **b)
+{
+    pmPeak *A = *(pmPeak **)a;
+    pmPeak *B = *(pmPeak **)b;
+
+    psF32 diff;
+
+    diff = A->smoothFlux - B->smoothFlux;
+    if (diff < FLT_EPSILON) {
+        return (-1);
+    } else if (diff > FLT_EPSILON) {
+        return (+1);
+    }
+    return (0);
+}
+int pmPeaksSortBySmoothFluxDescend (const void **a, const void **b)
+{
+    pmPeak *A = *(pmPeak **)a;
+    pmPeak *B = *(pmPeak **)b;
+
+    psF32 diff;
+
+    diff = A->smoothFlux - B->smoothFlux;
+    if (diff < FLT_EPSILON) {
+        return (+1);
+    } else if (diff > FLT_EPSILON) {
+        return (-1);
+    }
+    return (0);
+}
+
+// // sort by SN (descending)
+// int pmPeakSortBySN (const void **a, const void **b)
+// {
+//     pmPeak *A = *(pmPeak **)a;
+//     pmPeak *B = *(pmPeak **)b;
+// 
+//     psF32 fA = A->flux;
+//     psF32 fB = B->flux;
+//     if (isnan (fA)) fA = 0;
+//     if (isnan (fB)) fB = 0;
+// 
+//     psF32 diff = fA - fB;
+//     if (diff > FLT_EPSILON) return (-1);
+//     if (diff < FLT_EPSILON) return (+1);
+//     return (0);
+// }
+// 
+// // sort by Y (ascending)
+// int pmPeakSortByY (const void **a, const void **b)
+// {
+//     pmPeak *A = *(pmPeak **)a;
+//     pmPeak *B = *(pmPeak **)b;
+// 
+//     psF32 fA = A->y;
+//     psF32 fB = B->y;
+// 
+//     psF32 diff = fA - fB;
+//     if (diff > FLT_EPSILON) return (+1);
+//     if (diff < FLT_EPSILON) return (-1);
+//     return (0);
+// }
 
 /******************************************************************************
@@ -554 +608 @@
     return(list);
 }
-
-// return a new array of peaks which are in the valid region and below threshold
-// XXX this function is unused and probably could be dropped
-psArray *pmPeaksSubset(
-    psArray *peaks,
-    psF32 maxValue,
-    const psRegion valid)
-{
-    psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(peaks, NULL);
-
-    psArray *output = psArrayAllocEmpty (200);
-
-    psTrace ("psModules.objects", 3, "list size is %ld\n", peaks->n);
-
-    for (int i = 0; i < peaks->n; i++) {
-        pmPeak *tmpPeak = (pmPeak *) peaks->data[i];
-        if (tmpPeak->value > maxValue)
-            continue;
-        if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))
-            continue;
-        psArrayAdd (output, 200, tmpPeak);
-    }
-    psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__);
-    return(output);
-}

trunk/psModules/src/objects/pmPeaks.h

@@ -35 +35 @@
     PM_PEAK_EDGE,                       ///< Peak on edge.
     PM_PEAK_FLAT,                       ///< Peak has equal-value neighbors.
+    PM_PEAK_SUSPECT_SATURATION,         ///< Peak is probably saturated
     PM_PEAK_UNDEF                       ///< Undefined.
 } pmPeakType;
@@ -45 +46 @@
  *  associated with the source:
  *
+ *  There are 3 values which define the amplitude of the peak and which may be used to sort the
+ *  peaks: 
+ *  * detValue - the peak in the detection image (nominally, the S/N)
+ *  * rawFlux - the peak in the unsmoothed image
+ *  * smoothFlux - the peak in the smoothed image
+ * 
+ *  For a given image, peaks do necesarily not have the same sequence for these three values.
+ *  Depending on the analysis, it may make sense to sort by one or the other of these values
  */
 typedef struct
@@ -55 +64 @@
     float dx;                           ///< bicube fit error on peak coord (x)
     float dy;                           ///< bicube fit error on peak coord (y)
-    float value;                        ///< level in detection image
-    float flux;                         ///< level in unsmoothed sci image
-    float SN;                           ///< S/N implied by detection level
+    float detValue;                     ///< peak flux in detection image (= S/N)
+    float rawFlux;                      ///< peak flux in unsmoothed signal image
+    float rawFluxStdev;                 ///< peak stdev in unsmoothed signal image
+    float smoothFlux;                   ///< peak flux in smoothed signal image
+    float smoothFluxStdev;              ///< peak stdev in smoothed signal image
+    // float SNestimated;                  ///< S/N estimated from the detection image
     bool assigned;                      ///< is peak assigned to a source?
     pmPeakType type;                    ///< Description of peak.
@@ -136 +148 @@
 );
 
-int pmPeaksCompareAscend (const void **a, const void **b);
-int pmPeaksCompareDescend (const void **a, const void **b);
-
-int pmPeakSortBySN (const void **a, const void **b);
-int pmPeakSortByY (const void **a, const void **b);
+int pmPeaksSortByDetValueAscend (const void **a, const void **b);
+int pmPeaksSortByDetValueDescend (const void **a, const void **b);
+int pmPeaksSortByRawFluxAscend (const void **a, const void **b);
+int pmPeaksSortByRawFluxDescend (const void **a, const void **b);
+int pmPeaksSortBySmoothFluxAscend (const void **a, const void **b);
+int pmPeaksSortBySmoothFluxDescend (const void **a, const void **b);
 
 /// @}

trunk/psModules/src/objects/pmPhotObj.c

@@ -67 +67 @@
         return false; 
     }
-    if (!finite(source->peak->xf)) {
+    if (!isfinite(source->peak->xf)) {
         psError(PS_ERR_UNKNOWN, true, "NAN peak coordinate");
         return false; 
     }
-    if (!finite(source->peak->yf)) {
+    if (!isfinite(source->peak->yf)) {
         psError(PS_ERR_UNKNOWN, true, "NAN peak coordinate");
         return false; 
@@ -81 +81 @@
         object->x  = source->peak->xf;
         object->y  = source->peak->yf;
-        object->SN = source->peak->SN;
+        object->flux = source->peak->rawFlux;
     } else {
-        object->SN = PS_MAX(object->SN, source->peak->SN);
+        object->flux = PS_MAX(object->flux, source->peak->rawFlux);
     }
     psArrayAdd (object->sources, 1, source);
@@ -89 +89 @@
 }
 
-// sort by SN (descending)
-int pmPhotObjSortBySN (const void **a, const void **b)
+// sort by flux (descending)
+int pmPhotObjSortByFlux (const void **a, const void **b)
 {
     pmPhotObj *objA = *(pmPhotObj **)a;
     pmPhotObj *objB = *(pmPhotObj **)b;
 
-    psF32 fA = objA->SN;
-    psF32 fB = objB->SN;
+    psF32 fA = objA->flux;
+    psF32 fB = objB->flux;
 
     psF32 diff = fA - fB;

trunk/psModules/src/objects/pmPhotObj.h

@@ -34 +34 @@
     float x;
     float y;
-    float SN;                           // max of peak->SN for all matched sources
+    float flux;                         // max of peak->rawFlux for all matched sources
 } pmPhotObj;
 
@@ -41 +41 @@
 bool pmPhotObjAddSource(pmPhotObj *object, pmSource *source);
 
-int pmPhotObjSortBySN (const void **a, const void **b);
+int pmPhotObjSortByFlux (const void **a, const void **b);
 int pmPhotObjSortByX (const void **a, const void **b);
 

trunk/psModules/src/objects/pmSource.c

@@ -171 +171 @@
     // peak has the same values as the original
     if (in->peak != NULL) {
-        source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->value, in->peak->type);
+        source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->detValue, in->peak->type);
         source->peak->xf = in->peak->xf;
         source->peak->yf = in->peak->yf;
-        source->peak->flux = in->peak->flux;
-        source->peak->SN = in->peak->SN;
+        source->peak->rawFlux         = in->peak->rawFlux;
+        source->peak->rawFluxStdev    = in->peak->rawFluxStdev;
+        source->peak->smoothFlux      = in->peak->smoothFlux;
+        source->peak->smoothFluxStdev = in->peak->smoothFluxStdev;
+        // source->peak->SN = in->peak->SN;
     }
 
@@ -328 +331 @@
 *****************************************************************************/
 
-// XXX EAM include a S/N cutoff in selecting the sources?
-// XXX this function should probably accept the values, not a recipe. wrap with a
-// psphot-specific function which applies the recipe values
-// only apply selection to sources within specified region
 pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 AR_MAX)
 {
@@ -429 +428 @@
 
         psfClump.nSigma = stats->sampleStdev;
+        psfClump.nTotal = nValid;
 
         if (savedImage) {
@@ -465 +465 @@
         psStats *stats  = NULL;
 
-        // select the single highest peak
-        psArraySort (peaks, pmPeaksCompareDescend);
+        // select the single highest peak (note that we only have detValue, not rawFlux, etc
+        psArraySort (peaks, pmPeaksSortByDetValueDescend);
         clump = peaks->data[0];
-        psTrace ("psModules.objects", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->value);
+        psTrace ("psModules.objects", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->detValue);
 
         // XXX store the mean sigma?
         float meanSigma = psfClump.nSigma;
-        psfClump.nStars = clump->value;
-        psfClump.nSigma = clump->value / meanSigma;
+        psfClump.nStars = clump->detValue;
+        psfClump.nSigma = clump->detValue / meanSigma;
 
         // define section window for clump
@@ -643 +643 @@
             }
 
+            // check for insignificant sources or excessively low-surface brightness
+            float coreSN = source->moments->KronCore / source->moments->KronCoreErr;
+            float coreKR = source->moments->KronCore / source->moments->KronFlux;
+
+            // XXX these values need to be in the recipe...
+            if (false && isfinite(coreSN) && (coreSN < 5.0)) {
+                source->type = PM_SOURCE_TYPE_DEFECT;
+                source->mode |= PM_SOURCE_MODE_DEFECT;
+                Ncr ++;
+                continue;
+            }
+            if (false && isfinite(coreKR) && (coreKR < 0.1)) {
+                source->type = PM_SOURCE_TYPE_DEFECT;
+                source->mode |= PM_SOURCE_MODE_DEFECT;
+                Ncr ++;
+                continue;
+            }
+
             // likely unsaturated extended source (too large to be stellar)
             if (sigX > clump.X + 3*clump.dX || sigY > clump.Y + 3*clump.dY) {
@@ -651 +669 @@
 
             // the rest are probable stellar objects
+            // the vectors below are accumulated to give user feedback on the S/N ranges
             starsn_moments->data.F32[starsn_moments->n] = source->moments->SN;
             starsn_moments->n ++;
-            starsn_peaks->data.F32[starsn_peaks->n] = source->peak->SN;
+            starsn_peaks->data.F32[starsn_peaks->n] = sqrt(source->peak->detValue);
             starsn_peaks->n ++;
             Nstar ++;
@@ -1149 +1168 @@
 }
 
+// this function decides if the source position should be based on the peak or the moments.
+// this is only used if we know we should not use a model fit position (eg, no model, or no
+// model yet)
+bool pmSourcePositionUseMoments(pmSource *source) {
+
+    if (!source->moments) return false;          // can't if there are no moments
+    if (!source->moments->nPixels) return false; // can't if the moments were not measured
+    if (source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE) return false; // can't if the moments failed...
+
+    if (source->mode & PM_SOURCE_MODE_SATSTAR) return true; // moments are best for SATSTARs
+
+    float dX = source->moments->Mx - source->peak->xf;
+    float dY = source->moments->My - source->peak->yf;
+    float dR = hypot(dX, dY);
+    
+    // only use the moments position if the moment-peak offset is small or the star is saturated
+    if (dR > 1.5) return false;
+
+    return true;
+}
+
 // sort by SN (descending)
-int pmSourceSortBySN (const void **a, const void **b)
+int pmSourceSortByFlux (const void **a, const void **b)
 {
     pmSource *A = *(pmSource **)a;
     pmSource *B = *(pmSource **)b;
 
-    psF32 fA = (A->peak == NULL) ? 0 : A->peak->SN;
-    psF32 fB = (B->peak == NULL) ? 0 : B->peak->SN;
+    psF32 fA = (A->peak == NULL) ? 0 : A->peak->rawFlux;
+    psF32 fB = (B->peak == NULL) ? 0 : B->peak->rawFlux;
     if (isnan (fA)) fA = 0;
     if (isnan (fB)) fB = 0;
@@ -1166 +1206 @@
 }
 
+// sort by SN (descending)
+int pmSourceSortByParentFlux (const void **a, const void **b)
+{
+    pmSource *Ao = *(pmSource **)a;
+    pmSource *Bo = *(pmSource **)b;
+    pmSource *A  = Ao->parent;
+    pmSource *B  = Bo->parent;
+
+    psF32 fA = (A->peak == NULL) ? 0 : A->peak->rawFlux;
+    psF32 fB = (B->peak == NULL) ? 0 : B->peak->rawFlux;
+    if (isnan (fA)) fA = 0;
+    if (isnan (fB)) fB = 0;
+
+    psF32 diff = fA - fB;
+    if (diff > FLT_EPSILON) return (-1);
+    if (diff < FLT_EPSILON) return (+1);
+    return (0);
+}
+
 // sort by Y (ascending)
 int pmSourceSortByY (const void **a, const void **b)
@@ -1201 +1260 @@
     pmSource *A = *(pmSource **)a;
     pmSource *B = *(pmSource **)b;
+
+    int iA = A->seq;
+    int iB = B->seq;
+
+    int diff = iA - iB;
+    if (diff > 0) return (+1);
+    if (diff < 0) return (-1);
+    return (0);
+}
+
+// sort by Seq (ascending)
+int pmSourceSortByParentSeq (const void **a, const void **b)
+{
+    pmSource *Ao = *(pmSource **)a;
+    pmSource *Bo = *(pmSource **)b;
+    pmSource *A  = Ao->parent;
+    pmSource *B  = Bo->parent;
 
     int iA = A->seq;

trunk/psModules/src/objects/pmSource.h

@@ -121 +121 @@
     float dY;
     int nStars;
+    int nTotal;
     float nSigma;
 }
@@ -286 +287 @@
 bool pmSourceCachePSF (pmSource *source, psImageMaskType maskVal);
 
-int  pmSourceSortBySN (const void **a, const void **b);
+bool pmSourcePositionUseMoments(pmSource *source);
+
 int  pmSourceSortByY (const void **a, const void **b);
 int  pmSourceSortByX (const void **a, const void **b);
 int  pmSourceSortBySeq (const void **a, const void **b);
+int  pmSourceSortByParentSeq (const void **a, const void **b);
+int  pmSourceSortByFlux (const void **a, const void **b);
+int  pmSourceSortByParentFlux (const void **a, const void **b);
 
 pmSourceMode pmSourceModeFromString (const char *name);

trunk/psModules/src/objects/pmSourceExtendedPars.c

@@ -55 +55 @@
     psFree(radial->flux);
     psFree(radial->fluxErr);
+    psFree(radial->fluxStdev);
     psFree(radial->fill);
 }
@@ -65 +66 @@
     radial->flux = NULL;
     radial->fluxErr = NULL;
+    radial->fluxStdev = NULL;
     radial->fill = NULL;
     return radial;

trunk/psModules/src/objects/pmSourceExtendedPars.h

@@ -23 +23 @@
 typedef struct {
     psVector *flux;                     // fluxes measured at above radii
+    psVector *fluxStdev;                // scatter (standard deviation) of flux
     psVector *fluxErr;                  // formal error on the fluxes (sqrt\sum(variance))
     psVector *fill;                     // angles corresponding to above radial profiles
@@ -63 +64 @@
     float petrosianR50;
     float petrosianR50Err;
+    float petrosianFill;
 } pmSourceExtendedPars;
 

trunk/psModules/src/objects/pmSourceFitModel.c

@@ -40 +40 @@
 #include "pmSourceDiffStats.h"
 #include "pmSource.h"
+#include "pmSourcePhotometry.h"
 #include "pmSourceFitModel.h"
 
@@ -59 +60 @@
     opt->maxChisqDOF = NAN;
     opt->poissonErrors = true;
+    opt->saveCovariance = false;
 
     return opt;
@@ -231 +233 @@
         psTrace ("psModules.objects", 4, "%f +/- %f", params->data.F32[i], dparams->data.F32[i]);
     }
+    if (options->saveCovariance) {
+        model->covar = psMemIncrRefCounter(covar);
+    }
+    model->nIter = myMin->iter;
+    model->nPar = nParams;
+
     psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
 
     // save the resulting chisq, nDOF, nIter
-    model->chisq = myMin->value;
-    model->nIter = myMin->iter;
-    model->nPix  = y->n;
-    model->nDOF  = y->n - nParams;
-    model->chisqNorm = model->chisq / model->nDOF;
+    // NOTE: if (!options->poissonErrors) chisq will be wrong : recalculate
+    if (options->poissonErrors) {
+        model->chisq = myMin->value;
+        model->nPix  = y->n;
+        model->nDOF  = y->n - model->nPar;
+        model->chisqNorm = model->chisq / model->nDOF;
+    } else {
+        pmSourceChisqUnsubtracted (source, model, maskVal);
+    }
+
+    // set the model success or failure status
     model->flags |= PM_MODEL_STATUS_FITTED;
     if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;

trunk/psModules/src/objects/pmSourceFitModel.h

@@ -31 +31 @@
     float maxChisqDOF;                  ///< convergence criterion
     float weight;                       ///< use this weight for constant-weight fits
+    float covarFactor;                  ///< covariance factor for calculating the chisq
     bool poissonErrors;                 ///< use poisson errors for fits?
+    bool saveCovariance;
 } pmSourceFitOptions;
 

trunk/psModules/src/objects/pmSourceFitPCM.c

@@ -38 +38 @@
 #include "pmSourceDiffStats.h"
 #include "pmSource.h"
+#include "pmSourcePhotometry.h"
 #include "pmSourceFitModel.h"
 #include "pmPCMdata.h"
@@ -70 +71 @@
         psTrace ("psModules.objects", 4, "%f +/- %f", params->data.F32[i], dparams->data.F32[i]);
     }
+    if (fitOptions->saveCovariance) {
+        psFree(pcm->modelConv->covar);
+        pcm->modelConv->covar = psMemIncrRefCounter(covar);
+    }
     psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
 
@@ -79 +84 @@
         }
     }
+    pcm->modelConv->nIter = myMin->iter;
+    pcm->modelConv->nPar = pcm->nPar;
 
     // save the resulting chisq, nDOF, nIter
-    pcm->modelConv->chisq = myMin->value;
-    pcm->modelConv->nIter = myMin->iter;
-    pcm->modelConv->nPix = pcm->nPix;
-    pcm->modelConv->nDOF = pcm->nDOF;
-    pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF;
+    if (fitOptions->poissonErrors) {
+        pcm->modelConv->chisq = myMin->value;
+        pcm->modelConv->nPix = pcm->nPix;
+        pcm->modelConv->nDOF = pcm->nDOF;
+        pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF;
+    } else {
+        // xxx this is wrong because it does not convolve with the psf
+        pmSourceChisqUnsubtracted (source, pcm->modelConv, maskVal);
+    }
+
+    // set the model success or failure status
     pcm->modelConv->flags |= PM_MODEL_STATUS_FITTED;
     if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE;

trunk/psModules/src/objects/pmSourceFitSet.c

@@ -39 +39 @@
 #include "pmSourceDiffStats.h"
 #include "pmSource.h"
+#include "pmSourcePhotometry.h"
 
 #include "pmSourceFitModel.h"
@@ -296 +297 @@
 
 // set the model parameters for this fit set
-bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam,
-                           const psVector *param, pmSource *source,
-                           psMinimization *myMin, int nPix, bool fitStatus)
+bool pmSourceFitSetValues (pmSourceFitSetData *set, 
+                           const psVector *dparam, const psVector *param, const psImage *covar, 
+                           pmSource *source, psMinimization *myMin, int nPix, 
+                           bool fitStatus, pmSourceFitOptions *options, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(set, false);
@@ -311 +313 @@
 
     int n = 0;
+    int nStart = 0;
     for (int i = 0; i < set->paramSet->n; i++) {
 
@@ -320 +323 @@
             psTrace ("psModules.objects", 4, "%f +/- %f", param->data.F32[n], dparam->data.F32[n]);
         }
+        if (options->saveCovariance) {
+            // we only save the covar matrix for this object with itself (ignore cross terms between objects)
+            model->covar = psImageAlloc(model->params->n, model->params->n, PS_TYPE_F32);
+            for (int ix = 0; ix < model->params->n; ix++) {
+                for (int iy = 0; iy < model->params->n; iy++) {
+                    model->covar->data.F32[iy][ix] = covar->data.F32[nStart+iy][nStart+ix];
+                }
+            }
+        }
+        nStart += model->params->n;
         psTrace ("psModules.objects", 4, " src %d", i);
+
+        model->nIter = myMin->iter;
+        // model->nPar is set by pmSourceFitSetMasks
 
         // save the resulting chisq, nDOF, nIter
         // these are not unique for any one source
-        model->chisq = myMin->value;
-        model->nIter = myMin->iter;
-        model->nDOF  = nPix - model->params->n;
+        if (options->poissonErrors) {
+            model->chisq = myMin->value;
+            model->nPix  = nPix;
+            model->nDOF  = nPix - model->nPar;
+            model->chisqNorm = model->chisq / model->nDOF;
+        } else {
+            pmSourceChisqUnsubtracted (source, model, maskVal);
+        }
 
         // set the model success or failure status
@@ -380 +401 @@
     for (int i = 0; i < set->paramSet->n; i++) {
         psVector *paramOne = set->paramSet->data[i];
+        pmModel  *modelOne = set->modelSet->data[i];
 
         switch (mode) {
@@ -387 +409 @@
                 if (j == PM_PAR_I0) continue;
                 constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + j] = 1;
+                modelOne->nPar = 1;
             }
             break;
@@ -396 +419 @@
                 if (j == PM_PAR_I0) continue;
                 constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + j] = 1;
+                modelOne->nPar = 3;
             }
             break;
@@ -401 +425 @@
             // EXT model fits all params (except sky)
             constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + PM_PAR_SKY] = 1;
+            modelOne->nPar = paramOne->n - 1;
             break;
           default:
@@ -550 +575 @@
     }
 
-// parameter errors from the covariance matrix
+    // parameter errors from the covariance matrix
     psVector *dparams = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32);
     for (int i = 0; i < dparams->n; i++) {
@@ -558 +583 @@
     }
 
-// get the Gauss-Newton distance for fixed model parameters
+    // get the Gauss-Newton distance for fixed model parameters
     if (constraint->paramMask != NULL) {
         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
@@ -580 +605 @@
     }
 
-    pmSourceFitSetValues (thisSet, dparams, params, source, myMin, y->n, fitStatus);
+    pmSourceFitSetValues (thisSet, dparams, params, covar, source, myMin, y->n, fitStatus, options, maskVal);
     psTrace ("psModules.objects", 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nPix: %ld\n", onPic, fitStatus, myMin->iter, myMin->value, y->n);
 

trunk/psModules/src/objects/pmSourceFitSet.h

@@ -40 +40 @@
 bool pmSourceFitSetJoin (psVector *deriv, psVector *param, pmSourceFitSetData *set);
 bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param);
-bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, const psVector *param, pmSource *source, psMinimization *myMin, int nPix, bool fitStatus);
+
+bool pmSourceFitSetValues (pmSourceFitSetData *set, 
+                           const psVector *dparam, const psVector *param, const psImage *covar, 
+                           pmSource *source, psMinimization *myMin, int nPix, 
+                           bool fitStatus, pmSourceFitOptions *options, psImageMaskType maskVal);
 
 psF32 pmSourceFitSetFunction(psVector *deriv, const psVector *param, const psVector *x);
@@ -57 +61 @@
     psArray *modelSet,                  ///< model to be fitted
     pmSourceFitOptions *options,        ///< define options for fitting process
-    psImageMaskType maskVal             ///< Vale to mask
+    psImageMaskType maskVal             ///< Value to mask
 
 );

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c

@@ -49 +49 @@
 
 #include "pmSourceIO.h"
+#include "pmSourceOutputs.h"
 
 // panstarrs-style FITS table output (header + table in 1st extension)
@@ -62 +63 @@
     PS_ASSERT_PTR_NON_NULL(extname, false);
 
-    int i;
     psArray *table;
     psMetadata *row;
-    psF32 *PAR, *dPAR;
-    psEllipseAxes axes;
-    psF32 xPos, yPos;
-    psF32 xErr, yErr;
-    psF32 errMag, chisq, apRadius;
-    psS32 nPix, nDOF;
 
     pmChip *chip = readout->parent->parent;
-    pmFPA  *fpa  = chip->parent;
-
-    bool status1 = false;
-    bool status2 = false;
-    float magOffset = NAN;
-    float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
-    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
-    if (!isfinite(zeropt)) {
-        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
-    }
-    if (status1 && status2 && (exptime > 0.0)) {
-        magOffset = zeropt + 2.5*log10(exptime);
-    }
-    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
+
 
     // if the sequence is defined, write these in seq order; otherwise
@@ -94 +75 @@
         if (source->seq == -1) {
           // let's write these out in S/N order
-          sources = psArraySort (sources, pmSourceSortBySN);
+          sources = psArraySort (sources, pmSourceSortByFlux);
         } else {
           sources = psArraySort (sources, pmSourceSortBySeq);
@@ -100 +81 @@
     }
 
+    short nImageOverlap; 
+    float magOffset; 
+    float zeroptErr; 
+    float fwhmMajor; 
+    float fwhmMinor;
+    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
+
     table = psArrayAllocEmpty (sources->n);
 
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
     // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
-    for (i = 0; i < sources->n; i++) {
+    for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
 
@@ -115 +103 @@
         }
 
-        // no difference between PSF and non-PSF model
-        pmModel *model = source->modelPSF;
-
-        if (model != NULL) {
-            PAR = model->params->data.F32;
-            dPAR = model->dparams->data.F32;
-            xPos = PAR[PM_PAR_XPOS];
-            yPos = PAR[PM_PAR_YPOS];
-            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
-              xErr = dPAR[PM_PAR_XPOS];
-              yErr = dPAR[PM_PAR_YPOS];
-            } else {
-              // in linear-fit mode, there is no error on the centroid
-              xErr = source->peak->dx;
-              yErr = source->peak->dy;
-            }
-            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0);
-            } else {
-                axes.major = NAN;
-                axes.minor = NAN;
-                axes.theta = NAN;
-            }
-            chisq = model->chisq;
-            nDOF = model->nDOF;
-            nPix = model->nPix;
-            apRadius = source->apRadius;
-            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
-        } else {
-            xPos = source->peak->xf;
-            yPos = source->peak->yf;
-            xErr = source->peak->dx;
-            yErr = source->peak->dy;
-            axes.major = NAN;
-            axes.minor = NAN;
-            axes.theta = NAN;
-            chisq = NAN;
-            nDOF = 0;
-            nPix = 0;
-            apRadius = NAN;
-            errMag = NAN;
-        }
-
-        float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
-        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
-        psS16 nImageOverlap = 1;
-
-        psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
-        float posAngle = 0.0;
-        float pltScale = 0.0;
-        pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
+        // set the 'best' values for various output fields:
+        pmSourceOutputs outputs;
+        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
+
+        pmSourceOutputsMoments moments;
+        pmSourceOutputsSetMoments (&moments, source);
 
         pmSourceDiffStats diffStats;
@@ -175 +118 @@
         row = psMetadataAlloc ();
         psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr); // XXX this is only measured for non-linear fits
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr); // XXX this is only measured for non-linear fits
-        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr); // XXX this is only measured for non-linear fits
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr); // XXX this is only measured for non-linear fits
+        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
+        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfFlux);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfFluxErr);
         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
-        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
+        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
-        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
+        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
+        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
         psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
         psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
-
-        // distinguish moments measure from window vs S/N > XX ??
-        float mxx = source->moments ? source->moments->Mxx : NAN;
-        float mxy = source->moments ? source->moments->Mxy : NAN;
-        float myy = source->moments ? source->moments->Myy : NAN;
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
-
-        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                    diffStats.nGood);
-        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                      diffStats.fRatio);
-        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                      diffStats.nRatioBad);
-        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                     diffStats.nRatioMask);
-        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",             diffStats.nRatioAll);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                     diffStats.nGood);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                       diffStats.fRatio);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                       diffStats.nRatioBad);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                      diffStats.nRatioMask);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",              diffStats.nRatioAll);
 
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
@@ -331 +270 @@
         // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
         source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
-        source->peak->flux = peakFlux;
+        source->peak->rawFlux = peakFlux;
+        source->peak->smoothFlux = peakFlux;
         source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
         source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];
-        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
-          source->peak->SN = 1.0 / source->errMag;
-        } else {
-          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
-        }
 
         source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
@@ -353 +288 @@
 
         source->moments = pmMomentsAlloc ();
+        source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
+        source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
+
         source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
         source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
@@ -395 +333 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -568 +506 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c

@@ -49 +49 @@
 
 #include "pmSourceIO.h"
+#include "pmSourceOutputs.h"
 
 // panstarrs-style FITS table output (header + table in 1st extension)
@@ -62 +63 @@
     PS_ASSERT_PTR_NON_NULL(extname, false);
 
-    int i;
     psArray *table;
     psMetadata *row;
-    psF32 *PAR, *dPAR;
-    psEllipseAxes axes;
-    psF32 xPos, yPos;
-    psF32 xErr, yErr;
-    psF32 chisq, apRadius;
-    psS32 nPix, nDOF;
 
     pmChip *chip = readout->parent->parent;
-    pmFPA  *fpa  = chip->parent;
-
-    bool status1 = false;
-    bool status2 = false;
-    float magOffset = NAN;
-    float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
-    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
-    if (!isfinite(zeropt)) {
-        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
-    }
-    if (status1 && status2 && (exptime > 0.0)) {
-        magOffset = zeropt + 2.5*log10(exptime);
-    }
-    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
-
-    // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors
-    float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ");
-    if (!status1) {
-        fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1");
-        if (!status1) {
-            fwhmMajor = 5.0; // XXX just a guess!
-        }
-    }
-    float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN");
-    if (!status1) {
-        fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2");
-        if (!status1) {
-            fwhmMinor = 5.0; // XXX just a guess!
-        }
-    }
 
     // if the sequence is defined, write these in seq order; otherwise
@@ -110 +74 @@
         if (source->seq == -1) {
           // let's write these out in S/N order
-          sources = psArraySort (sources, pmSourceSortBySN);
+          sources = psArraySort (sources, pmSourceSortByFlux);
         } else {
           sources = psArraySort (sources, pmSourceSortBySeq);
@@ -118 +82 @@
     table = psArrayAllocEmpty (sources->n);
 
+    short nImageOverlap; 
+    float magOffset; 
+    float zeroptErr; 
+    float fwhmMajor; 
+    float fwhmMinor;
+    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
+
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
     // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
-    for (i = 0; i < sources->n; i++) {
+    for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
 
@@ -131 +102 @@
         }
 
-        // no difference between PSF and non-PSF model
-        pmModel *model = source->modelPSF;
-
-        if (model != NULL) {
-            PAR = model->params->data.F32;
-            dPAR = model->dparams->data.F32;
-            xPos = PAR[PM_PAR_XPOS];
-            yPos = PAR[PM_PAR_YPOS];
-            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
-              xErr = dPAR[PM_PAR_XPOS];
-              yErr = dPAR[PM_PAR_YPOS];
-            } else {
-              // in linear-fit mode, there is no error on the centroid
-                xErr = fwhmMajor * source->errMag / 2.35;
-                yErr = fwhmMinor * source->errMag / 2.35;
-            }
-            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0);
-            } else {
-                axes.major = NAN;
-                axes.minor = NAN;
-                axes.theta = NAN;
-            }
-            chisq = model->chisq;
-            nDOF = model->nDOF;
-            nPix = model->nPix;
-            apRadius = source->apRadius;
-        } else {
-            bool useMoments = true;
-            useMoments = (useMoments && source->moments);          // can't if there are no moments
-            useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
-            useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
-
-            if (useMoments) {
-                xPos = source->moments->Mx;
-                yPos = source->moments->My;
-                xErr = fwhmMajor * source->errMag / 2.35;
-                yErr = fwhmMinor * source->errMag / 2.35;
-            } else {
-                xPos = source->peak->xf;
-                yPos = source->peak->yf;
-                xErr = source->peak->dx;
-                yErr = source->peak->dy;
-            }
-            axes.major = NAN;
-            axes.minor = NAN;
-            axes.theta = NAN;
-            chisq = NAN;
-            nDOF = 0;
-            nPix = 0;
-            apRadius = NAN;
-        }
-
-        float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
-        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
-        psS16 nImageOverlap = 1;
-
-        psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
-        float posAngle = 0.0;
-        float pltScale = 0.0;
-        pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
+        // set the 'best' values for various output fields:
+        pmSourceOutputs outputs;
+        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
+
+        pmSourceOutputsMoments moments;
+        pmSourceOutputsSetMoments (&moments, source);
 
         pmSourceDiffStats diffStats;
@@ -201 +117 @@
         row = psMetadataAlloc ();
         psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr); // XXX this is only measured for non-linear fits
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr); // XXX this is only measured for non-linear fits
-        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr); // XXX this is only measured for non-linear fits
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr); // XXX this is only measured for non-linear fits
+        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
+        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
@@ -214 +130 @@
         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in real aperture",                 source->apMagRaw);
-        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
+        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
         psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",     source->apFlux);
         psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",                        source->apFluxErr);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
-        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
+        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
+        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
         psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
         psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
-
-        // distinguish moments measure from window vs S/N > XX ??
-        float mxx = source->moments ? source->moments->Mxx : NAN;
-        float mxy = source->moments ? source->moments->Mxy : NAN;
-        float myy = source->moments ? source->moments->Myy : NAN;
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
-
-        float Mrf  = source->moments ? source->moments->Mrf : NAN;
-        float Mrh  = source->moments ? source->moments->Mrh : NAN;
-        float Krf  = source->moments ? source->moments->KronFlux : NAN;
-        float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
-
-        float Kinner = source->moments ? source->moments->KronFinner : NAN;
-        float Kouter = source->moments ? source->moments->KronFouter : NAN;
-
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
-
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",                     Kinner);
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",                     Kouter);
-
-        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                    diffStats.nGood);
-        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                      diffStats.fRatio);
-        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                      diffStats.nRatioBad);
-        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                     diffStats.nRatioMask);
-        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",             diffStats.nRatioAll);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",                      moments.Kinner);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",                      moments.Kouter);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                     diffStats.nGood);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                       diffStats.fRatio);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                       diffStats.nRatioBad);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                      diffStats.nRatioMask);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",              diffStats.nRatioAll);
 
         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",          diffStats.Rp);
@@ -385 +288 @@
         // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
         source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
-        source->peak->flux = peakFlux;
+        source->peak->rawFlux = peakFlux;
+        source->peak->smoothFlux = peakFlux;
         source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
         source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];
-        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
-          source->peak->SN = 1.0 / source->errMag;
-        } else {
-          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
-        }
 
         source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
@@ -408 +307 @@
 
         source->moments = pmMomentsAlloc ();
+        source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
+        source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
+
         source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
         source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
@@ -456 +358 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -629 +531 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

@@ -50 +50 @@
 
 #include "pmSourceIO.h"
+#include "pmSourceOutputs.h"
 
 // panstarrs-style FITS table output (header + table in 1st extension)
@@ -66 +67 @@
     psArray *table;
     psMetadata *row;
-    psF32 *PAR, *dPAR;
-    psEllipseAxes axes;
-    psF32 xPos, yPos;
-    psF32 xErr, yErr;
-    psF32 errMag, chisq, apRadius;
-    psS32 nPix, nDOF;
 
     pmChip *chip = readout->parent->parent;
-    pmFPA  *fpa  = chip->parent;
-
-    bool status1 = false;
-    bool status2 = false;
-    float magOffset = NAN;
-    float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
-    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
-    if (!isfinite(zeropt)) {
-        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
-    }
-    if (status1 && status2 && (exptime > 0.0)) {
-        magOffset = zeropt + 2.5*log10(exptime);
-    }
-    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
-
-    // if the sequence is defined, write these in seq order; otherwise
-    // write them in S/N order:
+
+    // if the sequence is defined, write these in seq order; otherwise write them in S/N order.
+    // Careful: if we are working with child sources, then we need to sort by the parent info,
+    // not our info
     if (sources->n > 0) {
-        pmSource *source = (pmSource *) sources->data[0];
+        pmSource *source = sources->data[0];
         if (source->seq == -1) {
-            // let's write these out in S/N order
-            sources = psArraySort (sources, pmSourceSortBySN);
+            sources = psArraySort (sources, pmSourceSortByFlux);
         } else {
             sources = psArraySort (sources, pmSourceSortBySeq);
@@ -103 +84 @@
     table = psArrayAllocEmpty (sources->n);
 
-# if (0)
-    // we use this just to define the output vectors (which must be present for all objects)
-    bool status = false;
-    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
-    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
-    psAssert (radMax, "this must have been defined and tested earlier!");
-    psAssert (radMax->n, "this must have been defined and tested earlier!");
-    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
-
-    // write the radial profile apertures to header
-    for (int i = 0; i < radMax->n; i++) {
-      sprintf (keyword1, "RMIN_%02d", i);
-      sprintf (keyword2, "RMAX_%02d", i);
-      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
-      psMetadataAddF32 (imageHeader, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
-    }
-# endif
+    short nImageOverlap; 
+    float magOffset; 
+    float zeroptErr; 
+    float fwhmMajor; 
+    float fwhmMinor;
+    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
 
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
     // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
     for (int i = 0; i < sources->n; i++) {
-        pmSource *source = (pmSource *) sources->data[i];
-
-        // If source->seq is -1, source was generated in this analysis.  If source->seq is
-        // not -1, source was read from elsewhere: in the latter case, preserve the source
-        // ID.  source.seq is used instead of source.id since the latter is a const
-        // generated on Alloc, and would thus be wrong for read in sources.
+        // this is the source associated with this image
+        pmSource *thisSource = sources->data[i];
+
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+
+        // If source->seq is -1, source is unique and generated in this analysis.  If
+        // source->seq is not -1, source was read from elsewhere or tied to other source (eg
+        // from another image): in the latter case, preserve the source ID.  source.seq is used
+        // instead of source.id since the latter is a const generated on Alloc, and would thus
+        // be wrong for read in sources.
         if (source->seq == -1) {
             source->seq = i;
         }
 
-        // no difference between PSF and non-PSF model
-        pmModel *model = source->modelPSF;
-
-        if (model != NULL) {
-            PAR = model->params->data.F32;
-            dPAR = model->dparams->data.F32;
-            xPos = PAR[PM_PAR_XPOS];
-            yPos = PAR[PM_PAR_YPOS];
-            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
-                xErr = dPAR[PM_PAR_XPOS];
-                yErr = dPAR[PM_PAR_YPOS];
-            } else {
-                // in linear-fit mode, there is no error on the centroid
-                xErr = source->peak->dx;
-                yErr = source->peak->dy;
-            }
-            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0);
-            } else {
-                axes.major = NAN;
-                axes.minor = NAN;
-                axes.theta = NAN;
-            }
-            chisq = model->chisq;
-            nDOF = model->nDOF;
-            nPix = model->nPix;
-            apRadius = source->apRadius;
-            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
-        } else {
-            xPos = source->peak->xf;
-            yPos = source->peak->yf;
-            xErr = source->peak->dx;
-            yErr = source->peak->dy;
-            axes.major = NAN;
-            axes.minor = NAN;
-            axes.theta = NAN;
-            chisq = NAN;
-            nDOF = 0;
-            nPix = 0;
-            apRadius = NAN;
-            errMag = NAN;
-        }
-
-        float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
-        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
-        psS16 nImageOverlap = 1;
-
-        psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
-        float posAngle = 0.0;
-        float pltScale = 0.0;
-        pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
+        // set the 'best' values for various output fields:
+        pmSourceOutputs outputs;
+        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
+
+        pmSourceOutputsMoments moments;
+        pmSourceOutputsSetMoments (&moments, source);
 
         row = psMetadataAlloc ();
         psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr); // XXX this is only measured for non-linear fits
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr); // XXX this is only measured for non-linear fits
-        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
+        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
+        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
-        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
-        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
+        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
-        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
+        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
+        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
         psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
         psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
-
-        // distinguish moments measure from window vs S/N > XX ??
-        float Mxx = source->moments ? source->moments->Mxx : NAN;
-        float Mxy = source->moments ? source->moments->Mxy : NAN;
-        float Myy = source->moments ? source->moments->Myy : NAN;
-
-        float Mrf  = source->moments ? source->moments->Mrf : NAN;
-        float Mrh  = source->moments ? source->moments->Mrh : NAN;
-        float Krf  = source->moments ? source->moments->KronFlux : NAN;
-        float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
-
-        float Kinner = source->moments ? source->moments->KronFinner : NAN;
-        float Kouter = source->moments ? source->moments->KronFouter : NAN;
-
-        float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN;
-        float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN;
-        float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN;
-        float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN;
-
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      Mxx);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      Mxy);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      Myy);
-
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                    M_c3);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                    M_s3);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                   M_c4);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                   M_s4);
-
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
-
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kinner);
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kouter);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
 
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                     source->mode);
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
-
-# if (0)
-        // XXX if we have raw radial apertures, write them out here
-        psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
-        psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
-        psVector *radFill    = psVectorAlloc(radMax->n, PS_TYPE_F32);
-        psVectorInit (radFlux,    NAN);
-        psVectorInit (radFluxErr, NAN);
-        psVectorInit (radFill,    NAN);
-        if (!source->radial) goto empty_annuli;
-        if (!source->radial->flux) goto empty_annuli;
-        if (!source->radial->fill) goto empty_annuli;
-        psAssert (source->radial->flux->n <= radFlux->n, "inconsistent vector lengths");
-        psAssert (source->radial->fill->n <= radFlux->n, "inconsistent vector lengths");
-
-        // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
-        for (int j = 0; j < source->radial->flux->n; j++) {
-            radFlux->data.F32[j]    = source->radial->flux->data.F32[j];
-            radFluxErr->data.F32[j] = source->radial->fluxErr->data.F32[j];
-            radFill->data.F32[j]    = source->radial->fill->data.F32[j];
-        }
-
-    empty_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_FILL",     PS_DATA_VECTOR, "fill factor of annuli", radFill);
-        psFree (radFlux);
-        psFree (radFluxErr);
-        psFree (radFill);
-# endif
 
         // XXX not sure how to get this : need to load Nimages with weight?
@@ -406 +287 @@
         // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
         source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
-        source->peak->flux = peakFlux;
+        source->peak->rawFlux = peakFlux;
+        source->peak->smoothFlux = peakFlux;
         source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
         source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];
-        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
-          source->peak->SN = 1.0 / source->errMag;
-        } else {
-          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
-        }
 
         source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
@@ -429 +306 @@
 
         source->moments = pmMomentsAlloc ();
+        source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
+        source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
+
         source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
         source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
@@ -500 +380 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -522 +402 @@
     // we write out all sources, regardless of quality.  the source flags tell us the state
     for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
+        // this is the source associated with this image
+        pmSource *thisSource = sources->data[i];
+
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
 
         // skip sources without measurements
@@ -556 +440 @@
         }
         psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO",       PS_DATA_F32, "Axial Ratio of radial profile",              AxialRatio);
-        psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",                  AxialTheta);
+        psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",            AxialTheta);
 
         // Petrosian measurements
@@ -567 +451 @@
                 float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
                 float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude", mag);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", magErr);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",           PS_DATA_F32, "Petrosian Magnitude", mag);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",       PS_DATA_F32, "Petrosian Magnitude Error", magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",        PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR",    PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr);
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50);
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err);
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90);
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL",          PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill);
             } else {
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
-                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",           PS_DATA_F32, "Petrosian Magnitude",       NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",       PS_DATA_F32, "Petrosian Magnitude Error", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",        PS_DATA_F32, "Petrosian Radius",          NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR",    PS_DATA_F32, "Petrosian Radius Error",    NAN);
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", NAN);
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); 
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", NAN);
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); 
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL",          PS_DATA_F32, "Petrosian Fill Factor", NAN);
             }
         }
@@ -672 +558 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
     int nParamMax = 0;
     for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
+        // this is the source associated with this image
+        pmSource *thisSource = sources->data[i];
+
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+
         if (source->modelFits == NULL) continue;
         for (int j = 0; j < source->modelFits->n; j++) {
@@ -691 +582 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
+        pmSource *thisSource = sources->data[i];
+
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
 
         // XXX if no model fits are saved, write out modelEXT?
@@ -747 +641 @@
 
                 if (k < model->params->n) {
-                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]);
                 } else {
-                    psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", NAN);
                 }
             }
 
-            // XXX other parameters which may be set.
-            // XXX flag / value to define the model
-            // XXX write out the model type, fit status flags
-
+            // optionally, write out the covariance matrix values
+            // XXX do I need to pad this to match the biggest covar matrix?
+            if (model->covar) {
+                for (int iy = 0; iy < model->covar->numCols; iy++) {
+                    for (int ix = iy; ix < model->covar->numCols; ix++) {
+                        snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix);
+                        psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]);
+
+                    }
+                }                   
+            }
             psArrayAdd (table, 100, row);
             psFree (row);
@@ -790 +691 @@
 bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
 {
-
+    bool status = false;
     psArray *table;
     psMetadata *row;
@@ -796 +697 @@
     char keyword1[80], keyword2[80];
 
+    // perform full non-linear fits / extended source analysis?
+    if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
+        psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n");
+        return true;
+    }
+
     // create a header to hold the output data
     psMetadata *outhead = psMetadataAlloc ();
@@ -803 +710 @@
 
     // we use this just to define the output vectors (which must be present for all objects)
-    bool status = false;
     psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
     psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
@@ -818 +724 @@
     }
 
+    // the FWHM values are available if we measured a psf-matched convolved set
     psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
-    if (!fwhmValues) {
-        psError (PM_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
-        return false;
-    }
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -832 +735 @@
     for (int i = 0; i < sources->n; i++) {
 
-        pmSource *source = sources->data[i];
+        // this is the source associated with this image
+        pmSource *thisSource = sources->data[i];
+
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
 
         // skip sources without radial aper measurements (or insufficient)
-        if (source->radialAper == NULL) continue;
-        psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
-
-        for (int entry = 0; entry < fwhmValues->n; entry++) {
+        if (source->radialAper == NULL) continue;
+
+        // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
+
+        for (int entry = 0; entry < source->radialAper->n; entry++) {
 
             // choose the convolved EXT model, if available, otherwise the simple one
@@ -844 +752 @@
             assert (radialAper);
 
-            bool useMoments = true;
-            useMoments = (useMoments && source->moments);          // can't if there are no moments
-            useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
-            useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
-
-            if (useMoments) {
+            if (pmSourcePositionUseMoments(source)) {
                 xPos = source->moments->Mx;
                 yPos = source->moments->My;
@@ -863 +766 @@
             psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
             psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                        fwhmValues->data.F32[entry]);
+            if (fwhmValues) {
+                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                    fwhmValues->data.F32[entry]);
+            } else {
+                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "image is not FWHM-matched",              NAN);
+            }
 
             // XXX if we have raw radial apertures, write them out here
-            psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
-            psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
-            psVector *radFill    = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFlux      = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFluxErr   = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFill      = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32);
             psVectorInit (radFlux,    NAN);
             psVectorInit (radFluxErr, NAN);
@@ -879 +787 @@
             // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
             for (int j = 0; j < radialAper->flux->n; j++) {
-                radFlux->data.F32[j]    = radialAper->flux->data.F32[j];
-                radFluxErr->data.F32[j] = radialAper->fluxErr->data.F32[j];
-                radFill->data.F32[j]    = radialAper->fill->data.F32[j];
+                radFlux->data.F32[j]      = radialAper->flux->data.F32[j];
+                radFluxErr->data.F32[j]   = radialAper->fluxErr->data.F32[j];
+                radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j];
+                radFill->data.F32[j]      = radialAper->fill->data.F32[j];
             }
 
         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_FILL",     PS_DATA_VECTOR, "fill factor of annuli", radFill);
+            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);
             psFree (radFlux);
             psFree (radFluxErr);
+            psFree (radFluxStdev);
             psFree (radFill);
 

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c

@@ -49 +49 @@
 
 #include "pmSourceIO.h"
+#include "pmSourceOutputs.h"
 
 // panstarrs-style FITS table output (header + table in 1st extension)
@@ -62 +63 @@
     psArray *table;
     psMetadata *row;
-    int i;
-    psF32 *PAR, *dPAR;
-    psEllipseAxes axes;
-    psF32 xPos, yPos;
-    psF32 xErr, yErr;
-    psF32 errMag, chisq, apRadius;
-    psS32 nPix, nDOF;
 
     pmChip *chip = readout->parent->parent;
-    pmFPA  *fpa  = chip->parent;
-
-    bool status1 = false;
-    bool status2 = false;
-    float magOffset = NAN;
-    float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
-    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
-    if (!isfinite(zeropt)) {
-        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
-    }
-    if (status1 && status2 && (exptime > 0.0)) {
-        magOffset = zeropt + 2.5*log10(exptime);
-    }
-    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
 
     // if the sequence is defined, write these in seq order; otherwise
@@ -92 +72 @@
         if (source->seq == -1) {
           // let's write these out in S/N order
-          sources = psArraySort (sources, pmSourceSortBySN);
+          sources = psArraySort (sources, pmSourceSortByFlux);
         } else {
           sources = psArraySort (sources, pmSourceSortBySeq);
@@ -100 +80 @@
     table = psArrayAllocEmpty (sources->n);
 
+    short nImageOverlap; 
+    float magOffset; 
+    float zeroptErr; 
+    float fwhmMajor; 
+    float fwhmMinor;
+    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
+
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
     // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
-    for (i = 0; i < sources->n; i++) {
+    for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
 
@@ -113 +100 @@
         }
 
-        // no difference between PSF and non-PSF model
-        pmModel *model = source->modelPSF;
-
-        if (model != NULL) {
-            PAR = model->params->data.F32;
-            dPAR = model->dparams->data.F32;
-            xPos = PAR[PM_PAR_XPOS];
-            yPos = PAR[PM_PAR_YPOS];
-            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
-              xErr = dPAR[PM_PAR_XPOS];
-              yErr = dPAR[PM_PAR_YPOS];
-            } else {
-              // in linear-fit mode, there is no error on the centroid
-              xErr = source->peak->dx;
-              yErr = source->peak->dy;
-            }
-            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0);
-            } else {
-                axes.major = NAN;
-                axes.minor = NAN;
-                axes.theta = NAN;
-            }
-            chisq = model->chisq;
-            nDOF = model->nDOF;
-            nPix = model->nPix;
-            apRadius = source->apRadius;
-            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
-        } else {
-            xPos = source->peak->xf;
-            yPos = source->peak->yf;
-            xErr = source->peak->dx;
-            yErr = source->peak->dy;
-            axes.major = NAN;
-            axes.minor = NAN;
-            axes.theta = NAN;
-            chisq = NAN;
-            nDOF = 0;
-            nPix = 0;
-            apRadius = NAN;
-            errMag = NAN;
-        }
-
-        float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
-        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
-        psS16 nImageOverlap = 1;
-
-        psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
-        float posAngle = 0.0;
-        float pltScale = 0.0;
-        pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
+        // set the 'best' values for various output fields:
+        pmSourceOutputs outputs;
+        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
+
+        pmSourceOutputsMoments moments;
+        pmSourceOutputsSetMoments (&moments, source);
 
         row = psMetadataAlloc ();
         psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr); // XXX this is only measured for non-linear fits
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr); // XXX this is only measured for non-linear fits
-        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr); // XXX this is only measured for non-linear fits
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr); // XXX this is only measured for non-linear fits
+        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                outputs.ra);
+        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               outputs.dec);
+        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
+        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
-        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
+        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
         psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
         psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
-
-        // distinguish moments measure from window vs S/N > XX ??
-        float mxx = source->moments ? source->moments->Mxx : NAN;
-        float mxy = source->moments ? source->moments->Mxy : NAN;
-        float myy = source->moments ? source->moments->Myy : NAN;
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
 
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
@@ -317 +255 @@
         // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
         source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
-        source->peak->flux = peakFlux;
+        source->peak->rawFlux = peakFlux;
+        source->peak->smoothFlux = peakFlux;
         source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
         source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];
-        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
-          source->peak->SN = 1.0 / source->errMag;
-        } else {
-          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
-        }
 
         source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
@@ -339 +273 @@
 
         source->moments = pmMomentsAlloc ();
+        source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
+        source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
+
         source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
         source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
@@ -371 +308 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -549 +486 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

@@ -49 +49 @@
 
 #include "pmSourceIO.h"
+#include "pmSourceOutputs.h"
 
 // panstarrs-style FITS table output (header + table in 1st extension)
@@ -62 +63 @@
     psArray *table;
     psMetadata *row;
-    psF32 *PAR, *dPAR;
-    psEllipseAxes axes;
-    psF32 xPos, yPos;
-    psF32 xErr, yErr;
-    psF32 errMag, chisq, apRadius;
-    psS32 nPix, nDOF;
 
     pmChip *chip = readout->parent->parent;
-    pmFPA  *fpa  = chip->parent;
-
-    bool status1 = false;
-    bool status2 = false;
-    float magOffset = NAN;
-    float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
-    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
-    if (!isfinite(zeropt)) {
-        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
-    }
-    if (status1 && status2 && (exptime > 0.0)) {
-        magOffset = zeropt + 2.5*log10(exptime);
-    }
-    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
 
     // if the sequence is defined, write these in seq order; otherwise
@@ -91 +72 @@
         if (source->seq == -1) {
             // let's write these out in S/N order
-            sources = psArraySort (sources, pmSourceSortBySN);
+            sources = psArraySort (sources, pmSourceSortByFlux);
         } else {
             sources = psArraySort (sources, pmSourceSortBySeq);
@@ -98 +79 @@
 
     table = psArrayAllocEmpty (sources->n);
+
+    short nImageOverlap; 
+    float magOffset; 
+    float zeroptErr; 
+    float fwhmMajor; 
+    float fwhmMinor;
+    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
 
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
@@ -112 +100 @@
         }
 
-        // no difference between PSF and non-PSF model
-        pmModel *model = source->modelPSF;
-
-        if (model != NULL) {
-            PAR = model->params->data.F32;
-            dPAR = model->dparams->data.F32;
-            xPos = PAR[PM_PAR_XPOS];
-            yPos = PAR[PM_PAR_YPOS];
-            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
-                xErr = dPAR[PM_PAR_XPOS];
-                yErr = dPAR[PM_PAR_YPOS];
-            } else {
-                // in linear-fit mode, there is no error on the centroid
-                xErr = source->peak->dx;
-                yErr = source->peak->dy;
-            }
-            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0);
-            } else {
-                axes.major = NAN;
-                axes.minor = NAN;
-                axes.theta = NAN;
-            }
-            chisq = model->chisq;
-            nDOF = model->nDOF;
-            nPix = model->nPix;
-            apRadius = source->apRadius;
-            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
-        } else {
-            xPos = source->peak->xf;
-            yPos = source->peak->yf;
-            xErr = source->peak->dx;
-            yErr = source->peak->dy;
-            axes.major = NAN;
-            axes.minor = NAN;
-            axes.theta = NAN;
-            chisq = NAN;
-            nDOF = 0;
-            nPix = 0;
-            apRadius = NAN;
-            errMag = NAN;
-        }
-
-        float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
-        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
-        psS16 nImageOverlap = 1;
-
-        psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
-        float posAngle = 0.0;
-        float pltScale = 0.0;
-        pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
+        // set the 'best' values for various output fields:
+        pmSourceOutputs outputs;
+        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
+
+        pmSourceOutputsMoments moments;
+        pmSourceOutputsSetMoments (&moments, source);
 
         row = psMetadataAlloc ();
         psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr); // XXX this is only measured for non-linear fits
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr); // XXX this is only measured for non-linear fits
-        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr); // XXX this is only measured for non-linear fits
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr); // XXX this is only measured for non-linear fits
+        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
+        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
-        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
+        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
-        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
+        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
+        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
         psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
         psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeightNotBad);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
-
-        // distinguish moments measure from window vs S/N > XX ??
-        float mxx = source->moments ? source->moments->Mxx : NAN;
-        float mxy = source->moments ? source->moments->Mxy : NAN;
-        float myy = source->moments ? source->moments->Myy : NAN;
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
 
         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
@@ -322 +261 @@
         // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
         source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
-        source->peak->flux = peakFlux;
+        source->peak->rawFlux = peakFlux;
+        source->peak->smoothFlux = peakFlux;
         source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
         source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];
-        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
-          source->peak->SN = 1.0 / source->errMag;
-        } else {
-          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
-        }
 
         source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
@@ -344 +279 @@
 
         source->moments = pmMomentsAlloc ();
+        source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
+        source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
+
         source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
         source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
@@ -392 +330 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -611 +549 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c

@@ -49 +49 @@
 
 #include "pmSourceIO.h"
+#include "pmSourceOutputs.h"
 
 // panstarrs-style FITS table output (header + table in 1st extension)
@@ -62 +63 @@
     psArray *table;
     psMetadata *row;
-    psF32 *PAR, *dPAR;
-    psEllipseAxes axes;
-    psF32 xPos, yPos;
-    psF32 xErr, yErr;
-    psF32 chisq, apRadius;
-    psS32 nPix, nDOF;
 
     pmChip *chip = readout->parent->parent;
-    pmFPA  *fpa  = chip->parent;
-
-    bool status1 = false;
-    bool status2 = false;
-    float magOffset = NAN;
-    float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
-    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
-    if (!isfinite(zeropt)) {
-        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
-    }
-    if (status1 && status2 && (exptime > 0.0)) {
-        magOffset = zeropt + 2.5*log10(exptime);
-    }
-    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
-
-    // we need a measure of the image quality (FWHM) for this image, in order to get the positional errors
-    float fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MAJ");
-    if (!status1) {
-        fwhmMajor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW1");
-        if (!status1) {
-            fwhmMajor = 5.0; // XXX just a guess!
-        }
-    }
-    float fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "FWHM_MIN");
-    if (!status1) {
-        fwhmMinor = psMetadataLookupF32(&status1, readout->analysis, "IQ_FW2");
-        if (!status1) {
-            fwhmMinor = 5.0; // XXX just a guess!
-        }
-    }
 
     // if the sequence is defined, write these in seq order; otherwise
@@ -107 +72 @@
         if (source->seq == -1) {
             // let's write these out in S/N order
-            sources = psArraySort (sources, pmSourceSortBySN);
+            sources = psArraySort (sources, pmSourceSortByFlux);
         } else {
             sources = psArraySort (sources, pmSourceSortBySeq);
@@ -114 +79 @@
 
     table = psArrayAllocEmpty (sources->n);
+
+    short nImageOverlap; 
+    float magOffset; 
+    float zeroptErr; 
+    float fwhmMajor; 
+    float fwhmMinor;
+    pmSourceOutputsCommonValues (&nImageOverlap, &magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
 
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
@@ -128 +100 @@
         }
 
-        // no difference between PSF and non-PSF model
-        pmModel *model = source->modelPSF;
-
-        if (model != NULL) {
-            PAR = model->params->data.F32;
-            dPAR = model->dparams->data.F32;
-            xPos = PAR[PM_PAR_XPOS];
-            yPos = PAR[PM_PAR_YPOS];
-            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
-                xErr = dPAR[PM_PAR_XPOS];
-                yErr = dPAR[PM_PAR_YPOS];
-            } else {
-                xErr = fwhmMajor * source->errMag / 2.35;
-                yErr = fwhmMinor * source->errMag / 2.35;
-            }
-            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0);
-            } else {
-                axes.major = NAN;
-                axes.minor = NAN;
-                axes.theta = NAN;
-            }
-            chisq = model->chisq;
-            nDOF = model->nDOF;
-            nPix = model->nPix;
-            apRadius = source->apRadius;
-        } else {
-            bool useMoments = true;
-            useMoments = (useMoments && source->moments);          // can't if there are no moments
-            useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
-            useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
-
-            if (useMoments) {
-                xPos = source->moments->Mx;
-                yPos = source->moments->My;
-                xErr = fwhmMajor * source->errMag / 2.35;
-                yErr = fwhmMinor * source->errMag / 2.35;
-            } else {
-                xPos = source->peak->xf;
-                yPos = source->peak->yf;
-                xErr = source->peak->dx;
-                yErr = source->peak->dy;
-            }
-            axes.major = NAN;
-            axes.minor = NAN;
-            axes.theta = NAN;
-            chisq = NAN;
-            nDOF = 0;
-            nPix = 0;
-            apRadius = NAN;
-        }
-
-        float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
-        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
-        psS16 nImageOverlap = 1;
-
-        psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
-        float posAngle = 0.0;
-        float pltScale = 0.0;
-        pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
+        // set the 'best' values for various output fields:
+        pmSourceOutputs outputs;
+        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
+
+        pmSourceOutputsMoments moments;
+        pmSourceOutputsSetMoments (&moments, source);
 
         row = psMetadataAlloc ();
         psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
-        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr);
-        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr);
-        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
+        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
+        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);
@@ -203 +121 @@
         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
         psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
-        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
-
-        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
         
         // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
-        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
-        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
-
-        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
+        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
         psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
         psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
 
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
-
-        // distinguish moments measure from window vs S/N > XX ??
-        float Mxx = source->moments ? source->moments->Mxx : NAN;
-        float Mxy = source->moments ? source->moments->Mxy : NAN;
-        float Myy = source->moments ? source->moments->Myy : NAN;
-
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      Mxx);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      Mxy);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      Myy);
-
-        float M_c3 = source->moments ? 1.0*source->moments->Mxxx - 3.0*source->moments->Mxyy : NAN;
-        float M_s3 = source->moments ? 3.0*source->moments->Mxxy - 1.0*source->moments->Myyy : NAN;
-        float M_c4 = source->moments ? 1.0*source->moments->Mxxxx - 6.0*source->moments->Mxxyy + 1.0*source->moments->Myyyy : NAN;
-        float M_s4 = source->moments ? 4.0*source->moments->Mxxxy - 4.0*source->moments->Mxyyy : NAN;
-
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                    M_c3);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                    M_s3);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                   M_c4);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                   M_s4);
-
-        float Mrf  = source->moments ? source->moments->Mrf : NAN;
-        float Mrh  = source->moments ? source->moments->Mrh : NAN;
-        float Krf  = source->moments ? source->moments->KronFlux : NAN;
-        float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
-
-        float Kinner = source->moments ? source->moments->KronFinner : NAN;
-        float Kouter = source->moments ? source->moments->KronFouter : NAN;
-
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                       Mrf);
-        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                        Mrh);
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Krf);
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                          dKrf);
-
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kinner);
-        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                     Kouter);
-
-        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                     source->mode);
-        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
+        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
+        // 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);
+        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
+        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
         psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
 
@@ -384 +284 @@
         // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
         source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
-        source->peak->flux = peakFlux;
+        source->peak->rawFlux = peakFlux;
+        source->peak->smoothFlux = peakFlux;
         source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
         source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];
-        if (isfinite (source->errMag) && (source->errMag > 0.0)) {
-          source->peak->SN = 1.0 / source->errMag;
-        } else {
-          source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
-        }
+
+        // we no longer sort by S/N, only flux
+        // if (isfinite (source->errMag) && (source->errMag > 0.0)) {
+        //   source->peak->SN = 1.0 / source->errMag;
+        // } else {
+        //   source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N
+        // }
 
         source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
@@ -407 +310 @@
 
         source->moments = pmMomentsAlloc ();
+        source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
+        source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
+
         source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
         source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
@@ -477 +383 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -649 +555 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams

trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c

@@ -91 +91 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -136 +136 @@
 
         float calMag = source->psfMag + magOffset;
-        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
+        float peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
         psS16 nImageOverlap = 1;
 
@@ -294 +294 @@
 
         source->peak       = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
-        source->peak->flux = peakFlux;
+        source->peak->rawFlux = peakFlux;
+        source->peak->smoothFlux = peakFlux;
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];
@@ -357 +358 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -576 +577 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams

trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c

@@ -101 +101 @@
         }
 
-        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
+        float peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
         psS16 nImageOverlap = 1;
         psS32 ID = 0; // XXX need to figure out how to generate this
@@ -220 +220 @@
 
         source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
-        source->peak->flux = peakFlux;
+        source->peak->rawFlux = peakFlux;
+        source->peak->smoothFlux = peakFlux;
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];

trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

@@ -73 +73 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -117 +117 @@
         }
 
-        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
+        float peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
         psS16 nImageOverlap = 1;
 
@@ -263 +263 @@
 
         source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
-        source->peak->flux = peakFlux;
+        source->peak->rawFlux = peakFlux;
+        source->peak->smoothFlux = peakFlux;
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];
@@ -307 +308 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     table = psArrayAllocEmpty (sources->n);
@@ -480 +481 @@
 
     // let's write these out in S/N order
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams

trunk/psModules/src/objects/pmSourceIO_RAW.c

@@ -126 +126 @@
                  source[0].apMag, source[0].type, source[0].mode,
                  logChi, logChiNorm,
-                 source[0].peak->SN,
+                 source[0].peak->rawFlux,
                  source[0].apRadius,
                  source[0].pixWeightNotBad,
@@ -188 +188 @@
                  log10(model[0].chisq/model[0].nDOF),
                  log10(model[0].chisqNorm/model[0].nDOF),
-                 source[0].peak->SN,
+                 source[0].peak->rawFlux,
                  source[0].apRadius,
                  source[0].pixWeightNotBad,
@@ -234 +234 @@
 
         fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n",
-                 source->peak->x, source->peak->y, source->peak->value,
+                 source->peak->x, source->peak->y, source->peak->detValue,
                  moment->Mx, moment->My,
                  moment->Mxx, moment->Myy,
@@ -270 +270 @@
             continue;
         fprintf (f, "%5d %5d  %8.1f  %7.1f %7.1f  %6.3f %6.3f  %10.1f %7.1f %7.1f %7.1f  %4d %2d %#5x\n",
-                 source->peak->x, source->peak->y, source->peak->value,
+                 source->peak->x, source->peak->y, source->peak->detValue,
                  source->moments->Mx, source->moments->My,
                  source->moments->Mxx, source->moments->Myy,
@@ -299 +299 @@
             continue;
         fprintf (f, "%5d %5d  %7.1f %7.2f %7.2f %7.2f\n",
-                 peak->x, peak->y, peak->value, peak->SN, peak->xf, peak->yf);
-    }
-    fclose (f);
-    return true;
-}
-
+                 peak->x, peak->y, peak->detValue, peak->rawFlux, peak->xf, peak->yf);
+    }
+    fclose (f);
+    return true;
+}
+

trunk/psModules/src/objects/pmSourceIO_SMPDATA.c

@@ -205 +205 @@
 
         source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
-        source->peak->flux = peakFlux;
 
         sources->data[i] = source;

trunk/psModules/src/objects/pmSourceMasks.h

@@ -46 +46 @@
     PM_SOURCE_MODE2_DIFF_WITH_DOUBLE = 0x00000002, ///< diff source matched to positive detections in both images
     PM_SOURCE_MODE2_MATCHED          = 0x00000004, ///< diff source matched to positive detections in both images
+
+    PM_SOURCE_MODE2_ON_SPIKE         = 0x00000008, ///< > 25% of (PSF-weighted) pixels land on diffraction spike
+    PM_SOURCE_MODE2_ON_STARCORE      = 0x00000010, ///< > 25% of (PSF-weighted) pixels land on starcore
+    PM_SOURCE_MODE2_ON_BURNTOOL      = 0x00000020, ///< > 25% of (PSF-weighted) pixels land on burntool
+    PM_SOURCE_MODE2_ON_CONVPOOR      = 0x00000040, ///< > 25% of (PSF-weighted) pixels land on convpoor
 } pmSourceMode2;
 

trunk/psModules/src/objects/pmSourceMoments.c

@@ -82 +82 @@
     psF32 Sum = 0.0;
     psF32 Var = 0.0;
+    psF32 SumCore = 0.0;
+    psF32 VarCore = 0.0;
     psF32 R2 = PS_SQR(radius);
     psF32 minSN2 = PS_SQR(minSN);
@@ -248 +250 @@
 
     int nKronPix = 0;
+    int nCorePix = 0;
+    int nInner = 0;
+    int nOuter = 0;
     Sum = Var = 0.0;
     float SumInner = 0.0;
@@ -255 +260 @@
 
         psF32 yDiff = row - yCM;
-        if (fabs(yDiff) > radKron) continue;
+        if (fabs(yDiff) > radKouter) continue;
 
         psF32 *vPix = source->pixels->data.F32[row];
@@ -272 +277 @@
 
             psF32 xDiff = col - xCM;
-            if (fabs(xDiff) > radKron) continue;
+            if (fabs(xDiff) > radKouter) continue;
 
             // radKron is just a function of (xDiff, yDiff)
@@ -293 +298 @@
             }
 
+            // 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 ++;
             }
         }
     }
-    source->moments->KronFlux = Sum;
-    source->moments->KronFinner = SumInner;
-    source->moments->KronFouter = SumOuter;
-    source->moments->KronFluxErr = sqrt(Var);
+    // *** 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",
@@ -313 +330 @@
 
     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->flux, source->peak->SN, source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, source->moments->nPixels);
+             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);

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -23 +23 @@
 #include "pmFPAMaskWeight.h"
 
+#include "pmConfigMask.h"
 #include "pmTrend2D.h"
 #include "pmResiduals.h"
@@ -49 +50 @@
 static float AP_MIN_SN = 0.0;
 
-bool pmSourceMagnitudesInit (psMetadata *config)
-{
-    PS_ASSERT_PTR_NON_NULL(config, false);
+// make this a bit more clever and dynamic
+static psImageMaskType maskSuspect  = 0;
+static psImageMaskType maskSpike    = 0;
+static psImageMaskType maskStarCore = 0;
+static psImageMaskType maskBurntool = 0;
+static psImageMaskType maskConvPoor = 0;
+
+bool pmSourceMagnitudesInit (pmConfig *config, psMetadata *recipe)
+{
+    PS_ASSERT_PTR_NON_NULL(recipe, false);
     bool status;
 
-    float limit = psMetadataLookupF32 (&status, config, "AP_MIN_SN");
+    // we are going to test specially against these poor values
+    if (config) {
+        maskSpike    = pmConfigMaskGet("SPIKE", config);
+        maskStarCore = pmConfigMaskGet("STARCORE", config);
+        maskBurntool = pmConfigMaskGet("BURNTOOL", config);
+        maskConvPoor = pmConfigMaskGet("CONV.POOR", config);
+        maskSuspect  = maskSpike | maskStarCore | maskBurntool | maskConvPoor;
+    }
+
+    float limit = psMetadataLookupF32 (&status, recipe, "AP_MIN_SN");
     if (status) {
         AP_MIN_SN = limit;
@@ -77 +94 @@
 // XXX masked region should be (optionally) elliptical
 // if mode is PM_SOURCE_PHOT_PSFONLY, we skip all other magnitudes
-bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal)
+bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal, float radius)
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -166 +183 @@
     // measure the contribution of included pixels
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (&source->pixWeightNotBad, &source->pixWeightNotPoor, model, source->maskObj, maskVal, markVal);
+        pmSourcePixelWeight (source, model, source->maskObj, maskVal, radius);
     }
 
@@ -342 +359 @@
 
 // return source aperture magnitude
-bool pmSourcePixelWeight (float *pixWeightNotBad, float *pixWeightNotPoor, pmModel *model, psImage *mask, psImageMaskType maskVal, psImageMaskType markVal)
-{
-    PS_ASSERT_PTR_NON_NULL(pixWeightNotBad, false);
-    PS_ASSERT_PTR_NON_NULL(pixWeightNotPoor, false);
+bool pmSourcePixelWeight (pmSource *source, pmModel *model, psImage *mask, psImageMaskType maskVal, float radius)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    source->pixWeightNotBad = NAN;
+    source->pixWeightNotPoor = NAN;
+
     PS_ASSERT_PTR_NON_NULL(mask, false);
     PS_ASSERT_PTR_NON_NULL(model, false);
@@ -355 +374 @@
     float value;
 
+    float spikeSum = 0;
+    float starcoreSum = 0;
+    float burntoolSum = 0;
+    float convpoorSum = 0;
+
     int Xo, Yo, dP;
     int dX, DX, NX;
     int dY, DY, NY;
 
-    *pixWeightNotBad = 0.0;
-    *pixWeightNotPoor = 0.0;
+    float radius2 = PS_SQR(radius);
 
     // we only care about the value of the object model, not the local sky
@@ -387 +410 @@
     NY = mask->numRows;
 
+    psImageMaskType maskBad = maskVal;
+    maskBad &= ~maskSuspect;
+
     // measure modelSum and validSum.  this function is applied to a sources' subimage.  the
     // value of DX is chosen (see above) to cover the full possible size of the subimage if it
     // were not by an edge; ie, if the source is cut in half by an image edge, we correctly
     // count the virtual pixels off the edge in normalizing the value of the pixWeight
+
+    // we skip any pixels [real or virtual] outside of the specified radius (nominally the aperture radius)
     for (int ix = -DX; ix < DX + 1; ix++) {
+        if (ix > radius) continue;
         int mx = ix + dX;
         for (int iy = -DY; iy < DY + 1; iy++) {
+            if (iy > radius) continue;
+            if (ix*ix + iy*iy > radius2) continue;
             int my = iy + dY;
 
@@ -409 +440 @@
             if (my >= NY) continue;
 
-            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)) {
+            // count pixels which are masked only with bad pixels
+            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskBad)) {
                 notBadSum += value;
             }
-            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & ~markVal)) {
+
+            // count pixels which are masked with an mask bit (bad or poor)
+            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)) {
                 notPoorSum += value;
             }
+
+            // count pixels which are masked with an mask bit (bad or poor)
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskSpike) {
+                spikeSum += value;
+            }
+            // count pixels which are masked with an mask bit (bad or poor)
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskStarCore) {
+                starcoreSum += value;
+            }
+            // count pixels which are masked with an mask bit (bad or poor)
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskBurntool) {
+                burntoolSum += value;
+            }
+            // count pixels which are masked with an mask bit (bad or poor)
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskConvPoor) {
+                convpoorSum += value;
+            }
         }
     }
     psFree (coord);
 
-    *pixWeightNotBad  = notBadSum  / modelSum;
-    *pixWeightNotPoor = notPoorSum / modelSum;
+    source->pixWeightNotBad  = notBadSum  / modelSum;
+    source->pixWeightNotPoor = notPoorSum / modelSum;
+
+    if ((spikeSum/modelSum) > 0.25) {
+        source->mode2 |= PM_SOURCE_MODE2_ON_SPIKE;
+    }
+    if ((starcoreSum/modelSum) > 0.25) {
+        source->mode2 |= PM_SOURCE_MODE2_ON_STARCORE;
+    }
+    if ((burntoolSum/modelSum) > 0.25) {
+        source->mode2 |= PM_SOURCE_MODE2_ON_BURNTOOL;
+    }
+    if ((convpoorSum/modelSum) > 0.25) {
+        source->mode2 |= PM_SOURCE_MODE2_ON_CONVPOOR;
+    }
+
+    if (isfinite(source->pixWeightNotBad) && isfinite(source->pixWeightNotPoor)) {
+        psAssert (source->pixWeightNotBad <= source->pixWeightNotPoor, "error: all bad pixels should also be poor");
+    }
 
     return (true);
@@ -427 +495 @@
 # define FLUX_LIMIT 3.0
 
-// measure stats that may be using in difference images for distinguishing real sources from bad residuals
+// measure stats that may be used in difference images for distinguishing real sources from bad residuals
 bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
 {
@@ -673 +741 @@
 # endif
 
-// determine chisq, etc for linear normalization-only fit
-bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor)
+// determine chisq, nPix, nDOF, chisqNorm : model->nPar must be set
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(model, false);
@@ -689 +757 @@
             if (variance->data.F32[j][i] <= 0)
                 continue;
-            dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]);
+            dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i];
             Npix ++;
         }
     }
     model->nPix = Npix;
-    model->nDOF = Npix - 1;
+    model->nDOF = Npix - model->nPar;
     model->chisq = dC;
+    model->chisqNorm = dC / model->nDOF;
 
     return (true);
 }
 
+
+// return source aperture magnitude
+bool pmSourceChisqUnsubtracted (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;
+
+    // the model function returns the source flux at a position
+    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+    psVector *params = model->params;
+    psImage  *image = source->pixels;
+    psImage  *mask = source->maskObj;
+    psImage  *variance = source->variance;
+
+    int dX = image->col0;
+    int dY = image->row0;
+
+    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;
+
+            coord->data.F32[0] = (psF32) ix + dX + 0.5;
+            coord->data.F32[1] = (psF32) iy + dY + 0.5;
+
+            // for the full model, add all points
+            float value = model->modelFunc (NULL, params, coord);
+
+            // fprintf (stderr, "%d, %d : %f, %f : %f - %f : %f\n", 
+            // ix, iy, coord->data.F32[0], coord->data.F32[1], image->data.F32[iy][ix], value, dC);
+
+            dC += PS_SQR (image->data.F32[iy][ix] - value) / variance->data.F32[iy][ix];
+            Npix ++;
+        }
+    }
+    model->nPix = Npix;
+    model->nDOF = Npix - model->nPar;
+    model->chisq = dC;
+    model->chisqNorm = dC / model->nDOF;
+
+    psFree (coord);
+    return (true);
+}
 
 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)

trunk/psModules/src/objects/pmSourcePhotometry.h

@@ -64 +64 @@
 );
 
-bool pmSourceMagnitudesInit (psMetadata *config);
-bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal);
+bool pmSourceMagnitudesInit (pmConfig *config, psMetadata *recipe);
+bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode, psImageMaskType maskVal, psImageMaskType markVal, float radius);
 
-bool pmSourcePixelWeight (float *pixWeightNotBad, float *pixWeightNotPoor, pmModel *model, psImage *mask, psImageMaskType maskVal, psImageMaskType markVal);
+bool pmSourcePixelWeight (pmSource *source, pmModel *model, psImage *mask, psImageMaskType maskVal, float radius);
 
-bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor);
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal);
+bool pmSourceChisqUnsubtracted (pmSource *source, pmModel *model, psImageMaskType maskVal);
 
 bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);

trunk/psModules/src/objects/pmSourceVisual.c

@@ -56 +56 @@
 
     if (!pmVisualTestLevel("psphot.psf.metric", 2)) return true;
-
-    if (kapa1 == -1) {
-        kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots");
-        if (kapa1 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    if (!pmVisualInitWindow (&kapa1, "pmSource:plots")) return false;
 
     KapaClearSections (kapa1);
@@ -140 +132 @@
 
     if (!pmVisualTestLevel("psphot.psf.subpix", 3)) return true;
-
-    if (kapa1 == -1) {
-        kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots");
-        if (kapa1 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    if (!pmVisualInitWindow (&kapa1, "pmSource:plots")) return false;
 
     KapaClearSections (kapa1);
     KapaInitGraph (&graphdata);
+    section.bg = KapaColorByName ("none"); // XXX probably should be 'none'
 
     int n;
@@ -302 +287 @@
 
     if (!pmVisualTestLevel("psphot.psf.fits", 2)) return true;
-
-    if (kapa2 == -1) {
-        kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images");
-        if (kapa2 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    if (!pmVisualInitWindow (&kapa2, "pmSource:images")) return false;
 
     // create images 1/10 scale:
@@ -385 +362 @@
     if (!source->pixels) return false;
     if (!source->modelFlux) return false;
-
-    if (kapa2 == -1) {
-        kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images");
-        if (kapa2 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    if (!pmVisualInitWindow (&kapa2, "pmSource:images")) return false;
 
     // KapaClearSections (kapa2);
@@ -428 +397 @@
     if (!plotPSF) return true;
     if (!pmVisualTestLevel("psphot.psf.resid", 2)) return true;
-
-    if (kapa1 == -1) {
-        kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots");
-        if (kapa1 == -1) {
-            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
-            pmVisualSetVisual(false);
-            return false;
-        }
-    }
+    if (!pmVisualInitWindow (&kapa1, "pmSource:plots")) return false;
 
     KapaClearPlots (kapa1);
     KapaInitGraph (&graphdata);
+    section.bg = KapaColorByName ("none"); // XXX probably should be 'none'
 
     float Xmin = +1e32;

@@ -28 +28 @@
 ChangeLog
 psmodules-*.tar.*
+a.out.dSYM