Changeset 31451

Timestamp:

May 5, 2011, 11:02:53 AM (15 years ago)

Author:

eugene

Message:

merge updates from eam branch 20110404

Location:

trunk/psModules

Files:

• src/camera/pmFPAMaskWeight.c (modified) (1 diff)
• src/imcombine/pmPSFEnvelope.c (modified) (1 diff)
• src/imcombine/pmSubtractionStamps.c (modified) (1 diff)
• src/objects/Makefile.am (modified) (1 diff)
• src/objects/models/pmModel_DEV.c (modified) (1 diff)
• src/objects/models/pmModel_EXP.c (modified) (1 diff)
• src/objects/models/pmModel_PGAUSS.c (modified) (1 diff)
• src/objects/models/pmModel_PS1_V1.c (modified) (1 diff)
• src/objects/models/pmModel_QGAUSS.c (modified) (1 diff)
• src/objects/models/pmModel_RGAUSS.c (modified) (1 diff)
• src/objects/models/pmModel_SERSIC.c (modified) (1 diff)
• src/objects/pmGrowthCurve.c (modified) (3 diffs)
• src/objects/pmGrowthCurve.h (modified) (1 diff)
• src/objects/pmGrowthCurveGenerate.c (modified) (3 diffs)
• src/objects/pmGrowthCurveGenerate.h (copied) (copied from branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.h )
• src/objects/pmPCMdata.c (modified) (2 diffs)
• src/objects/pmPSF.c (modified) (3 diffs)
• src/objects/pmPSF.h (modified) (5 diffs)
• src/objects/pmPSF_IO.c (modified) (11 diffs)
• src/objects/pmPSFtryFitPSF.c (modified) (2 diffs)
• src/objects/pmPeaks.c (modified) (3 diffs)
• src/objects/pmPeaks.h (modified) (2 diffs)
• src/objects/pmSource.c (modified) (9 diffs)
• src/objects/pmSource.h (modified) (4 diffs)
• src/objects/pmSourceIO_CMF_PS1_DV1.c (modified) (2 diffs)
• src/objects/pmSourceIO_CMF_PS1_DV2.c (modified) (2 diffs)
• src/objects/pmSourceIO_CMF_PS1_SV1.c (modified) (3 diffs)
• src/objects/pmSourceIO_CMF_PS1_V1.c (modified) (2 diffs)
• src/objects/pmSourceIO_CMF_PS1_V2.c (modified) (2 diffs)
• src/objects/pmSourceIO_CMF_PS1_V3.c (modified) (4 diffs)
• src/objects/pmSourceIO_CMP.c (modified) (3 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) (2 diffs)
• src/objects/pmSourceIO_SMPDATA.c (modified) (2 diffs)
• src/objects/pmSourceIO_SX.c (modified) (1 diff)
• src/objects/pmSourceMasks.h (modified) (2 diffs)
• src/objects/pmSourceMatch.c (modified) (2 diffs)
• src/objects/pmSourceMoments.c (modified) (19 diffs)
• src/objects/pmSourceOutputs.c (modified) (2 diffs)
• src/objects/pmSourcePhotometry.c (modified) (14 diffs)
• src/objects/pmSourcePhotometry.h (modified) (1 diff)
• src/objects/pmSourceUtils.h (modified) (1 diff)
• src/psmodules.h (modified) (1 diff)
• test/objects/tap_pmSource.c (modified) (1 diff)
• test/objects/tap_pmSourceIO_PS1_DEV_0.c (modified) (2 diffs)
• test/objects/tap_pmSourceIO_PS1_DEV_1.c (modified) (3 diffs)
• test/objects/tap_pmSourceIO_SMPDATA.c (modified) (2 diffs)
• test/objects/tap_pmSourcePhotometry.c (modified) (1 diff)

trunk/psModules/src/camera/pmFPAMaskWeight.c

@@ -587 +587 @@
                 double imageValue, varianceValue; // Image and variance value from interpolation
                 psImageMaskType maskValue = 0; // Mask value from interpolation
-                psImageInterpolateStatus status = psImageInterpolate(&imageValue, &varianceValue, &maskValue, x, y, interp);
+
+                // interpolate to pixel center (index + 0.5)
+                psImageInterpolateStatus status = psImageInterpolate(&imageValue, &varianceValue, &maskValue, x + 0.5, y + 0.5, interp);
                 if (status == PS_INTERPOLATE_STATUS_ERROR || status == PS_INTERPOLATE_STATUS_OFF) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to interpolate readout at %d,%d", x, y);

trunk/psModules/src/imcombine/pmPSFEnvelope.c

@@ -380 +380 @@
 
         // measure the source moments: tophat windowing, no pixel S/N cutoff
