Changeset 33638 for branches/eam_branches/ipp-20111122/psModules/src/objects

Timestamp:

Mar 30, 2012, 2:49:37 PM (14 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ipp-20111122/psModules/src/objects

Files:

• models/pmModel_QGAUSS.c (modified) (1 diff)
• models/pmModel_SERSIC.c (modified) (1 diff)
• pmFootprintCullPeaks.c (modified) (1 diff)
• pmModelUtils.c (modified) (1 diff)
• pmPCMdata.c (modified) (1 diff)
• pmPSFtryFitEXT.c (modified) (2 diffs)
• pmSourceIO.c (modified) (4 diffs)
• pmSourceIO.h (modified) (1 diff)
• pmSourceIO_CMF_PS1_SV1.c (modified) (3 diffs)
• pmSourcePhotometry.c (modified) (1 diff)

branches/eam_branches/ipp-20111122/psModules/src/objects/models/pmModel_QGAUSS.c

@@ -402 +402 @@
     assert (psf->params->n > PM_PAR_YPOS);
     assert (psf->params->n > PM_PAR_XPOS);
+
+    if (! isfinite(Io)) {
+        fprintf(stderr, "non-finite Io passed to PM_MODEL_PARAMS_FROM_PSF\n");
+        return false;
+    }
 
     PAR[PM_PAR_SKY]  = 0.0;

branches/eam_branches/ipp-20111122/psModules/src/objects/models/pmModel_SERSIC.c

@@ -192 +192 @@
     psF32 z0 = PAR[PM_PAR_I0]*f1;
     psF32 f0 = PAR[PM_PAR_SKY] + z0;
+
+    if (!isfinite(z0)) {
+        fprintf(stderr, "z0 is not finite for %f %f %f %f %f.  Parameters: \n", X, Y, radius, z, f1);
+        fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4],
+            PAR[5], PAR[6], PAR[7]);
+    }
 
     assert (isfinite(f2));

branches/eam_branches/ipp-20111122/psModules/src/objects/pmFootprintCullPeaks.c

@@ -178 +178 @@
             psArray *myFP = pmFootprintsFind(subImg, threshold, 5);
             if (!myFP) {
-                psWarning ("missing footprint?");
+                psWarning ("missing footprint? threshold: %.f", threshold);
                 continue;
             }
             if (!myFP->n) {
-                psWarning ("empty footprint?");
+                psWarning ("empty footprint? threshold: %.f", threshold);
                 psFree (myFP);
                 continue;

branches/eam_branches/ipp-20111122/psModules/src/objects/pmModelUtils.c

@@ -144 +144 @@
 
     *Io = source->peak->rawFlux;
+
+#ifndef ALLOW_NONFINITE_PEAK
+    // Gene says fail of peak !finite
+    if (!isfinite(*Io)) return false;
+#else 
+    // This is the way it used to be. Somtimes an infinite value Io made it's way down the pipeline
+    // causing assertion failures
     if (!isfinite(*Io) && !source->moments) return false;
 
     *Io = source->moments->Peak;
     if (!isfinite(*Io)) return false;
+#endif
 
     return true;

branches/eam_branches/ipp-20111122/psModules/src/objects/pmPCMdata.c

@@ -136 +136 @@
             sum += value;
         }
+    }
+
+    if (!(sum > 0.0)) {
+        // Crazy PSF image print out some debugging information ...
+        fprintf(stderr, "invalid kernel sum %f found by pmPCMkernelFromPSF\n", sum);    for (int j = psf->yMin; j <= psf->yMax; j++) {
+            fprintf(stderr, "Row %d\n", j);
+            for (int i = psf->xMin; i <= psf->xMax; i++) {
+                double value = source->psfImage->data.F32[y0 + j][x0 + i];
+                fprintf(stderr, "  %d %f\n", i, value);
+            }
+        }
+        fflush(stderr);
+        // ... but avoid the asssertion two lines down by escaping
+        goto escape;
     }
     assert (sum > 0.0);

branches/eam_branches/ipp-20111122/psModules/src/objects/pmPSFtryFitEXT.c

@@ -73 +73 @@
             continue;
         }
