/** @file  pmSourceIO.c
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-18 02:44:19 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <pslib.h>

#include "pmConfig.h"
#include "pmErrorCodes.h"
#include "pmDetrendDB.h"

#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPALevel.h"
#include "pmFPAview.h"
#include "pmFPAfile.h"

#include "pmTrend2D.h"
#include "pmResiduals.h"
#include "pmGrowthCurve.h"
#include "pmSpan.h"
#include "pmFootprintSpans.h"
#include "pmFootprint.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmModelFuncs.h"
#include "pmModel.h"
#include "pmModelUtils.h"
#include "pmModelClass.h"
#include "pmSourceMasks.h"
#include "pmSourceExtendedPars.h"
#include "pmSourceDiffStats.h"
#include "pmSource.h"
#include "pmSourceFitModel.h"
#include "pmPSF.h"
#include "pmPSFtry.h"

#include "pmSourceIO.h"
#include "pmSourceOutputs.h"

// panstarrs-style FITS table output (header + table in 1st extension)
// this format consists of a header derived from the image header
// followed by a zero-size matrix, followed by the table data

// NOTE: this output function is intended for psphotStack analysis: it includes per-psf radial fluxes 
// XXX currently, the 'read' function is NOT consistent with the 'write' function (does not read radial fluxes)

bool pmSourcesWrite_CMF_PS1_SV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
{
    PS_ASSERT_PTR_NON_NULL(fits, false);
    PS_ASSERT_PTR_NON_NULL(sources, false);
    PS_ASSERT_PTR_NON_NULL(extname, false);

    psArray *table;
    psMetadata *row;

    pmChip *chip = readout->parent->parent;

    // if the sequence is defined, write these in seq order; otherwise write them in S/N order.
    // Careful: if we are working with child sources, then we need to sort by the parent info,
    // not our info
    if (sources->n > 0) {
        pmSource *source = sources->data[0];
        if (source->seq == -1) {
	    sources = psArraySort (sources, pmSourceSortByFlux);
        } else {
	    sources = psArraySort (sources, pmSourceSortBySeq);
        }
    }

    table = psArrayAllocEmpty (sources->n);

    float magOffset; 
    float zeroptErr; 
    float fwhmMajor; 
    float fwhmMinor;
    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);

    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    for (int i = 0; i < sources->n; i++) {
	// this is the source associated with this image
        pmSource *thisSource = sources->data[i];

	// this is the "real" version of this source 
	pmSource *source = thisSource->parent ? thisSource->parent : thisSource;

        // If source->seq is -1, source is unique and generated in this analysis.  If
        // source->seq is not -1, source was read from elsewhere or tied to other source (eg
        // from another image): in the latter case, preserve the source ID.  source.seq is used
        // instead of source.id since the latter is a const generated on Alloc, and would thus
        // be wrong for read in sources.
        if (source->seq == -1) {
	    source->seq = i;
        }

	// set the 'best' values for various output fields:
	pmSourceOutputs outputs;
	pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);

	pmSourceOutputsMoments moments;
	pmSourceOutputsSetMoments (&moments, source);

        row = psMetadataAlloc ();
        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->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);
        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);

        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);

        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);

        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);

        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);

        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
        psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);

        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                     source->mode);
        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                     source->mode2);

        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
        psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);

        psArrayAdd (table, 100, row);
        psFree (row);

        // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not
        // subtracted

        // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image
        // edge; 3) any pixels in the 3x3 peak region are masked;
    }

    // XXX why do we make a copy here to be supplemented with the masks?  why not do this in the calling function?
    psMetadata *header = psMetadataCopy(NULL, tableHeader);
    pmSourceMasksHeader(header);

    if (table->n == 0) {
        if (!psFitsWriteBlank(fits, header, extname)) {
            psError(psErrorCodeLast(), false, "Unable to write blank sources file.");
            psFree(table);
            psFree(header);
            return false;
        }
        psFree(table);
        psFree(header);
        return true;
    }

    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    if (!psFitsWriteTable(fits, header, table, extname)) {
        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
        psFree(table);
        psFree(header);
        return false;
    }
    psFree(table);
    psFree(header);

    return true;
}

// read in a readout from the fits file
psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header)
{
    PS_ASSERT_PTR_NON_NULL(fits, false);
    PS_ASSERT_PTR_NON_NULL(header, false);

    bool status;
    psF32 *PAR, *dPAR;
    psEllipseAxes axes;

    // define PSF model type
    // XXX need to carry the extra model parameters
    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");

    char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
    if (PSF_NAME != NULL) {
        modelType = pmModelClassGetType (PSF_NAME);
    }
    assert (modelType > -1);

    // We get the size of the table, and allocate the array of sources first because the table
    // is large and ephemeral --- when the table gets blown away, whatever is allocated after
    // the table is read blocks the free.  In fact, it's better to read the table row by row.
    long numSources = psFitsTableSize(fits); // Number of sources in table
    psArray *sources = psArrayAlloc(numSources); // Array of sources, to return

    // convert the table to the pmSource entriesa
    for (int i = 0; i < numSources; i++) {
        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
        if (!row) {
            psError(psErrorCodeLast(), false, "Unable to read row %d of sources", i);
            psFree(sources);
            return NULL;
        }

        pmSource *source = pmSourceAlloc ();
        pmModel *model = pmModelAlloc (modelType);
        source->modelPSF  = model;
        source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags

        // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models.
        PAR = model->params->data.F32;
        dPAR = model->dparams->data.F32;

        source->seq       = psMetadataLookupU32 (&status, row, "IPP_IDET");
        PAR[PM_PAR_XPOS]  = psMetadataLookupF32 (&status, row, "X_PSF");
        PAR[PM_PAR_YPOS]  = psMetadataLookupF32 (&status, row, "Y_PSF");
        dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG");
        dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG");
        axes.major        = psMetadataLookupF32 (&status, row, "PSF_MAJOR");
        axes.minor        = psMetadataLookupF32 (&status, row, "PSF_MINOR");
        axes.theta        = psMetadataLookupF32 (&status, row, "PSF_THETA");

        PAR[PM_PAR_SKY]   = psMetadataLookupF32 (&status, row, "SKY");
        dPAR[PM_PAR_SKY]  = psMetadataLookupF32 (&status, row, "SKY_SIGMA");
        source->sky       = PAR[PM_PAR_SKY];
        source->skyErr    = dPAR[PM_PAR_SKY];

        // XXX use these to determine PAR[PM_PAR_I0]?
        source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
        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->psfMagErr : NAN;

        pmPSF_AxesToModel (PAR, axes, modelType);

        float peakMag     = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
        float peakFlux    = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN;

        // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
        source->peak->rawFlux = peakFlux;
        source->peak->smoothFlux = peakFlux;
        source->peak->xf   = PAR[PM_PAR_XPOS]; // more accurate position
        source->peak->yf   = PAR[PM_PAR_YPOS]; // more accurate position
        source->peak->dx   = dPAR[PM_PAR_XPOS];
        source->peak->dy   = dPAR[PM_PAR_YPOS];

        source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
        source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
        source->apRadius  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");

        // note that some older versions used PSF_PROBABILITY: this was not well defined.
        model->chisq      = psMetadataLookupF32 (&status, row, "PSF_CHISQ");
        model->nDOF       = psMetadataLookupS32 (&status, row, "PSF_NDOF");
        model->nPix       = psMetadataLookupS32 (&status, row, "PSF_NPIX");

        source->moments = pmMomentsAlloc ();
        source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
        source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf

        source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
        source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
        source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY");

        source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
        source->moments->Mrh         = psMetadataLookupF32 (&status, row, "MOMENTS_RH");
        source->moments->KronFlux    = psMetadataLookupF32 (&status, row, "KRON_FLUX");
        source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR");

        source->moments->KronFinner  = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER");
        source->moments->KronFouter  = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER");

	// XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
	// we are storing enough information so the output will be consistent with the input
        source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
        source->moments->Mxxy = 0.0;
        source->moments->Mxyy = 0.0;
        source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S");

        source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C");
        source->moments->Mxxxy = 0.0;
        source->moments->Mxxyy = 0.0;
        source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
        source->moments->Myyyy = 0.0;

        source->mode = psMetadataLookupU32 (&status, row, "FLAGS");
        source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2");
        assert (status);

        sources->data[i] = source;
        psFree(row);
    }

    return sources;
}

bool pmSourcesWrite_CMF_PS1_SV1_XSRC (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
{

    bool status;
    psArray *table;
    psMetadata *row;
    psF32 *PAR, *dPAR;
    psF32 xPos, yPos;
    psF32 xErr, yErr;
    int nRow = -1;
    char keyword1[80], keyword2[80];

    // create a header to hold the output data
    psMetadata *outhead = psMetadataAlloc ();

    // write the links to the image header
    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);

    pmChip *chip = readout->parent->parent;
    pmFPA  *fpa  = chip->parent;

    // zero point corrections
    bool status1 = false;
    bool status2 = false;
    float magOffset = 0.0;
    float exptime   = psMetadataLookupF32(&status1, fpa->concepts, "FPA.EXPOSURE");
    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
    if (!isfinite(zeropt)) {
        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
    }
    if (status1 && status2 && (exptime > 0.0)) {
        magOffset = zeropt + 2.5*log10(exptime);
    }

    // let's write these out in S/N order
    sources = psArraySort (sources, pmSourceSortByFlux);

    table = psArrayAllocEmpty (sources->n);

    // which extended source analyses should we perform?
    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");

    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    psAssert (radMin->n == radMax->n, "inconsistent annular bins");

    // write the radial profile apertures to header
    for (int i = 0; i < radMax->n; i++) {
      sprintf (keyword1, "RMIN_%02d", i);
      sprintf (keyword2, "RMAX_%02d", i);
      psMetadataAddF32 (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]);
    }

    // we write out all sources, regardless of quality.  the source flags tell us the state
    for (int i = 0; i < sources->n; i++) {
	// this is the source associated with this image
        pmSource *thisSource = sources->data[i];

	// this is the "real" version of this source 
	pmSource *source = thisSource->parent ? thisSource->parent : thisSource;

        // skip sources without measurements
        if (source->extpars == NULL) continue;

        // we require a PSF model fit (ignore the real crud)
        pmModel *model = source->modelPSF;
        if (model == NULL) continue;

        // XXX I need to split the extended models from the extended aperture measurements
        PAR = model->params->data.F32;
        dPAR = model->dparams->data.F32;
        xPos = PAR[PM_PAR_XPOS];
        yPos = PAR[PM_PAR_YPOS];
        xErr = dPAR[PM_PAR_XPOS];
        yErr = dPAR[PM_PAR_YPOS];

        row = psMetadataAlloc ();

        // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
        psMetadataAdd (row, PS_LIST_TAIL, "X_EXT",            PS_DATA_F32, "EXT model x coordinate",                     xPos);
        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT",            PS_DATA_F32, "EXT model y coordinate",                     yPos);
        psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG",        PS_DATA_F32, "Sigma in EXT x coordinate",                  xErr);
        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);

	float AxialRatio = NAN;
	float AxialTheta = NAN;
	pmSourceExtendedPars *extpars = source->extpars;
	if (extpars) {
	    AxialRatio = extpars->axes.minor / extpars->axes.major;
	    AxialTheta = extpars->axes.theta;
	}
        psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO",       PS_DATA_F32, "Axial Ratio of radial profile",              AxialRatio);
        psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",            AxialTheta);

        // Petrosian measurements
        // XXX insert header data: petrosian ref radius, flux ratio
	// XXX check flags to see if Pet was measured
        if (doPetrosian) {
	    pmSourceExtendedPars *extpars = source->extpars;
            if (extpars) {
		// XXX note that this mag is either calibrated or instrumental depending on existence of zero point 
		float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
		float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        	 PS_DATA_F32, "Petrosian Magnitude", mag);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    	 PS_DATA_F32, "Petrosian Magnitude Error", magErr);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     	 PS_DATA_F32, "Petrosian Radius (pix)", extpars->petrosianRadius);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", 	 PS_DATA_F32, "Petrosian Radius Error (pix)", extpars->petrosianRadiusErr);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", extpars->petrosianR50);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)", extpars->petrosianR50Err);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL",          PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill);
            } else {
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        	 PS_DATA_F32, "Petrosian Magnitude",       NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    	 PS_DATA_F32, "Petrosian Magnitude Error", NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     	 PS_DATA_F32, "Petrosian Radius",          NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", 	 PS_DATA_F32, "Petrosian Radius Error",    NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50",     PS_DATA_F32, "Petrosian R50 (pix)", NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_50_ERR", PS_DATA_F32, "Petrosian R50 Error (pix)",NAN); 
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN); 
                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL",          PS_DATA_F32, "Petrosian Fill Factor", NAN);
            }
        }

        // Flux Annuli (if we have extended source measurements, we have these.  only optionally save them)
        if (doAnnuli) {
	    psVector *radSB   = psVectorAlloc(radMin->n, PS_TYPE_F32);
	    psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32);
	    psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32);
	    psVectorInit (radSB, NAN);
	    psVectorInit (radFlux, NAN);
	    psVectorInit (radFill, NAN);
	    if (!source->extpars) goto empty_annuli;
	    if (!source->extpars->radProfile) goto empty_annuli;
	    if (!source->extpars->radProfile->binSB) goto empty_annuli;
	    psAssert (source->extpars->radProfile->binSum, "programming error");
	    psAssert (source->extpars->radProfile->binFill, "programming error");
	    psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths");
	    psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths");
	    psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths");

	    // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
	    for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) {
		radSB->data.F32[j]   = source->extpars->radProfile->binSB->data.F32[j];
		radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j];
		radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j];
	    }

	empty_annuli:
	    psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB);
	    psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux);
	    psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill);
	    psFree (radSB);
	    psFree (radFlux);
	    psFree (radFill);
	}
	if (nRow < 0) {
	    nRow = row->list->n;
	} else {
	    psAssert (nRow == row->list->n, "inconsistent row lengths");
	}
	psArrayAdd (table, 100, row);
	psFree (row);
    }
    
    if (table->n == 0) {
	if (!psFitsWriteBlank (fits, outhead, extname)) {
	    psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
	    psFree(outhead);
	    psFree(table);
	    return false;
	}
	psFree (outhead);
	psFree (table);
	return true;
    }
    
    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    if (!psFitsWriteTable (fits, outhead, table, extname)) {
	psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
	psFree (outhead);
    psFree(table);
    return false;
    }
    psFree (outhead);
    psFree (table);
    
    return true;
}

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)
{

    psArray *table;
    psMetadata *row;
    psF32 *PAR, *dPAR;
    psEllipseAxes axes;
    psF32 xPos, yPos;
    psF32 xErr, yErr;
    char name[64];

    // create a header to hold the output data
    psMetadata *outhead = psMetadataAlloc ();

    // write the links to the image header
    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);

    // let's write these out in S/N order
    sources = psArraySort (sources, pmSourceSortByFlux);

    // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
    int nParamMax = 0;
    for (int i = 0; i < sources->n; i++) {
	// this is the source associated with this image
        pmSource *thisSource = sources->data[i];

	// this is the "real" version of this source 
	pmSource *source = thisSource->parent ? thisSource->parent : thisSource;

        if (source->modelFits == NULL) continue;
        for (int j = 0; j < source->modelFits->n; j++) {
            pmModel *model = source->modelFits->data[j];
            assert (model);
            nParamMax = PS_MAX (nParamMax, model->params->n);
        }
    }

    table = psArrayAllocEmpty (sources->n);

    // we write out all sources, regardless of quality.  the source flags tell us the state
    for (int i = 0; i < sources->n; i++) {

        pmSource *thisSource = sources->data[i];

	// this is the "real" version of this source 
	pmSource *source = thisSource->parent ? thisSource->parent : thisSource;

        // XXX if no model fits are saved, write out modelEXT?
        if (source->modelFits == NULL) continue;

        // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
        for (int j = 0; j < source->modelFits->n; j++) {

            // choose the convolved EXT model, if available, otherwise the simple one
            pmModel *model = source->modelFits->data[j];
            assert (model);

	    // skip models which were not actually fitted
	    if (model->flags & PM_MODEL_STATUS_BADARGS) continue;

            PAR = model->params->data.F32;
            dPAR = model->dparams->data.F32;
            xPos = PAR[PM_PAR_XPOS];
            yPos = PAR[PM_PAR_YPOS];
            xErr = dPAR[PM_PAR_XPOS];
            yErr = dPAR[PM_PAR_YPOS];

            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);

            row = psMetadataAlloc ();

            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);

            psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
            psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));

            // XXX these should be major and minor, not 'x' and 'y'
            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width in x coordinate",                  axes.major);
            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width in y coordinate",                  axes.minor);
            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                      axes.theta);

            // write out the other generic parameters
            for (int k = 0; k < nParamMax; k++) {
                if (k == PM_PAR_I0) continue;
                if (k == PM_PAR_SKY) continue;
                if (k == PM_PAR_XPOS) continue;
                if (k == PM_PAR_YPOS) continue;
                if (k == PM_PAR_SXX) continue;
                if (k == PM_PAR_SXY) continue;
                if (k == PM_PAR_SYY) continue;

                snprintf (name, 64, "EXT_PAR_%02d", k);

                if (k < model->params->n) {
                    psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]);
                } else {
                    psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", NAN);
                }
            }

	    // optionally, write out the covariance matrix values
	    // XXX do I need to pad this to match the biggest covar matrix?
	    if (model->covar) {
		for (int iy = 0; iy < model->covar->numCols; iy++) {
		    for (int ix = iy; ix < model->covar->numCols; ix++) {
			snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix);
			psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]);

		    }
		}		    
	    }
            psArrayAdd (table, 100, row);
            psFree (row);
        }
    }

    if (table->n == 0) {
        if (!psFitsWriteBlank (fits, outhead, extname)) {
            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
            psFree(outhead);
            psFree(table);
            return false;
        }
        psFree (outhead);
        psFree (table);
        return true;
    }

    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    if (!psFitsWriteTable (fits, outhead, table, extname)) {
        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
        psFree (outhead);
        psFree(table);
        return false;
    }
    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;
}

// **** write out the radial flux values for the sources for a given matched-PSF image
// **** how do we distinguish the matched-PSF images from the non-matched version
bool pmSourcesWrite_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
{
    bool status = false;
    psArray *table;
    psMetadata *row;
    psF32 xPos, yPos;
    char keyword1[80], keyword2[80];

    // perform full non-linear fits / extended source analysis?
    if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
	psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n");
	return true;
    }

    // create a header to hold the output data
    psMetadata *outhead = psMetadataAlloc ();

    // write the links to the image header
    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "radial flux table extension", extname);

    // we use this just to define the output vectors (which must be present for all objects)
    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    psAssert (radMax, "this must have been defined and tested earlier!");
    psAssert (radMax->n, "this must have been defined and tested earlier!");
    psAssert (radMin->n == radMax->n, "inconsistent annular bins");

    // write the radial profile apertures to header
    for (int i = 0; i < radMax->n; i++) {
      sprintf (keyword1, "RMIN_%02d", i);
      sprintf (keyword2, "RMAX_%02d", i);
      psMetadataAddF32 (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]);
    }

    // the FWHM values are available if we measured a psf-matched convolved set
    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");

    // let's write these out in S/N order
    sources = psArraySort (sources, pmSourceSortByFlux);

    table = psArrayAllocEmpty (sources->n);

    // we write out all sources, regardless of quality.  the source flags tell us the state
    for (int i = 0; i < sources->n; i++) {

	// this is the source associated with this image
        pmSource *thisSource = sources->data[i];

	// this is the "real" version of this source 
	pmSource *source = thisSource->parent ? thisSource->parent : thisSource;

        // skip sources without radial aper measurements (or insufficient)
	if (source->radialAper == NULL) continue;

        // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");

	for (int entry = 0; entry < source->radialAper->n; entry++) {

	    // choose the convolved EXT model, if available, otherwise the simple one
	    pmSourceRadialApertures *radialAper = source->radialAper->data[entry];
	    assert (radialAper);

	    if (pmSourcePositionUseMoments(source)) {
		xPos = source->moments->Mx;
		yPos = source->moments->My;
	    } else {
		xPos = source->peak->xf;
		yPos = source->peak->yf;
	    }

	    row = psMetadataAlloc ();

	    // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
	    psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
	    psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
	    psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
	    if (fwhmValues) {
		psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                    fwhmValues->data.F32[entry]);
	    } else {
		psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "image is not FWHM-matched",              NAN);
	    }

	    // XXX if we have raw radial apertures, write them out here
	    psVector *radFlux      = psVectorAlloc(radMax->n, PS_TYPE_F32);
	    psVector *radFluxErr   = psVectorAlloc(radMax->n, PS_TYPE_F32);
	    psVector *radFill      = psVectorAlloc(radMax->n, PS_TYPE_F32);
	    psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32);
	    psVectorInit (radFlux,    NAN);
	    psVectorInit (radFluxErr, NAN);
	    psVectorInit (radFill,    NAN);
	    if (!radialAper->flux) goto write_annuli;
	    if (!radialAper->fill) goto write_annuli;
	    psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths");
	    psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths");

	    // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
	    for (int j = 0; j < radialAper->flux->n; j++) {
		radFlux->data.F32[j]      = radialAper->flux->data.F32[j];
		radFluxErr->data.F32[j]   = radialAper->fluxErr->data.F32[j];
		radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j];
		radFill->data.F32[j]      = radialAper->fill->data.F32[j];
	    }

	write_annuli:
	    psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",     	 PS_DATA_VECTOR, "flux within annuli",       radFlux);
	    psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR", 	 PS_DATA_VECTOR, "flux error in annuli",     radFluxErr);
	    psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation",  radFluxStdev);
	    psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",       PS_DATA_VECTOR, "fill factor of annuli",    radFill);
	    psFree (radFlux);
	    psFree (radFluxErr);
	    psFree (radFluxStdev);
	    psFree (radFill);

	    psArrayAdd (table, 100, row);
	    psFree (row);
	}
    }

    if (table->n == 0) {
        if (!psFitsWriteBlank (fits, outhead, extname)) {
            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
            psFree(outhead);
            psFree(table);
            return false;
        }
        psFree (outhead);
        psFree (table);
        return true;
    }

    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    if (!psFitsWriteTable (fits, outhead, table, extname)) {
        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
        psFree (outhead);
        psFree(table);
        return false;
    }
    psFree (outhead);
    psFree (table);
    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;
}