-        if (!pmSourceMoments(source, maxRadius, 0.0, 0.0, maskVal)) {
+        if (!pmSourceMoments(source, maxRadius, 0.0, 0.0, 0.0, maskVal)) {
             // Can't do anything about it; limp along as best we can
             psErrorClear();

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -1197 +1197 @@
        
         // fprintf (stderr, "%f,%f : %f %f : %f %f\n", source->peak->xf, source->peak->yf,
-        // source->psfMag, source->apMag, source->psfMag - source->apMag, source->errMag);
+        // source->psfMag, source->apMag, source->psfMag - source->apMag, source->psfMagErr);
 
         // XXX this is somewhat arbitrary...
-        if (source->errMag > 0.05) continue;
+        if (source->psfMagErr > 0.05) continue;
         if (fabs(source->psfMag - source->apMag) > 0.5) continue;
 

trunk/psModules/src/objects/Makefile.am

@@ -114 +114 @@
         pmTrend2D.h \
         pmGrowthCurve.h \
+        pmGrowthCurveGenerate.h \
         pmSourceMatch.h \
         pmDetEff.h \

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

@@ -268 +268 @@
     float f0 = 1.0;
     float f1, f2;
-    for (z = DZ; z < 50; z += DZ) {
+    for (z = DZ; z < 150; z += DZ) {
         f1 = exp(-pow(z,ALPHA));
         z += DZ;

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

@@ -255 +255 @@
     float f0 = 1.0;
     float f1, f2;
-    for (z = DZ; z < 50; z += DZ) {
+    for (z = DZ; z < 150; z += DZ) {
         f1 = exp(-sqrt(z));
         z += DZ;

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

@@ -239 +239 @@
     float f0 = 1.0;
     float f1, f2;
-    for (z = DZ; z < 50; z += DZ) {
+    for (z = DZ; z < 150; z += DZ) {
         f1 = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
         z += DZ;

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

@@ -261 +261 @@
     float f0 = 1.0;
     float f1, f2;
-    for (z = DZ; z < 50; z += DZ) {
+    for (z = DZ; z < 150; z += DZ) {
         f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         z += DZ;

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

@@ -262 +262 @@
     float f0 = 1.0;
     float f1, f2;
-    for (z = DZ; z < 50; z += DZ) {
+    for (z = DZ; z < 150; z += DZ) {
         f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         z += DZ;

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

@@ -251 +251 @@
     float f0 = 1.0;
     float f1, f2;
-    for (z = DZ; z < 50; z += DZ) {
+    for (z = DZ; z < 150; z += DZ) {
         f1 = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
         z += DZ;

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

@@ -316 +316 @@
     float f0 = 1.0;
     float f1, f2;
-    for (z = DZ; z < 50; z += DZ) {
+    for (z = DZ; z < 150; z += DZ) {
         // f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
         f1 = exp(-pow(z,PAR[PM_PAR_7]));

trunk/psModules/src/objects/pmGrowthCurve.c

@@ -73 +73 @@
     // Fractional pixel radii are not well defined; use integer pixel radii.  Use 1 pixel steps
     // until the scaling factor steps in intervals larger than 1 pixel
-    float Rlin = 1.0 / (fR - 1.0);
+    // float Rlin = 1.0 / (fR - 1.0);
 
     growth->radius = psVectorAllocEmpty (NPTS, PS_DATA_F32);
@@ -79 +79 @@
     // there will be NPTS radii + a few extras 
     float radius = minRadius;
-    while (radius < Rlin) {
+    while (radius < refRadius) {
         // fprintf (stderr, "r: %f\n", radius);
-        psVectorAppend (growth->radius, radius);
+        psVectorAppend (growth->radius, radius - 0.001);
         radius += 1.0;
     }    
+    growth->refBin = growth->radius->n - 1;
     while (radius < maxRadius) {
         // fprintf (stderr, "r: %f\n", radius);
-        psVectorAppend (growth->radius, radius);
+        psVectorAppend (growth->radius, radius - 0.001);
         radius *= fR;
         radius = (int) (radius + 0.5);
@@ -96 +97 @@
     growth->refRadius = refRadius;
     growth->maxRadius = maxRadius;
-    growth->apLoss = 0.0;
-    growth->fitMag = 0.0;
+    growth->fitMag = NAN;
+    growth->refMag = NAN;
+    growth->apRef  = NAN;
+    growth->apLoss = NAN;
+
     return growth;
 }

trunk/psModules/src/objects/pmGrowthCurve.h

@@ -22 +22 @@
     psF32 maxRadius;
     psF32 fitMag;
+    psF32 refMag;
     psF32 apRef;   // apMag[refRadius]
     psF32 apLoss;  // fitMag - apRef
+    int refBin;
 }
 pmGrowthCurve;

trunk/psModules/src/objects/pmGrowthCurveGenerate.c

@@ -50 +50 @@
 
 #include "pmSourcePhotometry.h"
-
-pmGrowthCurve *pmGrowthCurveForPosition (psImage *image, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType markVal, float xc, float yc);
+#include "pmGrowthCurveGenerate.h"
 
 /*****************************************************************************/
@@ -197 +196 @@
         // mask the given aperture and measure the apMag
         psImageKeepCircle (mask, xc, yc, radius, "OR", markVal);
-        if (!pmSourcePhotometryAper (&apMag, NULL, NULL, model, pixels, NULL, mask, maskVal)) {
+        if (!pmSourcePhotometryAper (NULL, &apMag, NULL, NULL, model, pixels, NULL, mask, maskVal)) {
             psFree (growth);
             psFree (view);
@@ -226 +225 @@
     return growth;
 }
+
+# define DEBUG 0
+# if (DEBUG)
+static FILE *fgr = NULL;
+# endif
+
+// we generate the growth curve for the center of the image with the specified psf model
+bool pmGrowthCurveGenerateFromSources (pmReadout *readout, pmPSF *psf, psArray *sources, bool INTERPOLATE_AP, psImageMaskType maskVal, psImageMaskType markVal)
+{
+    PS_ASSERT_PTR_NON_NULL(readout, false);
+    PS_ASSERT_PTR_NON_NULL(readout->image, false);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    pmSourcePhotometryMode photMode = INTERPOLATE_AP ? PM_SOURCE_PHOT_INTERP : 0;
+    
+    // measure the growth curve for each PSF source and average them together
+    psArray *growths = psArrayAllocEmpty (100);
+
+# if (DEBUG)
+    fgr = fopen ("growth.mags.dat", "w");
+# endif
+
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+
+        pmGrowthCurve *growth = pmGrowthCurveForSource (source, psf, photMode, maskVal, markVal);
+        if (!growth) continue;
+        
+        psArrayAdd (growths, 100, growth);
+        psFree (growth);
+    }
+    psAssert (growths->n, "cannot build growth curve (no valid PSF stars?)");
+
+# if (DEBUG)
+    fclose (fgr);
+# endif
+
+    // just use a simple sample median to get the 'best' value from each growth curve...
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    psVector *values = psVectorAlloc (growths->n, PS_DATA_F32);
+
+    // loop over a range of source fluxes
+    // no need to interpolate since we have forced the object center
+    // to 0.5, 0.5 above
+    for (int i = 0; i < psf->growth->radius->n; i++) {
+
+        // median the values for each radial bin
+        values->n = 0;
+        for (int j = 0; j < growths->n; j++) {
+            pmGrowthCurve *growth = growths->data[j];
+            if (!isfinite(growth->apMag->data.F32[i])) continue;
+            psVectorAppend (values, growth->apMag->data.F32[i] - growth->refMag);
+        }
+        if (values->n == 0) {
+            psf->growth->apMag->data.F32[i] = NAN;
+        } else {
+            if (!psVectorStats (stats, values, NULL, NULL, 0)) {
+                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+                return false;
+            }
+            psf->growth->apMag->data.F32[i] = stats->sampleMedian;
+        }
+    }
+
+    psf->growth->fitMag = psf->growth->apMag->data.F32[psf->growth->radius->n-1];
+    psf->growth->apRef = psVectorInterpolate (psf->growth->radius, psf->growth->apMag, psf->growth->refRadius);
+    psf->growth->apLoss = psf->growth->fitMag - psf->growth->apRef;
+
+    psLogMsg ("psphot.growth", 4, "GrowthCurve : apLoss : %f (fitMag - apMag @ ref : %f - %f)\n", psf->growth->apLoss, psf->growth->fitMag, psf->growth->apRef);
+
+    psFree (growths);
+    psFree (stats);
+    psFree (values);
+
+    return true;
+}
+
+pmGrowthCurve *pmGrowthCurveForSource (pmSource *source, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal, psImageMaskType markVal) {
+
+    float radius;
+
+    assert (psf->growth);
+
+    float minRadius = psf->growth->radius->data.F32[0];
+    pmGrowthCurve *growth = pmGrowthCurveAlloc (minRadius, psf->growth->maxRadius, psf->growth->refRadius);
+
+    // measure the fitMag for this source (for normalization)
+    // pmSourcePhotometryModel (&fitMag, NULL, source->psfModel);
+    growth->fitMag = source->psfMag;
+
+    float xc = source->peak->xf;
+    float yc = source->peak->yf;
+
+    // Loop over the range of radii
+    for (int i = 0; i < growth->radius->n; i++) {
+
+        radius = growth->radius->data.F32[i];
+
+        // mask the given aperture and measure the apMag
+        psImageKeepCircle (source->maskObj, xc, yc, radius, "OR", markVal);
+
+        if (!pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, radius)) {
+            psFree (growth);
+            return NULL;
+        }
+
+        // if (!pmSourcePhotometryAper (NULL, &apMag, NULL, NULL, NULL, source->pixels, NULL, source->maskObj, maskVal)) {
+        //     psFree (growth);
+        //     return NULL;
+        // }
+
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+
+        growth->apMag->data.F32[i] = source->apMag;
+    }
+    psAssert(growth->refBin >= 0, "invalid growth reference bin");
+    psAssert(growth->refBin < growth->apMag->n, "invalid growth reference bin");
+    growth->refMag = growth->apMag->data.F32[growth->refBin];
+
+    // Loop over the range of radii
+# if (DEBUG)
+    for (int i = 0; i < growth->radius->n; i++) {
+        fprintf (fgr, "%f %f  %f %f %f %f\n", xc, yc, growth->radius->data.F32[i], growth->apMag->data.F32[i], growth->fitMag, growth->refMag);
+    }
+# endif
+
+    return growth;
+}

trunk/psModules/src/objects/pmPCMdata.c

@@ -263 +263 @@
 bool pmPCMupdate(pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model) {
 
-    bool newWindow = (source->pixels->numRows != pcm->modelFlux->numRows) || (source->pixels->numCols != pcm->modelFlux->numCols);
+    bool sameWindow = (source->pixels->numRows == pcm->modelFlux->numRows);
+    sameWindow     &= (source->pixels->numCols == pcm->modelFlux->numCols);
+    sameWindow     &= (source->pixels->col0    == pcm->modelFlux->col0);
+    sameWindow     &= (source->pixels->row0    == pcm->modelFlux->row0);
 
     // re-count the number of unmasked pixels:
-    if (newWindow) {
+    if (!sameWindow) {
         for (psS32 i = 0; i < source->pixels->numRows; i++) {
             for (psS32 j = 0; j < source->pixels->numCols; j++) {
@@ -346 +349 @@
 
     // has the source pixel window changed?
-    if (newWindow) {
+    if (!sameWindow) {
 
         // adjust all supporting images:

trunk/psModules/src/objects/pmPSF.c

@@ -73 +73 @@
 
     options->type          = 0;
-
     options->stats         = NULL;
-    options->fitOptions    = NULL; // XXX this has to be set before calling pmPSF fit functions
 
     options->psfTrendMode  = PM_TREND_NONE;
@@ -90 +88 @@
 
     options->chiFluxTrend = true;
-
+    options->fitOptions    = NULL; // XXX this has to be set before calling pmPSF fit functions
+
+    options->fitRadius = NAN;
+    options->apRadius = NAN;
     return options;
 }
@@ -141 +142 @@
 
     psf->type     = options->type;
+
     psf->chisq    = 0.0;
     psf->ApResid  = 0.0;

trunk/psModules/src/objects/pmPSF.h

@@ -31 +31 @@
 struct pmPSF {
     pmModelType type;                   ///< PSF Model in use
-    psArray *params;                    ///< Model parameters (psPolynomial2D)
-    psStats *psfTrendStats;             ///< psf parameter trend clipping stats
-    pmTrend2DMode psfTrendMode;
-    psPolynomial1D *ChiTrend;           ///< Chisq vs flux fit (correction for systematic errors)
-    pmTrend2D *ApTrend;                 ///< ApResid vs (x,y)
-    pmTrend2D *FluxScale;               ///< Flux for PSF at (x,y) for normalization = 1.0
+
+    float chisq;                        ///< PSF goodness statistic (unused??)
     float ApResid;                      ///< apMag - psfMag (for PSF stars)
     float dApResid;                     ///< scatter of ApResid
     float skyBias;                      ///< implied residual sky offset from ApResid fit
     float skySat;                       ///< roll-over of ApResid fit
-    float chisq;                        ///< PSF goodness statistic (unused??)
     int nPSFstars;                      ///< number of stars used to measure PSF
     int nApResid;                       ///< number of stars used to measure ApResid
+
+    bool poissonErrorsPhotLMM;          ///< use poission errors for non-linear model fitting
+    bool poissonErrorsPhotLin;          ///< use poission errors for linear model fitting
+    bool poissonErrorsParams;           ///< use poission errors for model parameter fitting
+
+    pmTrend2D *ApTrend;                 ///< ApResid vs (x,y)
+    pmTrend2D *FluxScale;               ///< Flux for PSF at (x,y) for normalization = 1.0
+    psPolynomial1D *ChiTrend;           ///< Chisq vs flux fit (correction for systematic errors)
+
+    pmGrowthCurve *growth;              ///< apMag vs Radius
+    pmResiduals *residuals;             ///< normalized residual image (no spatial variation)
+
+    psArray *params;                    ///< Model parameters (psPolynomial2D)
+    psStats *psfTrendStats;             ///< psf parameter trend clipping stats
+
+    pmTrend2DMode psfTrendMode;
     int trendNx;
     int trendNy;
@@ -50 +61 @@
     int fieldXo;
     int fieldYo;
-    bool poissonErrorsPhotLMM;          ///< use poission errors for non-linear model fitting
-    bool poissonErrorsPhotLin;          ///< use poission errors for linear model fitting
-    bool poissonErrorsParams;           ///< use poission errors for model parameter fitting
-    pmGrowthCurve *growth;              ///< apMag vs Radius
-    pmResiduals *residuals;             ///< normalized residual image (no spatial variation)
 };
 
@@ -60 +66 @@
     pmModelType   type;
     psStats      *stats;                // psfTrend clipping stats
+
     pmTrend2DMode psfTrendMode;
     int           psfTrendNx;
@@ -67 +74 @@
     int           psfFieldXo;
     int           psfFieldYo;
+
     bool          poissonErrorsPhotLMM; ///< use poission errors for non-linear model fitting
     bool          poissonErrorsPhotLin; ///< use poission errors for linear model fitting
     bool          poissonErrorsParams; ///< use poission errors for model parameter fitting
+
+    bool          chiFluxTrend;         // Fit a trend in Chi2 as a function of flux?
+    pmSourceFitOptions *fitOptions;
+
     float         fitRadius;
     float         apRadius;
-    bool          chiFluxTrend;         // Fit a trend in Chi2 as a function of flux?
-    pmSourceFitOptions *fitOptions;
 } pmPSFOptions;
 
@@ -94 +104 @@
 double pmPSF_SXYtoModel (psF32 *fittedPar);
 
-bool pmGrowthCurveGenerate (pmReadout *readout, pmPSF *psf, bool ignore, psImageMaskType maskVal, psImageMaskType mark);
 pmPSF *pmPSFBuildSimple (char *typeName, float sxx, float syy, float sxy, ...);
 

trunk/psModules/src/objects/pmPSF_IO.c

@@ -63 +63 @@
 
 bool pmPSFmodelReadPSFClump (psMetadata *analysis, psMetadata *header);
+bool pmPSFmodelRead_ApTrend (pmPSF *psf, pmFPAfile *file);
+bool pmPSFmodelWrite_ApTrend (pmFPAfile *file, pmPSF *psf);
+
+bool pmPSFmodelRead_GrowthCurve (pmPSF *psf, pmFPAfile *file);
+bool pmPSFmodelWrite_GrowthCurve (pmFPAfile *file, pmPSF *psf);
 
 bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file)
@@ -197 +202 @@
         psError(psErrorCodeLast(), false, "Failed to write PSF for chip");
         return false;
+    }
+    return true;
+}
+
+// XXX we save the model term identifiers (item) as S32, but they probably should be more flexible
+bool pmTrend2DtoTable (psArray *table, pmTrend2D *trend, char *label, int item) {
+
+    if (trend == NULL) return true; 
+
+    if (trend->mode == PM_TREND_MAP) {
+        // write the image components into a table: this is needed because they may each be a different size
+        psImageMap *map = trend->map;
+        for (int ix = 0; ix < map->map->numCols; ix++) {
+            for (int iy = 0; iy < map->map->numRows; iy++) {
+                psMetadata *row = psMetadataAlloc ();
+                psMetadataAddS32 (row, PS_LIST_TAIL, label,        0, "", item);
+                psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
+                psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
+                psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
+
+                psArrayAdd (table, 100, row);
+                psFree (row);
+            }
+        }
+    } else {
+        // write the polynomial components into a table
+        psPolynomial2D *poly = trend->poly;
+        for (int ix = 0; ix <= poly->nX; ix++) {
+            for (int iy = 0; iy <= poly->nY; iy++) {
+                psMetadata *row = psMetadataAlloc ();
+                psMetadataAddS32 (row, PS_LIST_TAIL, label,        0, "", item);
+                psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
+                psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
+                psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->coeffMask[ix][iy]);
+
+                psArrayAdd (table, 100, row);
+                psFree (row);
+            }
+        }
+    }
+    return true;
+}
+
+// extra trend2D elements from a row
+bool pmTrend2DfromTableRow (pmTrend2D *trend, psMetadata *row) {
+
+    bool status = false;
+
+    int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
+    int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
+
+    if (trend->mode == PM_TREND_MAP) {
+        psImageMap *map = trend->map;
+        assert (map);
+        assert (map->map);
+        assert (map->error);
+        assert (xPow >= 0);
+        assert (yPow >= 0);
+        assert (xPow < map->map->numCols);
+        assert (yPow < map->map->numRows);
+        map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
+        map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
+    } else {
+        psPolynomial2D *poly = trend->poly;
+        assert (poly);
+        assert (xPow >= 0);
+        assert (yPow >= 0);
+        assert (xPow <= poly->nX);
+        assert (yPow <= poly->nY);
+        poly->coeff[xPow][yPow]     = psMetadataLookupF32 (&status, row, "VALUE");
+        poly->coeffErr[xPow][yPow]  = psMetadataLookupF32 (&status, row, "ERROR");
+        poly->coeffMask[xPow][yPow] = psMetadataLookupU8  (&status, row, "MASK");
     }
     return true;
@@ -403 +484 @@
 
             if (trend->mode == PM_TREND_MAP) {
-              nX = trend->map->map->numCols;
-              nY = trend->map->map->numRows;
+                nX = trend->map->map->numCols;
+                nY = trend->map->map->numRows;
             } else {
-              nX = trend->poly->nX;
-              nY = trend->poly->nY;
+                nX = trend->poly->nX;
+                nY = trend->poly->nY;
             }
             snprintf (name, 9, "PAR%02d_NX", i);
@@ -428 +509 @@
         psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
 
+        float PSF_APERTURE =  psMetadataLookupF32(&status, roAnalysis, "PSF_APERTURE");
+        if (status) {
+            psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_APERTURE", PS_DATA_F32, "aperture for psf objects", PSF_APERTURE);
+        }
+        float PSF_FIT_RADIUS =  psMetadataLookupF32(&status, roAnalysis, "PSF_FIT_RADIUS");
+        if (status) {
+            psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_DATA_F32, "aperture for psf objects", PSF_FIT_RADIUS);
+        }
+
         // build a FITS table of the PSF parameters
         psArray *psfTable = psArrayAllocEmpty (100);
         for (int i = 0; i < nPar; i++) {
             pmTrend2D *trend = psf->params->data[i];
-            if (trend == NULL) continue; // skip unset parameters (eg, XPOS)
-
-            if (trend->mode == PM_TREND_MAP) {
-                // write the image components into a table: this is needed because they may each be a different size
-                psImageMap *map = trend->map;
-                for (int ix = 0; ix < map->map->numCols; ix++) {
-                    for (int iy = 0; iy < map->map->numRows; iy++) {
-                        psMetadata *row = psMetadataAlloc ();
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
-                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
-                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
-                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
-
-                        psArrayAdd (psfTable, 100, row);
-                        psFree (row);
-                    }
-                }
-            } else {
-                // write the polynomial components into a table
-                psPolynomial2D *poly = trend->poly;
-                for (int ix = 0; ix <= poly->nX; ix++) {
-                    for (int iy = 0; iy <= poly->nY; iy++) {
-                        psMetadata *row = psMetadataAlloc ();
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
-                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
-                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
-                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
-                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->coeffMask[ix][iy]);
-
-                        psArrayAdd (psfTable, 100, row);
-                        psFree (row);
-                    }
-                }
-            }
+            pmTrend2DtoTable (psfTable, trend, "MODEL_TERM", i);
         }
 
@@ -547 +601 @@
     }
 
+    if (!pmPSFmodelWrite_ApTrend(file, psf)) {
+        psError(psErrorCodeLast(), false, "Unable to write PSF ApTrend");
+        return false;
+    }
+
+    if (!pmPSFmodelWrite_GrowthCurve(file, psf)) {
+        psError(psErrorCodeLast(), false, "Unable to write PSF Growth Curve");
+        return false;
+    }
+
     // write a representation of the psf model
     {
         psMetadata *header = psMetadataAlloc ();
-
-        if (0) {
-            // set some header keywords to make it clear there are no residuals?
-            if (!psFitsWriteBlank (file->fits, header, residName)) {
-                psError(psErrorCodeLast(), false, "Unable to write blank PSF residuals.");
-                psFree(residName);
-                psFree(header);
-                return false;
-            }
-            psFree (residName);
-            psFree (header);
-            return true;
-        }
 
         int DX = 65;
@@ -651 +702 @@
 }
 
+
+
 // if this file needs to have a PHU written out, write one
 bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config)
@@ -934 +987 @@
     psf->skyBias   = psMetadataLookupF32 (&status, header, "SKY_BIAS");
 
+    if (roAnalysis) {
+        float PSF_APERTURE =  psMetadataLookupF32(&status, header, "PSF_APERTURE");
+        if (status) {
+            psMetadataAddF32 (roAnalysis, PS_LIST_TAIL, "PSF_APERTURE", PS_DATA_F32, "aperture for psf objects", PSF_APERTURE);
+        }
+        float PSF_FIT_RADIUS =  psMetadataLookupF32(&status, header, "PSF_FIT_RADIUS");
+        if (status) {
+            psMetadataAddF32 (roAnalysis, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_DATA_F32, "aperture for psf objects", PSF_FIT_RADIUS);
+        }
+    } else {
+        psWarning ("unable to read PSF_APERTURE or PSF_FIT_RADIUS");
+    }
+
     // read the raw table data
     psArray *table = psFitsReadTable (file->fits);
@@ -945 +1011 @@
     for (int i = 0; i < table->n; i++) {
         psMetadata *row = table->data[i];
+
         int iPar = psMetadataLookupS32 (&status, row, "MODEL_TERM");
-        int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
-        int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
 
         pmTrend2D *trend = psf->params->data[iPar];
@@ -955 +1020 @@
         }
 
-        if (trend->mode == PM_TREND_MAP) {
-            psImageMap *map = trend->map;
-            assert (map);
-            assert (map->map);
-            assert (map->error);
-            assert (xPow >= 0);
-            assert (yPow >= 0);
-            assert (xPow < map->map->numCols);
-            assert (yPow < map->map->numRows);
-            map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
-            map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
-        } else {
-            psPolynomial2D *poly = trend->poly;
-            assert (poly);
-            assert (xPow >= 0);
-            assert (yPow >= 0);
-            assert (xPow <= poly->nX);
-            assert (yPow <= poly->nY);
-            poly->coeff[xPow][yPow]     = psMetadataLookupF32 (&status, row, "VALUE");
-            poly->coeffErr[xPow][yPow]  = psMetadataLookupF32 (&status, row, "ERROR");
-            poly->coeffMask[xPow][yPow] = psMetadataLookupU8  (&status, row, "MASK");
-        }
+        pmTrend2DfromTableRow(trend, row);
     }
     psFree (header);
@@ -1063 +1107 @@
     }
 
+    if (!pmPSFmodelRead_ApTrend (psf, file)) {
+        psError(psErrorCodeLast(), false, "Unable to read PSF ApTrend data.");
+        return false;
+    }
+
+    if (!pmPSFmodelRead_GrowthCurve(psf, file)) {
+        psError(psErrorCodeLast(), false, "Unable to read PSF Growth Curve");
+        return false;
+    }
+
     psMetadataAdd (chipAnalysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
     psFree (psf);
@@ -1070 +1124 @@
     psFree (header);
 
+    return true;
+}
+
+// write aperture trend to a FITS table
+bool pmPSFmodelWrite_ApTrend (pmFPAfile *file, pmPSF *psf) {
+
+    pmTrend2D *trend = psf->ApTrend;
+    if (trend == NULL) { 
+        psWarning ("no PSF ApTrend to write out, skipping");
+        return true; 
+    }
+
+    // we need to write a header for the table,
+    psMetadata *header = psMetadataAlloc();
+
+    int nX = 0, nY = 0;
+    if (trend->mode == PM_TREND_MAP) {
+        nX = trend->map->map->numCols;
+        nY = trend->map->map->numRows;
+    } else {
+        nX = trend->poly->nX;
+        nY = trend->poly->nY;
+    }
+    psMetadataAddS32 (header, PS_LIST_TAIL, "TREND_NX", 0, "", nX);
+    psMetadataAddS32 (header, PS_LIST_TAIL, "TREND_NY", 0, "", nY);
+    char *modeName = pmTrend2DModeToString (trend->mode);
+    psMetadataAddStr (header, PS_LIST_TAIL, "TREND_MD", 0, "", modeName);
+    psFree (modeName);
+
+    // build a FITS table of the ApTrend (only 1)
+    psArray *table = psArrayAllocEmpty (100);
+    pmTrend2DtoTable (table, trend, "APTREND", 0);
+
+    // write an empty FITS segment if we have no PSF information
+    if (table->n == 0) {
+        psError(PM_ERR_PROG, true, "No PSF data to write.");
+        psFree(table);
+        psFree(header);
+        return false;
+    } 
+
+    psTrace ("pmFPAfile", 5, "writing psf ApTrend data %s\n", "AP_TREND");
+    if (!psFitsWriteTable(file->fits, header, table, "AP_TREND")) {
+        psError(psErrorCodeLast(), false, "Error writing psf table data %s\n", "AP_TREND");
+        psFree(table);
+        psFree(header);
+        return false;
+    }
+
+    psFree (table);
+    psFree (header);
+    return true;
+}
+
+// read aperture trend to a FITS table
+bool pmPSFmodelRead_ApTrend (pmPSF *psf, pmFPAfile *file) {
+
+    bool status;
+
+    // move fits pointer to AP_TREND section
+    // advance to the table data extension
+    if (!psFitsMoveExtNameClean (file->fits, "AP_TREND")) {
+        psWarning ("no Aperture Trend data in PSF file, skipping");
+        return true;
+    }
+
+    psMetadata *header = psFitsReadHeader (NULL, file->fits);
+    if (!header) {
+        psError(psErrorCodeLast(), false, "Unable to read AP_TREND header.");
+        return false;
+    }
+        
+    // read the raw table data
+    psArray *table = psFitsReadTable (file->fits);
+    if (!table) {
+        psError(psErrorCodeLast(), false, "Unable to read AP_TREND table.");
+        psFree(header);
+        return false;
+    }
+
+    // XXX allow user to set this optionally?
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    psImageBinning *binning = psImageBinningAlloc();
+    binning->nXfine = psf->fieldNx;
+    binning->nYfine = psf->fieldNy;
+    binning->nXruff = psMetadataLookupS32 (&status, header, "TREND_NX");
+    binning->nYruff = psMetadataLookupS32 (&status, header, "TREND_NY");
+    psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
+    char *modeName  = psMetadataLookupStr (&status, header, "TREND_MD");
+    if (!status) {
+        psError(PM_ERR_PROG, true, "inconsistent PSF header: NX & NY defined for AP TREND, but not MD");
+        psFree (header);
+        psFree (stats);
+        psFree (table);
+        return false;
+    }
+    pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName);
+    if (psfTrendMode == PM_TREND_NONE) {
+        psfTrendMode = PM_TREND_POLY_ORD;
+    }
+
+    // measure Trend2D for the current spatial scale
+    pmTrend2D *apTrend = pmTrend2DNoImageAlloc (PM_TREND_MAP, binning, stats);
+
+    // fill in the matching psf->params entries
+    for (int i = 0; i < table->n; i++) {
+        psMetadata *row = table->data[i];
+        pmTrend2DfromTableRow(apTrend, row);
+    }
+    psf->ApTrend = apTrend;
+
+    psFree (binning);
+    psFree (header);
+    psFree (stats);
+    psFree (table);
+    return true;
+}
+
+// write aperture trend to a FITS table
+bool pmPSFmodelWrite_GrowthCurve (pmFPAfile *file, pmPSF *psf) {
+
+    pmGrowthCurve *growth = psf->growth;
+    if (growth == NULL) { 
+        psWarning ("no PSF Growth Curve to write out, skipping");
+        return true; 
+    }
+
+    // we need to write a header for the table,
+    psMetadata *header = psMetadataAlloc();
+
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_MIN_RAD", 0, "", growth->radius->data.F32[0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_MAX_RAD", 0, "", growth->maxRadius);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_REF_RAD", 0, "", growth->refRadius);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_AP_LOSS", 0, "", growth->apLoss);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_AP_REF",  0, "", growth->apRef);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "GROWTH_FIT_MAG", 0, "", growth->fitMag);
+
+    // build a FITS table of the ApTrend (only 1)
+    psArray *table = psArrayAllocEmpty (100);
+    for (int i = 0; i < growth->apMag->n; i++) {
+        psMetadata *row = psMetadataAlloc ();
+        psMetadataAddF32 (row, PS_LIST_TAIL, "RADIUS", 0, "", growth->radius->data.F32[i]);
+        psMetadataAddF32 (row, PS_LIST_TAIL, "AP_MAG", 0, "", growth->apMag->data.F32[i]);
+        psArrayAdd (table, 100, row);
+        psFree (row);
+    }
+
+    // write an empty FITS segment if we have no PSF information
+    if (table->n == 0) {
+        psError(PM_ERR_PROG, true, "No PSF data to write.");
+        psFree(table);
+        psFree(header);
+        return false;
+    } 
+
+    psTrace ("pmFPAfile", 5, "writing psf Growth Curve data %s\n", "GROWTH_CURVE");
+    if (!psFitsWriteTable(file->fits, header, table, "GROWTH_CURVE")) {
+        psError(psErrorCodeLast(), false, "Error writing psf table data %s\n", "GROWTH_CURVE");
+        psFree(table);
+        psFree(header);
+        return false;
+    }
+
+    psFree (table);
+    psFree (header);
+    return true;
+}
+
+// read aperture trend to a FITS table
+bool pmPSFmodelRead_GrowthCurve (pmPSF *psf, pmFPAfile *file) {
+
+    bool status;
+
+    // move fits pointer to AP_TREND section
+    // advance to the table data extension
+    if (!psFitsMoveExtNameClean (file->fits, "GROWTH_CURVE")) {
+        psWarning ("no Growth Curve data in PSF file, skipping");
+        return true;
+    }
+
+    psMetadata *header = psFitsReadHeader (NULL, file->fits);
+    if (!header) {
+        psError(psErrorCodeLast(), false, "Unable to read GROWTH_CURVE header.");
+        return false;
+    }
+        
+    // read the raw table data
+    psArray *table = psFitsReadTable (file->fits);
+    if (!table) {
+        psError(psErrorCodeLast(), false, "Unable to read GROWTH_CURVE table.");
+        psFree(header);
+        return false;
+    }
+
+    float minRadius = psMetadataLookupF32 (&status, header, "GROWTH_MIN_RAD"); if (!status) return false;
+    float maxRadius = psMetadataLookupF32 (&status, header, "GROWTH_MAX_RAD"); if (!status) return false;
+    float refRadius = psMetadataLookupF32 (&status, header, "GROWTH_REF_RAD"); if (!status) return false;
+
+    psf->growth = pmGrowthCurveAlloc(minRadius, maxRadius, refRadius);
+
+    psf->growth->apLoss = psMetadataLookupF32 (&status, header, "GROWTH_AP_LOSS"); if (!status) return false;
+    psf->growth->apRef  = psMetadataLookupF32 (&status, header, "GROWTH_AP_REF"); if (!status) return false;
+    psf->growth->fitMag = psMetadataLookupF32 (&status, header, "GROWTH_FIT_MAG"); if (!status) return false;
+
+    // fill in the matching psf->params entries
+    for (int i = 0; i < table->n; i++) {
+        psMetadata *row = table->data[i];
+        psf->growth->apMag->data.F32[i] = psMetadataLookupF32 (&status, row, "AP_MAG"); if (!status) return false;
+    }
+
+    psFree (header);
+    psFree (table);
     return true;
 }

trunk/psModules/src/objects/pmPSFtryFitPSF.c

@@ -122 +122 @@
         psfTry->fitMag->data.F32[i] = source->psfMag;
         psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
-        psfTry->metricErr->data.F32[i] = source->errMag;
+        psfTry->metricErr->data.F32[i] = source->psfMagErr;
 
         // XXX this did not work: modifies shape of psf too much
@@ -131 +131 @@
                  source->peak->xf, source->peak->yf, 
                  source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], 
-                 source->psfMag, source->apMag, source->errMag,
+                 source->psfMag, source->apMag, source->psfMagErr,
                  source->modelPSF->params->data.F32[PM_PAR_I0], 
                  source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], 

trunk/psModules/src/objects/pmPeaks.c

@@ -150 +150 @@
     tmp->x = x;
     tmp->y = y;
+    tmp->xf = x;
+    tmp->yf = y;
+    tmp->dx = NAN;
+    tmp->dy = NAN;
     tmp->detValue        = value;
     tmp->rawFlux         = value; // set this by default: it is up to the user to supply a better value
@@ -155 +159 @@
     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;
     tmp->assigned = false;
     tmp->type = type;
@@ -166 +167 @@
     psTrace("psModules.objects", 10, "---- %s() end ----\n", __func__);
     return(tmp);
+}
+
+// copy to an already allocated peak
+bool pmPeakCopy(pmPeak *out, pmPeak *in)
+{
+    out->x               = in->x;
+    out->y               = in->y;
+    out->xf              = in->xf;
+    out->yf              = in->yf;
+    out->dx              = in->dx;
+    out->dy              = in->dy;
+    out->detValue        = in->detValue;
+    out->rawFlux         = in->rawFlux;
+    out->rawFluxStdev    = in->rawFluxStdev;
+    out->smoothFlux      = in->smoothFlux;
+    out->smoothFluxStdev = in->smoothFluxStdev;
+    out->assigned        = in->assigned;
+    out->type            = in->type;
+
+    return true;
 }
 

trunk/psModules/src/objects/pmPeaks.h

@@ -69 +69 @@
     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.
@@ -89 +88 @@
 
 bool psMemCheckPeak(psPtr ptr);
+
+bool pmPeakCopy(pmPeak *out, pmPeak *in);
 
 /** pmPeaksInVector()

trunk/psModules/src/objects/pmSource.c

@@ -116 +116 @@
     source->psfImage = NULL;
     source->moments = NULL;
-    source->blends = NULL;
     source->modelPSF = NULL;
     source->modelEXT = NULL;
@@ -124 +123 @@
     source->mode2 = PM_SOURCE_MODE_DEFAULT;
     source->tmpFlags = 0;
-    source->extpars = NULL;
-    source->diffStats = NULL;
-    source->radialAper = NULL;
-    source->parent = NULL;
-
-    source->region = psRegionSet(NAN, NAN, NAN, NAN);
-    psMemSetDeallocator(source, (psFreeFunc) sourceFree);
 
     // default values are NAN
     source->psfMag           = NAN;
+    source->psfMagErr        = NAN;
     source->psfFlux          = NAN;
     source->psfFluxErr       = NAN;
     source->extMag           = NAN;    
-    source->errMag           = NAN;
     source->apMag            = NAN;
     source->apMagRaw         = NAN;
@@ -143 +135 @@
     source->apFlux           = NAN;
     source->apFluxErr        = NAN; 
-    source->sky              = NAN;
-    source->skyErr           = NAN;    
+
     source->pixWeightNotBad  = NAN;
     source->pixWeightNotPoor = NAN;
@@ -151 +142 @@
     source->crNsigma         = NAN;
     source->extNsigma        = NAN;
+    source->sky              = NAN;
+    source->skyErr           = NAN;    
+
+    source->region = psRegionSet(NAN, NAN, NAN, NAN);
+    source->blends = NULL;
+    source->extpars = NULL;
+    source->diffStats = NULL;
+    source->radialAper = NULL;
+    source->parent = NULL;
+    source->imageID = -1;
+
+    psMemSetDeallocator(source, (psFreeFunc) sourceFree);
 
     psTrace("psModules.objects", 10, "---- end ----\n");
@@ -172 +175 @@
     if (in->peak != NULL) {
         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->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;
+        pmPeakCopy(source->peak, in->peak);
     }
 
@@ -195 +192 @@
 
     // the maskObj is a unique mask array; create a new mask image
-    source->maskObj = in->maskObj ? psImageCopy (NULL, in->maskObj, PS_TYPE_IMAGE_MASK) : NULL;
-
-    source->type = in->type;
-    source->mode = in->mode;
-    source->imageID = in->imageID;
+    source->maskObj = in->maskObj   ? psImageCopy (NULL, in->maskObj, PS_TYPE_IMAGE_MASK) : NULL;
+
+    // NOTE : because of the const id element, we cannot just assign *source = *in
+
+    source->type             = in->type;
+    source->mode             = in->mode;
+    source->mode2            = in->mode2;
+    source->tmpFlags         = in->tmpFlags;
+    source->psfMag           = in->psfMag;
+    source->psfMagErr        = in->psfMagErr;
+    source->psfFlux          = in->psfFlux;
+    source->psfFluxErr       = in->psfFluxErr;
+    source->extMag           = in->extMag;
+    source->apMag            = in->apMag;
+    source->apMagRaw         = in->apMagRaw;
+    source->apRadius         = in->apRadius;
+    source->apFlux           = in->apFlux;
+    source->apFluxErr        = in->apFluxErr;
+    source->pixWeightNotBad  = in->pixWeightNotBad;
+    source->pixWeightNotPoor = in->pixWeightNotPoor;
+    source->psfChisq         = in->psfChisq;
+    source->crNsigma         = in->crNsigma;
+    source->extNsigma        = in->extNsigma;
+    source->sky              = in->sky;
+    source->skyErr           = in->skyErr;
+
+    source->region           = in->region;
 
     return(source);
@@ -266 +285 @@
     extend |= (mySource->maskObj == NULL);
     extend |= (mySource->maskView == NULL);
+
+    // if ((fabs(x - 2020) < 5) && (fabs(y - 366) < 5)) {
+    //  if (extend) {
+    //      fprintf (stderr, "extend T, %f, %f : %f, %f vs %f, %f : %f, %f\n", 
+    //               newRegion.x0, newRegion.y0, newRegion.x1, newRegion.y1, 
+    //               mySource->region.x0, mySource->region.y0, mySource->region.x1, mySource->region.y1);
+    //  } else {
+    //      fprintf (stderr, "extend F, %f, %f : %f, %f vs %f, %f : %f, %f\n", 
+    //               newRegion.x0, newRegion.y0, newRegion.x1, newRegion.y1, 
+    //               mySource->region.x0, mySource->region.y0, mySource->region.x1, mySource->region.y1);
+    //  }
+    // }
 
     if (extend) {
@@ -486 +517 @@
         // create vectors with Sx, Sy values in window
         // clip sources based on S/N
-        for (psS32 i = 0 ; i < sources->n ; i++)
+        for (psS32 i = 0; i < sources->n; i++)
         {
             pmSource *tmpSrc = (pmSource *) sources->data[i];
@@ -1175 +1206 @@
     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_MOMENTS_FAILURE) return false; // can't if the moments failed...
 
     if (source->mode & PM_SOURCE_MODE_SATSTAR) return true; // moments are best for SATSTARs

trunk/psModules/src/objects/pmSource.h

@@ -80 +80 @@
     pmSourceMode2 mode2;                ///< analysis flags set for object.
     pmSourceTmpF tmpFlags;              ///< internal-only flags
-    psArray *blends;                    ///< collection of sources thought to be confused with object
+
     float psfMag;                       ///< calculated from flux in modelPSF
+    float psfMagErr;                    ///< error in psfMag
     float psfFlux;                      ///< calculated from flux in modelPSF
-    float psfFluxErr;                   ///< calculated from flux in modelPSF
-    float extMag;                       ///< calculated from flux in modelEXT
-    float errMag;                       ///< error in psfMag OR extMag (depending on type)
+    float psfFluxErr;                   ///< error in psfFlux
+    float extMag;                       ///< calculated from flux in modelEXT -- NOTE this is not actually used
     float apMag;                        ///< apMag corresponding to psfMag or extMag (depending on type)
     float apMagRaw;                     ///< raw mag in given aperture
@@ -98 +98 @@
     float crNsigma;                     ///< Nsigma deviation from PSF to CR
     float extNsigma;                    ///< Nsigma deviation from PSF to EXT
-    float sky, skyErr;                  ///< The sky and its error at the center of the object
+    float sky;                          ///< The sky at the center of the object 
+    float skyErr;                       ///< The sky error at the center of the object
+
     psRegion region;                    ///< area on image covered by selected pixels
+    psArray *blends;                    ///< collection of sources thought to be confused with object
     pmSourceExtendedPars *extpars;      ///< extended source parameters
     pmSourceDiffStats *diffStats;       ///< extra parameters for difference detections
@@ -254 +257 @@
     float sigma,      ///< size of Gaussian window function (<= 0.0 -> skip window)
     float minSN,              ///< minimum pixel significance
+    float minKronRadius,      ///< minimum pixel significance
     psImageMaskType maskVal
 );
@@ -273 +277 @@
   float xGuess, float yGuess);
 
+float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM);
+
 pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source);
 

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c

@@ -125 +125 @@
         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);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
         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);
@@ -256 +256 @@
 
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
         PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
         pmPSF_AxesToModel (PAR, axes);

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c

@@ -124 +124 @@
         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);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
         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);
@@ -274 +274 @@
 
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
         PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
         pmPSF_AxesToModel (PAR, axes);

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

@@ -125 +125 @@
         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);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
         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);
@@ -273 +273 @@
         // XXX use these to determine PAR[PM_PAR_I0]?
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
         PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
         pmPSF_AxesToModel (PAR, axes);
@@ -720 +720 @@
       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]);
+      psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
+      psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
     }
 

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c

@@ -118 +118 @@
         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);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
         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",              outputs.apRadius);
@@ -241 +241 @@
         // XXX use these to determine PAR[PM_PAR_I0]?
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
         PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
         pmPSF_AxesToModel (PAR, axes);

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c

@@ -116 +116 @@
         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);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
         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",              outputs.apRadius);
@@ -247 +247 @@
         // XXX use these to determine PAR[PM_PAR_I0]?
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
         PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
         pmPSF_AxesToModel (PAR, axes);

trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V3.c

@@ -116 +116 @@
         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);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
         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);
@@ -161 +161 @@
         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);
+
+        // Do NOT write these : not consistent with the definition of PS1_V3 in Ohana/src/libautocode/dev/cmf-ps1-v3.d
         // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_FLUX",   PS_DATA_F32, "Kron Flux (in 1.0 R1)",                      moments.KronCore);
-        // psMetadataAdd (row, PS_LIST_TAIL, "KRON_CORE_ERROR",  PS_DATA_F32, "Kron Error (in 1.0 R1)",                     moments.KronCoreErr);
+        // 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);
@@ -270 +272 @@
         // XXX use these to determine PAR[PM_PAR_I0]?
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
         PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
         pmPSF_AxesToModel (PAR, axes);
@@ -292 +294 @@
 
         // we no longer sort by S/N, only flux
-        // if (isfinite (source->errMag) && (source->errMag > 0.0)) {
-        //   source->peak->SN = 1.0 / source->errMag;
+        // if (isfinite (source->psfMagErr) && (source->psfMagErr > 0.0)) {
+        //   source->peak->SN = 1.0 / source->psfMagErr;
         // } else {
         //   source->peak->SN = sqrt(source->peak->flux); // an alternate proxy: various functions sort by peak S/N

trunk/psModules/src/objects/pmSourceIO_CMP.c

@@ -137 +137 @@
         axes = pmPSF_ModelToAxes (PAR, 20.0);
 
-        float errMag = isfinite(source->errMag) ? source->errMag : 999;
+        float psfMagErr = isfinite(source->psfMagErr) ? source->psfMagErr : 999;
 
         psLineInit (line);
@@ -143 +143 @@
         psLineAdd (line, "%6.1f ",  PAR[PM_PAR_YPOS]);
         psLineAdd (line, "%6.3f ",  PS_MIN (99.0, source->psfMag + ZERO_POINT));
-        psLineAdd (line, "%03d ",   PS_MIN (999, (int)errMag));
+        psLineAdd (line, "%03d ",   PS_MIN (999, (int)psfMagErr));
         psLineAdd (line, "%2d ",    type);
         psLineAdd (line, "%3.1f ",  lsky);
@@ -280 +280 @@
             source->psfMag = atof (array->data[2]);
             source->extMag = atof (array->data[6]);
-            source->errMag = atof (array->data[3]) / 1000.0;
+            source->psfMagErr = atof (array->data[3]) / 1000.0;
             source->apMag  = atof (array->data[7]);
             axes.major     = atof (array->data[8]);

trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c

@@ -71 +71 @@
     psF32 xPos, yPos;
     psF32 xErr, yErr;
-    psF32 errMag, chisq;
+    psF32 psfMagErr, chisq;
 
     pmChip *chip = readout->parent->parent;
@@ -121 +121 @@
             chisq = model->chisq;
 
-            // need to determine the PSF photometry error: source->errMag is the error on the 'best' model mag.
-            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
+            // need to determine the PSF photometry error: source->psfMagErr is the error on the 'best' model mag.
+            psfMagErr = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
         } else {
             xPos = source->peak->xf;
@@ -132 +132 @@
             axes.theta = NAN;
             chisq = NAN;
-            errMag = NAN;
+            psfMagErr = NAN;
         }
 
@@ -167 +167 @@
         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr);
         psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "Calibrated Magnitude from PSF Fit",          calMag);
-        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "Calibrated Magnitude Error",                 errMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "Calibrated Magnitude Error",                 psfMagErr);
         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",        psfMagErr);
         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
@@ -284 +284 @@
         // XXX use these to determine PAR[PM_PAR_I0]?
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
         pmPSF_AxesToModel (PAR, axes);

trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c

@@ -113 +113 @@
         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr);
         psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             PS_MIN (99.0, source->psfMag));
-        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        PS_MIN (99.0, source->errMag));
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        PS_MIN (99.0, source->psfMagErr));
         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           PS_MIN (99.0, peakMag));
         psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
@@ -212 +212 @@
         // XXX use these to determine PAR[PM_PAR_I0]?
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
 
         pmPSF_AxesToModel (PAR, axes);

trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

@@ -70 +70 @@
     psF32 xPos, yPos;
     psF32 xErr, yErr;
-    psF32 errMag, chisq;
+    psF32 psfMagErr, chisq;
 
     // let's write these out in S/N order
@@ -103 +103 @@
             chisq = model->chisq;
 
-            // need to determine the PSF photometry error: source->errMag is the error on the 'best' model mag.
-            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
+            // need to determine the PSF photometry error: source->psfMagErr is the error on the 'best' model mag.
+            psfMagErr = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
         } else {
             xPos = source->peak->xf;
@@ -114 +114 @@
             axes.theta = NAN;
             chisq = NAN;
-            errMag = NAN;
+            psfMagErr = NAN;
         }
 
@@ -128 +128 @@
         psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr);
         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",        psfMagErr);
         psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
         psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
@@ -253 +253 @@
         // XXX use these to determine PAR[PM_PAR_I0]?
         source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
-        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->psfMagErr    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
-        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
 
         pmPSF_AxesToModel (PAR, axes);

trunk/psModules/src/objects/pmSourceIO_RAW.c

@@ -110 +110 @@
         fprintf (f, "%7.1f %7.1f  %7.1f %8.4f  %7.4f %7.4f  ",
                  PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], source->sky,
-                 source->psfMag, source->errMag, dPos);
+                 source->psfMag, source->psfMagErr, dPos);
 
         for (j = 4; j < model->params->n; j++) {
@@ -174 +174 @@
         fprintf (f, "%7.1f %7.1f  %7.1f %8.4f  %7.4f %7.4f  ",
                  PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], source->sky,
-                 source->extMag, source->errMag, dPos);
+                 source->extMag, source->psfMagErr, dPos);
 
         for (j = 4; j < model->params->n; j++) {

trunk/psModules/src/objects/pmSourceIO_SMPDATA.c

@@ -106 +106 @@
         psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX",   PS_DATA_F32, "", yPos);
         psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", PS_MIN (99.0, source->psfMag + ZERO_POINT));
-        psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", PS_MIN (999, 1000*source->errMag));
+        psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", PS_MIN (999, 1000*source->psfMagErr));
         psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", PS_MIN (99.0, source->extMag + ZERO_POINT));
         psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP",  PS_DATA_F32, "", PS_MIN (99.0, source->apMag + ZERO_POINT));
@@ -194 +194 @@
         source->psfMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT;
         source->extMag = psMetadataLookupF32 (&status, row, "MAG_GAL") - ZERO_POINT;
-        source->errMag = psMetadataLookupF32 (&status, row, "MAG_ERR") * 0.001;
+        source->psfMagErr = psMetadataLookupF32 (&status, row, "MAG_ERR") * 0.001;
         source->apMag  = psMetadataLookupF32 (&status, row, "MAG_AP")  - ZERO_POINT;
 

trunk/psModules/src/objects/pmSourceIO_SX.c

@@ -88 +88 @@
         psLineAdd (line, "%11.3f", PAR[PM_PAR_YPOS]);
         psLineAdd (line, "%9.4f",  source->psfMag);
-        psLineAdd (line, "%9.4f",  source->errMag);
+        psLineAdd (line, "%9.4f",  source->psfMagErr);
         psLineAdd (line, "%13.4f", source->sky);
         psLineAdd (line, "%9.2f",  axes.major);

trunk/psModules/src/objects/pmSourceMasks.h

@@ -14 +14 @@
     PM_SOURCE_MODE_PSFSTAR          = 0x00000040, ///< Source used to define PSF model
     PM_SOURCE_MODE_SATSTAR          = 0x00000080, ///< Source model peak is above saturation
-    PM_SOURCE_MODE_BLEND            = 0x00000100, ///< Source is a blend with other sourcers
+    PM_SOURCE_MODE_BLEND            = 0x00000100, ///< Source is a blend with other sources
     PM_SOURCE_MODE_EXTERNAL         = 0x00000200, ///< Source based on supplied input position
     PM_SOURCE_MODE_BADPSF           = 0x00000400, ///< Failed to get good estimate of object's PSF
@@ -51 +51 @@
     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
+
+    PM_SOURCE_MODE2_PASS1_SRC        = 0x00000080, ///< source detected in first pass analysis
 } pmSourceMode2;
 

trunk/psModules/src/objects/pmSourceMatch.c

@@ -88 +88 @@
         pmSource *source = sources->data[i]; // Source of interest
         if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) ||
-            !isfinite(source->errMag) || source->psfMag > SOURCE_FAINTEST) {
+            !isfinite(source->psfMagErr) || source->psfMag > SOURCE_FAINTEST) {
             continue;
         }
@@ -102 +102 @@
         (*y)->data.F32[num] = ySrc;
         (*mag)->data.F32[num] = source->psfMag;
-        (*magErr)->data.F32[num] = source->errMag;
+        (*magErr)->data.F32[num] = source->psfMagErr;
         (*indices)->data.S32[num] = i;
         num++;

trunk/psModules/src/objects/pmSourceMoments.c

@@ -66 +66 @@
 // if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
 
-bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
+bool pmSourceMoments(pmSource *source, float radius, float sigma, float minSN, float minKronRadius, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -74 +74 @@
 
     // this function assumes the sky has been well-subtracted for the image
-    psF32 sky = 0.0;
+    float sky = 0.0;
 
     if (source->moments == NULL) {
@@ -80 +80 @@
     }
 
-    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);
-    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
+    float Sum = 0.0;
+    float Var = 0.0;
+    float SumCore = 0.0;
+    float VarCore = 0.0;
+    float R2 = PS_SQR(radius);
+    float minSN2 = PS_SQR(minSN);
+    float rsigma2 = 0.5 / PS_SQR(sigma);
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -109 +109 @@
     // Xn  = SUM (x - xc)^n * (z - sky)
 
-    psF32 RF = 0.0;
-    psF32 RH = 0.0;
-    psF32 RS = 0.0;
-    psF32 XX = 0.0;
-    psF32 XY = 0.0;
-    psF32 YY = 0.0;
-    psF32 XXX = 0.0;
-    psF32 XXY = 0.0;
-    psF32 XYY = 0.0;
-    psF32 YYY = 0.0;
-    psF32 XXXX = 0.0;
-    psF32 XXXY = 0.0;
-    psF32 XXYY = 0.0;
-    psF32 XYYY = 0.0;
-    psF32 YYYY = 0.0;
+    float RFW = 0.0;
+    float RHW = 0.0;
+
+    float RF = 0.0;
+    float RH = 0.0;
+    float RS = 0.0;
+    float XX = 0.0;
+    float XY = 0.0;
+    float YY = 0.0;
+    float XXX = 0.0;
+    float XXY = 0.0;
+    float XYY = 0.0;
+    float YYY = 0.0;
+    float XXXX = 0.0;
+    float XXXY = 0.0;
+    float XXYY = 0.0;
+    float XYYY = 0.0;
+    float YYYY = 0.0;
 
     Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
+
+    // float dX = source->moments->Mx - source->peak->xf;
+    // float dY = source->moments->My - source->peak->yf;
+    // float dR = hypot(dX, dY);
+    // float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
+    // float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
+    float Xo = source->moments->Mx;
+    float Yo = source->moments->My;
 
     // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     // xCM, yCM from pixel coords to pixel index here.
-    psF32 xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
-    psF32 yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+    float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
+    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
-        psF32 yDiff = row - yCM;
+        float yDiff = row - yCM;
         if (fabs(yDiff) > radius) continue;
 
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
+        float *vPix = source->pixels->data.F32[row];
+        float *vWgt = source->variance->data.F32[row];
+
         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+        // psImageMaskType  *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
 
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -151 +164 @@
             if (isnan(*vPix)) continue;
 
-            psF32 xDiff = col - xCM;
+            float xDiff = col - xCM;
             if (fabs(xDiff) > radius) continue;
 
             // radius is just a function of (xDiff, yDiff)
-            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
             if (r2 > R2) continue;
 
-            psF32 fDiff = *vPix - sky;
-            psF32 pDiff = fDiff;
-            psF32 wDiff = *vWgt;
+            float fDiff = *vPix - sky;
+            float pDiff = fDiff;
+            float wDiff = *vWgt;
 
             // skip pixels below specified significance level.  this is allowed, but should be
@@ -171 +184 @@
             // weighting over weights the sky for faint sources
             if (sigma > 0.0) {
-                psF32 z = r2 * rsigma2;
+                float z = r2 * rsigma2;
                 assert (z >= 0.0);
-                psF32 weight  = exp(-z);
+                float weight  = exp(-z);
 
                 wDiff *= weight;
@@ -182 +195 @@
 
             // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
-            psF32 r = sqrt(r2);
-            psF32 rf = r * fDiff;
-            psF32 rh = sqrt(r) * fDiff;
-            psF32 rs = fDiff;
-
-            psF32 x = xDiff * pDiff;
-            psF32 y = yDiff * pDiff;
-
-            psF32 xx = xDiff * x;
-            psF32 xy = xDiff * y;
-            psF32 yy = yDiff * y;
-
-            psF32 xxx = xDiff * xx / r;
-            psF32 xxy = xDiff * xy / r;
-            psF32 xyy = xDiff * yy / r;
-            psF32 yyy = yDiff * yy / r;
-
-            psF32 xxxx = xDiff * xxx / r2;
-            psF32 xxxy = xDiff * xxy / r2;
-            psF32 xxyy = xDiff * xyy / r2;
-            psF32 xyyy = xDiff * yyy / r2;
-            psF32 yyyy = yDiff * yyy / r2;
+            float r = sqrt(r2);
+            float rf = r * fDiff;
+            float rh = sqrt(r) * fDiff;
+            float rs = fDiff;
+
+            float rfw = r * pDiff;
+            float rhw = sqrt(r) * pDiff;
+
+            float x = xDiff * pDiff;
+            float y = yDiff * pDiff;
+
+            float xx = xDiff * x;
+            float xy = xDiff * y;
+            float yy = yDiff * y;
+
+            float xxx = xDiff * xx / r;
+            float xxy = xDiff * xy / r;
+            float xyy = xDiff * yy / r;
+            float yyy = yDiff * yy / r;
+
+            float xxxx = xDiff * xxx / r2;
+            float xxxy = xDiff * xxy / r2;
+            float xxyy = xDiff * xyy / r2;
+            float xyyy = xDiff * yyy / r2;
+            float yyyy = yDiff * yyy / r2;
 
             RF  += rf;
             RH  += rh;
             RS  += rs;
+
+            RFW  += rfw;
+            RHW  += rhw;
 
             XX  += xx;
@@ -244 +263 @@
     source->moments->Myyyy = YYYY/Sum;
 
-    // Calculate the Kron magnitude (make this block optional?)
-    float radKinner = 1.0*source->moments->Mrf;
-    float radKron   = 2.5*source->moments->Mrf;
-    float radKouter = 4.0*source->moments->Mrf;
+    // if Mrf (first radial moment) is very small, we are getting into low-significance
+    // territory.  saturate at minKronRadius.  conversely, if Mrf is > radius, we are clearly
+    // making an error.  saturate at radius.
+    float kronRefRadius = MIN(radius, MAX(minKronRadius, source->moments->Mrf));
+
+    float radKinner = 1.0*kronRefRadius;
+    float radKron   = 2.5*kronRefRadius;
+    float radKouter = 4.0*kronRefRadius;
 
     int nKronPix = 0;
@@ -259 +282 @@
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
-        psF32 yDiff = row - yCM;
+        float yDiff = row - yCM;
         if (fabs(yDiff) > radKouter) continue;
 
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
+        float *vPix = source->pixels->data.F32[row];
+        float *vWgt = source->variance->data.F32[row];
+
         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+        // psImageMaskType  *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
 
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -276 +301 @@
             if (isnan(*vPix)) continue;
 
-            psF32 xDiff = col - xCM;
+            float xDiff = col - xCM;
             if (fabs(xDiff) > radKouter) continue;
 
             // radKron is just a function of (xDiff, yDiff)
-            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
-
-            psF32 pDiff = *vPix - sky;
-            psF32 wDiff = *vWgt;
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+
+            float pDiff = *vPix - sky;
+            float wDiff = *vWgt;
 
             // skip pixels below specified significance level.  this is allowed, but should be
@@ -290 +315 @@
             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
 
-            psF32 r  = sqrt(r2);
+            float r  = sqrt(r2);
             if (r < radKron) {
                 Sum += pDiff;
@@ -335 +360 @@
 }
 
-bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
+bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) { 
 
     // First Pass: calculate the first moments (these are subtracted from the coordinates below)
@@ -342 +367 @@
     // .. etc
 
-    psF32 sky = 0.0;
-
-    psF32 peakPixel = -PS_MAX_F32;
+    float sky = 0.0;
+
+    float peakPixel = -PS_MAX_F32;
     psS32 numPixels = 0;
-    psF32 Sum = 0.0;
-    psF32 Var = 0.0;
-    psF32 X1 = 0.0;
-    psF32 Y1 = 0.0;
-    psF32 R2 = PS_SQR(radius);
-    psF32 minSN2 = PS_SQR(minSN);
-    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
+    float Sum = 0.0;
+    float Var = 0.0;
+    float X1 = 0.0;
+    float Y1 = 0.0;
+    float R2 = PS_SQR(radius);
+    float minSN2 = PS_SQR(minSN);
+    float rsigma2 = 0.5 / PS_SQR(sigma);
 
     float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
@@ -358 +383 @@
 
     // we are guaranteed to have a valid pixel and variance at this location (right? right?)
-    // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
+    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
@@ -370 +395 @@
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
 
-        psF32 yDiff = row + 0.5 - yPeak;
+        float yDiff = row + 0.5 - yPeak;
         if (fabs(yDiff) > radius) continue;
 
-        psF32 *vPix = source->pixels->data.F32[row];
-        psF32 *vWgt = source->variance->data.F32[row];
+        float *vPix = source->pixels->data.F32[row];
+        float *vWgt = source->variance->data.F32[row];
+
         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
+        // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
 
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -387 +414 @@
             if (isnan(*vPix)) continue;
 
-            psF32 xDiff = col + 0.5 - xPeak;
+            float xDiff = col + 0.5 - xPeak;
             if (fabs(xDiff) > radius) continue;
 
             // radius is just a function of (xDiff, yDiff)
-            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
             if (r2 > R2) continue;
 
-            psF32 pDiff = *vPix - sky;
-            psF32 wDiff = *vWgt;
+            float pDiff = *vPix - sky;
+            float wDiff = *vWgt;
 
             // skip pixels below specified significance level.  for a PSFs, this
@@ -408 +435 @@
             // weighting over weights the sky for faint sources
             if (sigma > 0.0) {
-                psF32 z  = r2*rsigma2;
+                float z  = r2*rsigma2;
                 assert (z >= 0.0);
-                psF32 weight  = exp(-z);
+                float weight  = exp(-z);
 
                 wDiff *= weight;
@@ -419 +446 @@
             Sum += pDiff;
 
-            psF32 xWght = xDiff * pDiff;
-            psF32 yWght = yDiff * pDiff;
+            float xWght = xDiff * pDiff;
+            float yWght = yDiff * pDiff;
 
             X1  += xWght;
@@ -477 +504 @@
     return true;
 }
+
+float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) {
+
+    psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *src = sources->data[i]; // Source of interest
+        if (!src || !src->moments) {
+            continue;
+        }
+
+        if (src->mode & PM_SOURCE_MODE_BLEND) {
+            continue;
+        }
+
+        if (!src->moments->nPixels) continue;
+
+        if (src->moments->SN < PSF_SN_LIM) continue;
+
+        // XXX put in Mxx,Myy cut based on clump location
+
+        psVectorAppend(radii, src->moments->Mrf);
+    }
+
+    // find the peak in this image
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    if (!psVectorStats (stats, radii, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
+        psFree(stats);
+        return NAN;
+    }
+
+    float minRadius = stats->sampleMedian;
+
+    psFree(radii);
+    psFree(stats);
+    return minRadius;
+}
+

trunk/psModules/src/objects/pmSourceOutputs.c

@@ -98 +98 @@
         outputs->xPos = PAR[PM_PAR_XPOS];
         outputs->yPos = PAR[PM_PAR_YPOS];
-        if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
+        if ((source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) && !(source->mode & PM_SOURCE_MODE_EXTMODEL)) {
+            // we only do non-linear PSF fits for non-extended objects
             outputs->xErr = dPAR[PM_PAR_XPOS];
             outputs->yErr = dPAR[PM_PAR_YPOS];
         } else {
-            outputs->xErr = fwhmMajor * source->errMag / 2.35;
-            outputs->yErr = fwhmMinor * source->errMag / 2.35;
+            outputs->xErr = fwhmMajor * source->psfMagErr / 2.35;
+            outputs->yErr = fwhmMinor * source->psfMagErr / 2.35;
         }
-        if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
+        if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXY]) && isfinite(PAR[PM_PAR_SYY])) {
             axes = pmPSF_ModelToAxes (PAR, 20.0);
             outputs->psfMajor = axes.major;
@@ -125 +126 @@
             outputs->xPos = source->moments->Mx;
             outputs->yPos = source->moments->My;
-            outputs->xErr = fwhmMajor * source->errMag / 2.35;
-            outputs->yErr = fwhmMinor * source->errMag / 2.35;
+            outputs->xErr = fwhmMajor * source->psfMagErr / 2.35;
+            outputs->yErr = fwhmMinor * source->psfMagErr / 2.35;
         } else {
             outputs->xPos = source->peak->xf;

trunk/psModules/src/objects/pmSourcePhotometry.c

@@ -82 +82 @@
    - apMag  : only if S/N > AP_MIN_SN
    : is optionally corrected for curve-of-growth if:
-   - the source is a STAR (PSF)
    - the option is selected (mode & PM_SOURCE_PHOT_GROWTH)
    - psfMag : all sources with non-NULL modelPSF
    : is optionally corrected for aperture residual if:
-   - the source is a STAR (PSF)
    - the option is selected (mode & PM_SOURCE_PHOT_APCORR)
-
    - extMag : all sources with non-NULL modelEXT
 **/
@@ -100 +97 @@
 
     int status = false;
-    bool isPSF;
     float x, y;
-    float rflux;
     float SN;
-    pmModel *model;
 
     source->psfMag    = NAN;
     source->extMag    = NAN;
-    source->errMag    = NAN;
+    source->psfMagErr = NAN;
     source->apMag     = NAN;
     source->apMagRaw  = NAN;
@@ -114 +108 @@
     source->apFluxErr = NAN;
 
-    // we must have a valid model
-    // XXX allow aperture magnitudes for sources without a model
-    model = pmSourceGetModel (&isPSF, source);
-    if (model == NULL) {
-        psTrace ("psModules.objects", 3, "fail mag : no valid model");
+    // XXXXXX review:
+    // Select the 'best' model -- this is used for PSF_QF,_PERFECT & ???. isPSF is true if this
+    // object is a PSF (not extended).  We must have a valid model.  XXX NOTE: allow aperture
+    // magnitudes for sources without a model
+
+    // select the psf model
+    pmModel *modelPSF = source->modelPSF;
+    if (modelPSF == NULL) {
+        psTrace ("psModules.objects", 3, "fail mag : no valid PSF model");
         return false;
     }
 
-    // XXX handle negative flux, low-significance
-    if (model->dparams->data.F32[PM_PAR_I0] > 0) {
-        SN = fabs(model->params->data.F32[PM_PAR_I0] / model->dparams->data.F32[PM_PAR_I0]);
-        source->errMag = 1.0 / SN;
+    // get the error on the PSF model magnitude
+    if (modelPSF->dparams->data.F32[PM_PAR_I0] > 0) {
+        SN = fabs(modelPSF->params->data.F32[PM_PAR_I0] / modelPSF->dparams->data.F32[PM_PAR_I0]);
+        source->psfMagErr = 1.0 / SN;
     } else {
         SN = NAN;
-        source->errMag = NAN;
-    }
-    x = model->params->data.F32[PM_PAR_XPOS];
-    y = model->params->data.F32[PM_PAR_YPOS];
+        source->psfMagErr = NAN;
+    }
+    // the source position is used to recenter the aperture for ap photometry
+    x = modelPSF->params->data.F32[PM_PAR_XPOS];
+    y = modelPSF->params->data.F32[PM_PAR_YPOS];
 
     // measure PSF model photometry
-    // XXX TEST: do not use flux scale
-    // XXX NOTE: turn this back on?
-    if (0 && psf->FluxScale) {
-        // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
-        double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
-        psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
-        psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
-        source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
-        source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
-        source->psfMag = -2.5*log10(source->psfFlux);
-    } else {
-        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
-        source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
-    }
+    status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, modelPSF);
+    source->psfFluxErr = source->psfFlux * source->psfMagErr;
+
+# if (0)
+    // XXX NOTE: old code to use the flux scale.  test & turn this back on?  if so, need to save with psf model
+    // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
+    double fluxScale = pmTrend2DEval (psf->FluxScale, (float)source->peak->x, (float)source->peak->y);
+    psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
+    psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
+    source->psfFlux = fluxScale * modelPSF->params->data.F32[PM_PAR_I0];
+    source->psfFluxErr = fluxScale * modelPSF->dparams->data.F32[PM_PAR_I0];
+    source->psfMag = -2.5*log10(source->psfFlux);
+# endif
 
     if (mode == PM_SOURCE_PHOT_PSFONLY) {
@@ -153 +151 @@
     }
 
+    // get the EXT model photometry (all EXT models)
     // if we have a collection of model fits, check if one of them is a pointer to modelEXT
     if (source->modelFits) {
@@ -173 +172 @@
     }
 
-    // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
-    // XXX add a flag for "ap_mag is corrected?"
-    if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) {
+    // Correct psfMag to match aperture magnitude system (NOTE : Growth curve is already applied to ApTrend)
+    if ((mode & PM_SOURCE_PHOT_APCORR) && psf && psf->ApTrend) {
         // the source peak pixel is guaranteed to be on the image, and only minimally different from the source center
         double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y);
         source->psfMag += apTrend;
-    }
-
-    // measure the contribution of included pixels
+        source->psfFlux *= pow(10.0, -0.4*apTrend);
+        source->psfFluxErr *= pow(10.0, -0.4*apTrend);
+    }
+
+    // measure the contribution of included pixels to the PSF model fit
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (source, model, source->maskObj, maskVal, radius);
+        pmSourcePixelWeight (source, modelPSF, source->maskObj, maskVal, radius);
     }
 
@@ -217 +217 @@
         y += dy;
 
-        if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy,
-                              NAN, 0xff, PS_INTERPOLATE_BIQUADRATIC)) {
+        // if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_LANCZOS2)) {
+        // if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_BIQUADRATIC)) {
+        if (!psImageShiftMask(&flux, &mask, source->pixels, source->maskObj, maskVal, dx, dy, NAN, 0xff, PS_INTERPOLATE_BILINEAR)) {
             // Not much we can do about it
             psErrorClear();
@@ -232 +233 @@
 
     // measure object aperture photometry
-    status = pmSourcePhotometryAperSource (source, model, flux, variance, mask, maskVal);
+    status = pmSourcePhotometryAperSource (source, modelPSF, flux, variance, mask, maskVal);
     if (!status) {
         psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag");
@@ -241 +242 @@
     // detection limits (esp near bright neighbors)
     source->apMag = source->apMagRaw;
-    if (isfinite (source->apMag) && isPSF && psf) {
+    if (isfinite (source->apMag) && psf) {
         if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) {
-            source->apMag = source->apMagRaw + pmGrowthCurveCorrect (psf->growth, source->apRadius);
-            // XXX correct the apFlux?
-        }
-        if (mode & PM_SOURCE_PHOT_APCORR) {
-            // XXX this should be removed -- we no longer fit for the 'sky bias'
-            // XXX is this happening???
-            rflux   = pow (10.0, 0.4*source->psfMag);
-            psAssert (psf->skyBias == 0.0, "sky bias not 0");
-            psAssert (psf->skySat == 0.0, "sky sat not 0");
-            source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux;
+            float apOffset = pmGrowthCurveCorrect (psf->growth, source->apRadius);
+            source->apMag = source->apMagRaw + apOffset;
+            source->apFlux *= pow(10.0, -0.4*apOffset);
+            source->apFluxErr *= pow(10.0, -0.4*apOffset);
         }
     }
@@ -309 +304 @@
 
     bool status;
-    status = pmSourcePhotometryAper(&source->apMagRaw, &source->apFlux, &source->apFluxErr, model, image, variance, mask, maskVal);
-
+    int nPix = 0;
+    status = pmSourcePhotometryAper(&nPix, &source->apMagRaw, &source->apFlux, &source->apFluxErr, model, image, variance, mask, maskVal);
+    if (status) {
+        source->mode |= PM_SOURCE_MODE_AP_MAGS;
+    }
     return status;
 }
 
 // return source aperture magnitude
-bool pmSourcePhotometryAper (float *apMag, float *apFluxOut, float *apFluxErr, pmModel *model, psImage *image, psImage *variance, psImage *mask, psImageMaskType maskVal)
+bool pmSourcePhotometryAper (int *nPixOut, float *apMag, float *apFluxOut, float *apFluxErr, pmModel *model, psImage *image, psImage *variance, psImage *mask, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(apMag, false);
@@ -328 +326 @@
     float apFlux = 0;
     float apFluxVar = 0;
+    int nPix = 0;
 
     if (DO_SKY) {
@@ -345 +344 @@
             apFlux += imData[iy][ix] - sky;
             apFluxVar += varData[iy][ix];
-        }
-    }
+            nPix ++;
+        }
+    }
+    
     if (apFluxOut) *apFluxOut = apFlux;
     if (apFluxErr) *apFluxErr = sqrt(fabs(apFluxVar));