+        // If mask object does not exist, mark the source as bad.
+        // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes
+        // psphot to exit with a fault. 
+        if (source->maskObj == NULL) {
+            psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y);
+            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+            continue;
+        }
 
         source->modelEXT = pmSourceModelGuess (source, options->type);
@@ -89 +97 @@
 
         // clear object mask to define valid pixels
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
 
         // exclude the poor fits

branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO.c

@@ -57 +57 @@
 #define BLANK_HEADERS "BLANK.HEADERS"   // Name of metadata in camera configuration containing header names
                                         // for putting values into a blank PHU
+static bool pmReadoutReadXSRC(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex);
+static bool pmReadoutReadXFIT(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex);
+static bool pmReadoutReadXRAD(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex);
 
 // lookup the EXTNAME values used for table data and image header segments
@@ -961 +964 @@
         psString dataname = NULL;
         psString deteffname = NULL;
-        if (!pmSourceIOextnames(&headname, &dataname, &deteffname, NULL, NULL, NULL, file, view)) {
+        psString xsrcname = NULL;
+        psString xfitname = NULL;
+        psString xradname = NULL;
+
+        // determine the output table format. Assume if we need to output extendend source
+        // parameters that they may exist in the input. 
+        // XXX: Perhaps we should use different recipe values.
+        // I.E. EXTENDED_SOURCE_ANALYSIS_READ or something like that
+        psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT");
+        if (!status) {
+            psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data");
+            return false;
+        }
+        // if this is not TRUE, the output files only contain the psf measurements.
+        bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS");
+        bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS");
+        bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES");
+
+        if (!pmSourceIOextnames(&headname, &dataname, &deteffname, 
+                XSRC_OUTPUT ? &xsrcname : NULL, 
+                XFIT_OUTPUT ? &xfitname : NULL, 
+                XRAD_OUTPUT ? &xradname : NULL,
+                file, view)) {
             return false;
         }
@@ -1039 +1064 @@
             }
 
+            long *sourceIndex = NULL;
+            if (XSRC_OUTPUT || XFIT_OUTPUT || XRAD_OUTPUT) {
+                long seq_max = -1;
+                for (long i = sources->n -1; i >= 0; i--) {
+                    pmSource *source = sources->data[i];
+                    if (source->seq < 0) {
+                        // This can happen cmf files that have been corrupted
+                        psError(PS_ERR_IO, true, "seq < 0 for source %ld: Suspect %s is corrupt", i, file->origname);
+                        return false;
+                    }
+                    if (source->seq > seq_max) {
+                        seq_max = source->seq;
+                    }
+                }
+                sourceIndex = psAlloc((seq_max + 1) * sizeof(long));
+                for (long i = 0; i < seq_max; i++) {
+                    sourceIndex[i] = -1;
+                }
+                for (long i = 0; i < sources->n; i++) {
+                    pmSource *source = sources->data[i];
+                    sourceIndex[source->seq] = i;
+                }
+            }
+            if (XSRC_OUTPUT && xsrcname) {
+                if (!pmReadoutReadXSRC(file, exttype, hdu->header, xsrcname, sources, sourceIndex)) {
+                    // XXX: is this an error?
+                    psErrorClear();
+                }
+                psFree(xsrcname);
+            }
+            if (XFIT_OUTPUT && xfitname) {
+                if (!pmReadoutReadXFIT(file, exttype, hdu->header, xfitname, sources, sourceIndex)) {
+                    // XXX: is this an error?
+                    psErrorClear();
+                }
+                psFree(xfitname);
+            }
+            if (XRAD_OUTPUT && xradname) {
+                if (!pmReadoutReadXRAD(file, readout, exttype, hdu->header, xradname, sources, sourceIndex)) {
+                    // XXX: is this an error?
+                    psErrorClear();
+                }
+                psFree(xradname);
+            }
+            psFree(sourceIndex);
+
             if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) {
 #if 0
@@ -1165 +1236 @@
 }
 