+    if (nPixOut) *nPixOut = nPix;
 
     if (apFlux <= 0) {
@@ -412 +414 @@
     psImageMaskType maskBad = maskVal;
     maskBad &= ~maskSuspect;
+
+    psImageMaskType maskPoor = maskVal | maskSuspect;
 
     // measure modelSum and validSum.  this function is applied to a sources' subimage.  the
@@ -446 +450 @@
 
             // count pixels which are masked with an mask bit (bad or poor)
-            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskVal)) {
+            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[my][mx] & maskPoor)) {
                 notPoorSum += value;
             }
@@ -487 +491 @@
 
     if (isfinite(source->pixWeightNotBad) && isfinite(source->pixWeightNotPoor)) {
-        psAssert (source->pixWeightNotBad <= source->pixWeightNotPoor, "error: all bad pixels should also be poor");
+        psAssert (source->pixWeightNotPoor <= source->pixWeightNotBad, "error: all bad pixels should also be poor");
     }
 

trunk/psModules/src/objects/pmSourcePhotometry.h

@@ -45 +45 @@
 
 bool pmSourcePhotometryAper(
+    int *nPixOut,
     float *apMag,
     float *apFluxOut,

trunk/psModules/src/objects/pmSourceUtils.h

@@ -38 +38 @@
 
 /// @}
-# endif /* PM_MODEL_CLASS_H */
+# endif /* PM_SOURCE_UTILS_H */

trunk/psModules/src/psmodules.h

@@ -148 +148 @@
 #include <pmModelUtils.h>
 #include <pmSourcePhotometry.h>
+#include <pmGrowthCurveGenerate.h>
 #include <pmSourceVisual.h>
 #include <pmSourceMatch.h>

trunk/psModules/test/objects/tap_pmSource.c

@@ -91 +91 @@
         ok(isnan(src->psfMag), "pmSourceAlloc() pmSource->psfMag correctly");
         ok(isnan(src->extMag), "pmSourceAlloc() pmSource->extMag correctly");
-        ok(isnan(src->errMag), "pmSourceAlloc() pmSource->errMag correctly");
+        ok(isnan(src->psfMagErr), "pmSourceAlloc() pmSource->psfMagErr correctly");
         ok(isnan(src->apMag), "pmSourceAlloc() pmSource->apMag correctly");
         ok(isnan(src->sky), "pmSourceAlloc() pmSource->sky correctly");

trunk/psModules/test/objects/tap_pmSourceIO_PS1_DEV_0.c

@@ -167 +167 @@
             model->dparams->data.F32[PM_PAR_YPOS] = TEST_BASE_Y_ERR + (float) i;
             src->psfMag = TEST_BASE_PSF_MAG + (float) i;
-            src->errMag = TEST_BASE_ERR_MAG + (float) i;
+            src->psfMagErr = TEST_BASE_ERR_MAG + (float) i;
             src->peak->flux = TEST_BASE_PEAK_FLUX + (float) i;
             src->sky = TEST_BASE_SKY + (float) i;
@@ -215 +215 @@
              ok(src->psfMag == (TEST_BASE_PSF_MAG + (float) i), "pmSourcesRead_PS1_DEV_0() set src->psfMag correctly (is %.2f, should be %.2f)",
                 src->psfMag, (TEST_BASE_PSF_MAG + (float) i));
-             ok(src->errMag == (TEST_BASE_ERR_MAG + (float) i), "pmSourcesRead_PS1_DEV_0() set src->errMag correctly (is %.2f, should be %.2f)",
-                src->errMag, (TEST_BASE_ERR_MAG + (float) i));
+             ok(src->psfMagErr == (TEST_BASE_ERR_MAG + (float) i), "pmSourcesRead_PS1_DEV_0() set src->psfMagErr correctly (is %.2f, should be %.2f)",
+                src->psfMagErr, (TEST_BASE_ERR_MAG + (float) i));
 
              // XXX: Source code always sets src->modelPSF.  Is that right?

trunk/psModules/test/objects/tap_pmSourceIO_PS1_DEV_1.c

@@ -196 +196 @@
             model->dparams->data.F32[PM_PAR_YPOS] = TEST_BASE_Y_ERR + (float) i;
             src->psfMag = TEST_BASE_PSF_MAG + (float) i;
-            src->errMag = TEST_BASE_ERR_MAG + (float) i;
+            src->psfMagErr = TEST_BASE_ERR_MAG + (float) i;
             src->peak->flux = TEST_BASE_PEAK_FLUX + (float) i;
             src->sky = TEST_BASE_SKY + (float) i;
@@ -203 +203 @@
             sources->data[i] = (psPtr *) src;
             src->psfMag = TEST_BASE_PSF_MAG + (float) i;
-            src->errMag = TEST_BASE_ERR_MAG + (float) i;
+            src->psfMagErr = TEST_BASE_ERR_MAG + (float) i;
             src->sky = TEST_BASE_SKY + (float) i;
             src->skyErr = TEST_BASE_SKY_ERR + (float) i;
@@ -280 +280 @@
              ok(src->psfMag == (TEST_BASE_PSF_MAG + (float) i), "pmSourcesRead_PS1_DEV_1() set src->psfMag correctly (is %.2f, should be %.2f)",
                 src->psfMag, (TEST_BASE_PSF_MAG + (float) i));
-             ok(src->errMag == (TEST_BASE_ERR_MAG + (float) i), "pmSourcesRead_PS1_DEV_1() set src->errMag correctly (is %.2f, should be %.2f)",
-                src->errMag, (TEST_BASE_ERR_MAG + (float) i));
+             ok(src->psfMagErr == (TEST_BASE_ERR_MAG + (float) i), "pmSourcesRead_PS1_DEV_1() set src->psfMagErr correctly (is %.2f, should be %.2f)",
+                src->psfMagErr, (TEST_BASE_ERR_MAG + (float) i));
 
              // XXX: Source code always sets src->modelPSF.  Is that right?

trunk/psModules/test/objects/tap_pmSourceIO_SMPDATA.c

@@ -167 +167 @@
             model->dparams->data.F32[PM_PAR_YPOS] = TEST_BASE_Y_ERR + (float) i;
             src->psfMag = TEST_BASE_PSF_MAG + (float) i;
-            src->errMag = TEST_BASE_ERR_MAG + (float) i;
+            src->psfMagErr = TEST_BASE_ERR_MAG + (float) i;
             src->peak->flux = TEST_BASE_PEAK_FLUX + (float) i;
             src->sky = TEST_BASE_SKY + (float) i;
@@ -226 +226 @@
                 src->psfMag, (TEST_BASE_PSF_MAG + (float) i));
              float tmpF =  0.001 * PS_MIN(999, (1000 * (TEST_BASE_ERR_MAG + (float) i)));
-             ok(src->errMag == tmpF, "pmSourcesRead_SMPDATA() set src->errMag correctly (is %.2f, should be %.2f)",
-                src->errMag, tmpF);
+             ok(src->psfMagErr == tmpF, "pmSourcesRead_SMPDATA() set src->psfMagErr correctly (is %.2f, should be %.2f)",
+                src->psfMagErr, tmpF);
              tmpF = PS_MIN(99.0, (TEST_BASE_EXT_MAG + ZERO_POINT)) - ZERO_POINT;
              ok(src->extMag == tmpF, "pmSourcesRead_SMPDATA() set src->extMag correctly (is %.2f, should be %.2f)",

trunk/psModules/test/objects/tap_pmSourcePhotometry.c

@@ -103 +103 @@
     ok_float_tol(source->psfMag, fitMag, 0.0002, "source fitMag is %f", source->psfMag);
     ok_float_tol(source->apMag,  apMag, 0.0002, "source apMag is %f", source->apMag);
-    ok_float(source->errMag,   0.001, "source errMag is %f", source->errMag);
+    ok_float(source->psfMagErr,   0.001, "source psfMagErr is %f", source->psfMagErr);
     float refMag = source->apMag;