-
+// XXX: We might be able to macroize this and reuse for the other types
+
+static bool pmReadoutReadXSRC(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex) 
+{
+    if (!psFitsMoveExtName (file->fits, xsrcname)) {
+        psError(PS_ERR_UNKNOWN, false, "cannot find xsrc extension %s in %s", xsrcname, file->filename);
+        return false;
+    }
+
+    psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header
+    if (!tableHeader) psAbort("cannot read table header");
+
+    char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
+    if (!xtension) psAbort("cannot read table type");
+    if (strcmp (xtension, "BINTABLE")) {
+        psWarning ("no binary table in extension %s, skipping\n", xsrcname);
+        return false;
+    }
+
+    // XXX these are case-sensitive since the EXTYPE is case-sensitive
+    bool status = false;
+    if (file->type == PM_FPA_FILE_CMF) {
+        if (!strcmp (exttype, "PS1_SV1")) {
+            status  = pmSourcesRead_CMF_PS1_SV1_XSRC (file->fits, hduHeader, sources, sourceIndex);
+        }
+    }
+    psFree(tableHeader);
+    return status;
+}
+
+static bool pmReadoutReadXFIT(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString extname, psArray *sources, long *sourceIndex) 
+{
+    if (!psFitsMoveExtName (file->fits, extname)) {
+        psError(PS_ERR_UNKNOWN, false, "cannot find extension %s in %s", extname, file->filename);
+        return false;
+    }
+
+    psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header
+    if (!tableHeader) psAbort("cannot read table header");
+
+    char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
+    if (!xtension) psAbort("cannot read table type");
+    if (strcmp (xtension, "BINTABLE")) {
+        psWarning ("no binary table in extension %s, skipping\n", extname);
+        return false;
+    }
+
+    // XXX these are case-sensitive since the EXTYPE is case-sensitive
+    bool status = false;
+    if (file->type == PM_FPA_FILE_CMF) {
+        if (!strcmp (exttype, "PS1_SV1")) {
+            status  = pmSourcesRead_CMF_PS1_SV1_XFIT (file->fits, hduHeader, sources, sourceIndex);
+        }
+    }
+    psFree(tableHeader);
+    return status;
+}
+static bool pmReadoutReadXRAD(pmFPAfile *file, pmReadout *readout, char *exttype, psMetadata *hduHeader, psString extname, psArray *sources, long *sourceIndex) 
+{
+    if (!psFitsMoveExtName (file->fits, extname)) {
+        psError(PS_ERR_UNKNOWN, false, "cannot find extension %s in %s", extname, file->filename);
+        return false;
+    }
+
+    psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header
+    if (!tableHeader) psAbort("cannot read table header");
+
+    char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
+    if (!xtension) psAbort("cannot read table type");
+    if (strcmp (xtension, "BINTABLE")) {
+        psWarning ("no binary table in extension %s, skipping\n", extname);
+        return false;
+    }
+
+    // XXX these are case-sensitive since the EXTYPE is case-sensitive
+    bool status = false;
+    if (file->type == PM_FPA_FILE_CMF) {
+        if (!strcmp (exttype, "PS1_SV1")) {
+            status  = pmSourcesRead_CMF_PS1_SV1_XRAD (file->fits, readout, hduHeader, sources, sourceIndex);
+        }
+    }
+    psFree(tableHeader);
+    return status;
+}

branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO.h

@@ -93 +93 @@
 psArray *pmSourcesRead_CMF_PS1_V4 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header);
+bool pmSourcesRead_CMF_PS1_SV1_XSRC (psFits *fits, psMetadata *header, psArray *sources, long *);
+bool pmSourcesRead_CMF_PS1_SV1_XFIT (psFits *fits, psMetadata *header, psArray *sources, long *);
+bool pmSourcesRead_CMF_PS1_SV1_XRAD (psFits *fits, pmReadout *readout, psMetadata *header, psArray *sources, long *);
 psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_CMF_PS1_DV2 (psFits *fits, psMetadata *header);

branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c

@@ -539 +539 @@
 }
 
+bool pmSourcesRead_CMF_PS1_SV1_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(sources, false);
+
+    bool status;
+    long numSources = psFitsTableSize(fits); // Number of sources in table
+    if (numSources == 0) {
+        psError(psErrorCodeLast(), false, "XSRC Table contains no entries\n");
+        return false;
+    }
+
+    // petrosian mags are not saved, we need to calculate fluxes. For this we need exptime and zero point
+    float zeropt = psMetadataLookupF32(&status, hduHeader, "FPA.ZP");
+    float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME");
+    float magOffset = zeropt + 2.5*log10(exptime);
+
+    for (long i = 0; i < numSources; i++) {
+        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
+        if (!row) {
+            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
+            psFree(row);
+            return false;
+        }
+        // Find the source with this sequence number. 
+        // XXX: I am assuming that sources is sorted in order of seq
+        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
+        pmSource *source = NULL;
+#ifndef ASSUME_SORTED
+        long j = seq < sources->n ? seq : sources->n - 1;
+        for (; j >= 0; j--) {
+            source = sources->data[j];
+            if (source->seq == seq) {
+                break;
+            }
+        }
+#else
+        long j = sourceIndex[seq];
+        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
+        source = sources->data[j];
+#endif
+        if (!source) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
+            psFree(row);
+            return false;
+        }
+
+        if (!source->extpars) {
+            source->extpars = pmSourceExtendedParsAlloc ();
+        }
+        pmSourceExtendedPars *extpars = source->extpars;
+
+        // Assume that X_EXT Y_EXT and sigmas match the psf src so skip
+
+        // We don't have enough information to calculate the major and minor axis. Set major to 1. Should we scale this by
+        // psf size or something?
+        extpars->axes.major = 1.0;
+        extpars->axes.minor = extpars->axes.major * psMetadataLookupF32(&status, row, "F25_ARATIO");
+        extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA");
+
+        float mag = psMetadataLookupF32(&status, row, "PETRO_MAG");
+        float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR");
+        if (isfinite(mag)) {
+            extpars->petrosianFlux    = pow(10., (magOffset - mag) / 2.5);
+            if (isfinite(magErr)) {
+                extpars->petrosianFluxErr = extpars->petrosianFlux / magErr;
+            }
+        }
+
+        extpars->petrosianRadius   = psMetadataLookupF32(&status, row, "PETRO_RADIUS");
+        extpars->petrosianRadiusErr= psMetadataLookupF32(&status, row, "PETRO_RADIUS_ERR");
+        extpars->petrosianR50      = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50");
+        extpars->petrosianR50Err   = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50_ERR");
+        extpars->petrosianR90      = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90");
+        extpars->petrosianR90Err   = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90_ERR");
+        extpars->petrosianFill     = psMetadataLookupF32(&status, row, "PETRO_FILL");
+
+        psVector *radSB   = psMetadataLookupVector(&status, row, "PROF_SB");
+        psVector *radFlux = psMetadataLookupVector(&status, row, "PROF_FLUX");
+        psVector *radFill = psMetadataLookupVector(&status, row, "PROF_FILL");
+
+        if (radSB && radSB->n > 0) {
+            extpars->radProfile = pmSourceRadialProfileAlloc();
+            extpars->radProfile->binSB   = psMemIncrRefCounter(radSB);
+            extpars->radProfile->binSum   = psMemIncrRefCounter(radFlux);
+            extpars->radProfile->binFill = psMemIncrRefCounter(radFill);
+        }
+
+        psFree(row);
+    }
+
+    return true;
+}
+
 // XXX this layout is still the same as PS1_DEV_1
 bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname)
@@ -684 +778 @@
     psFree (outhead);
     psFree (table);
+    return true;
+}
+
+bool pmSourcesRead_CMF_PS1_SV1_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(sources, false);
+
+    bool status;
+    long numSources = psFitsTableSize(fits); // Number of sources in table
+    if (numSources == 0) {
+        psError(psErrorCodeLast(), false, "XFIT Table contains no entries\n");
+        return false;
+    }
+
+    for (long i = 0; i < numSources; i++) {
+        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
+        if (!row) {
+            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
+            psFree(row);
+            return false;
+        }
+        // Find the source with this sequence number. 
+        // XXX: I am assuming that sources is sorted in order of seq.
+        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
+        long j = seq < sources->n ? seq : sources->n - 1;
+        pmSource *source = NULL;
+        for (; j >= 0; j--) {
+            source = sources->data[j];
+            if (source->seq == seq) {
+                break;
+            }
+        }
+        if (!source) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
+            psFree(row);
+            return false;
+        }
+        if (!source->modelFits) {
+            // XXX: where to find the number of models to expect?
+            source->modelFits = psArrayAllocEmpty(5);
+        }
+        psString modelName = psMetadataLookupStr(&status, row, "MODEL_TYPE");
+        if (!modelName) {
+            psError(PS_ERR_UNKNOWN, true, "Failed to find model name for row %ld\n", i);
+            psFree(row);
+            return false;
+        }
+        pmModelType modelType = pmModelClassGetType(modelName);
+        if (modelType < 0) {
+            psError(PS_ERR_UNKNOWN, true, "Failed to find model type for %s\n", modelName);
+            psFree(row);
+            return false;
+        }
+        pmModel *model = pmModelAlloc(modelType);
+
+        psF32 *PAR = model->params->data.F32;
+        psF32 *dPAR = model->dparams->data.F32;
+
+        PAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT");
+        PAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT");
+        dPAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT_SIG");
+        dPAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT_SIG");
+
+        model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG");
+        model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG");
+
+        psEllipseAxes axes;
+        axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ");
+        axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN");
+        axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA");
+        if (!pmPSF_AxesToModel(PAR, axes, modelType)) {
+            // Do we need to fail here or can this happen?
+            psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model");
+            psFree(model);
+            psFree(row);
+            return false;
+        }
+        // XXX: clean this up
+        if (model->params->n > 7) {
+            PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07");
+        }
+        // read the covariance matrix
+        int nparams = model->params->n;
+        psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
+        for (int y = 0; y < nparams; y++) {
+            for (int x = 0; x < nparams; x++) {
+                char name[64];
+                snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
+                covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
+            }
+        }
+        model->covar = covar;
+
+        psArrayAdd(source->modelFits, 1, model);
+        psFree(model);
+
+        psFree(row);
+    }
+
     return true;
 }
@@ -831 +1025 @@
     return true;
 }
+
+bool pmSourcesRead_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(sources, false);
+
+    bool status;
+    long numSources = psFitsTableSize(fits); // Number of sources in table
+    if (numSources == 0) {
+        psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n");
+        return false;
+    }
+
+    long       seq_first = -1;
+    long       seq_last = -1;
+    psVector   *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
+    long       max_entries = -1;
+    long       num_entries = -1;
+
+    for (long i = 0; i < numSources; i++) {
+        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
+        if (!row) {
+            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
+            psFree(row);
+            return false;
+        }
+        // Find the source with this sequence number. 
+        // XXX: I am assuming that sources is sorted in order of seq.
+        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
+        long j = seq < sources->n ? seq : sources->n - 1;
+        pmSource *source = NULL;
+        for (; j >= 0; j--) {
+            source = sources->data[j];
+            if (source->seq == seq) {
+                break;
+            }
+        }
+        if (!source) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
+            psFree(row);
+            return false;
+        }
+        if (seq_first == -1) {
+            seq_first = seq;
+        }
+        if (seq == seq_first) {
+            psF32 value = psMetadataLookupF32(&status, row, "PSF_FWHM");
+            psVectorAppend(fwhmValues, value);
+        }
+        if (seq == seq_last) {
+            num_entries++;
+        } else {
+            num_entries = 1;
+            seq_last = seq;
+        }
+        if (num_entries > max_entries) {
+            max_entries = num_entries;
+        }
+
+        if (!source->radialAper) {
+            // XXX: where to find the number of models to expect?
+            source->radialAper = psArrayAllocEmpty(5);
+        }
+        pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc();
+
+        radialAper->flux = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX"));
+        radialAper->fluxStdev = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_STDEV"));
+        radialAper->fluxErr = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_ERR"));
+        radialAper->fill = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FILL"));
+
+        psArrayAdd(source->radialAper, 1, radialAper);
+
+        psFree(radialAper);
+        psFree(row);
+    }
+
+    // check for consistency between the length of fwhmValues and the maximum number of entries for each row
+    if (fwhmValues->n != max_entries) {
+        psError(PS_ERR_PROGRAMMING, true, "number of PSF_FWHM values found %ld does not match expected number: %ld\n",
+            fwhmValues->n, max_entries);
+        psAssert(0, "fixme");
+    }
+
+    if (!readout->analysis) {
+        readout->analysis = psMetadataAlloc();
+    }
+
+    psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
+    psFree(fwhmValues);
+
+    return true;
+}

branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourcePhotometry.c

@@ -569 +569 @@
         }
     }
-
     return (true);
 }