/** @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 "pmDetrendDB.h"

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

#include "pmSpan.h"
#include "pmFootprint.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmGrowthCurve.h"
#include "pmResiduals.h"
#include "pmTrend2D.h"
#include "pmPSF.h"
#include "pmModel.h"
#include "pmSource.h"
#include "pmModelClass.h"
#include "pmSourceIO.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 in 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;
    psF32 *PAR, *dPAR;
    psEllipseAxes axes;
    psF32 xPos, yPos;
    psF32 xErr, yErr;
    psF32 errMag, chisq, apRadius;
    psS32 nPix, nDOF;

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

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

    // if the sequence is defined, write these in seq order; otherwise
    // write them in S/N order:
    if (sources->n > 0) {
        pmSource *source = (pmSource *) sources->data[0];
        if (source->seq == -1) {
	    // let's write these out in S/N order
	    sources = psArraySort (sources, pmSourceSortBySN);
        } else {
	    sources = psArraySort (sources, pmSourceSortBySeq);
        }
    }

    table = psArrayAllocEmpty (sources->n);

    // we use this just to define the output vectors (which must be present for all objects)
    bool status = false;
    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!");

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

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

        // no difference between PSF and non-PSF model
        pmModel *model = source->modelPSF;

        if (model != NULL) {
            PAR = model->params->data.F32;
            dPAR = model->dparams->data.F32;
            xPos = PAR[PM_PAR_XPOS];
            yPos = PAR[PM_PAR_YPOS];
            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
		xErr = dPAR[PM_PAR_XPOS];
		yErr = dPAR[PM_PAR_YPOS];
            } else {
		// in linear-fit mode, there is no error on the centroid
		xErr = source->peak->dx;
		yErr = source->peak->dy;
            }
            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
                axes = pmPSF_ModelToAxes (PAR, 20.0);
            } else {
                axes.major = NAN;
                axes.minor = NAN;
                axes.theta = NAN;
            }
            chisq = model->chisq;
            nDOF = model->nDOF;
            nPix = model->nPix;
            apRadius = source->apRadius;
            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
        } else {
            xPos = source->peak->xf;
            yPos = source->peak->yf;
            xErr = source->peak->dx;
            yErr = source->peak->dy;
            axes.major = NAN;
            axes.minor = NAN;
            axes.theta = NAN;
            chisq = NAN;
            nDOF = 0;
            nPix = 0;
            apRadius = NAN;
            errMag = NAN;
        }

        float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
        psS16 nImageOverlap = 1;

        psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
        float posAngle = 0.0;
        float pltScale = 0.0;
        pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);

        row = psMetadataAlloc ();
        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr); // XXX this is only measured for non-linear fits
        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr); // XXX this is only measured for non-linear fits
        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
        psMetadataAdd (row, PS_LIST_TAIL, "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_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_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);

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

        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeight);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);

        // distinguish moments measure from window vs S/N > XX ??
        float mxx = source->moments ? source->moments->Mxx : NAN;
        float mxy = source->moments ? source->moments->Mxy : NAN;
        float myy = source->moments ? source->moments->Myy : NAN;
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);

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

	psVector *radFlux    = psVectorAlloc(radMax->n, PS_TYPE_F32);
	psVector *radFluxErr = psVectorAlloc(radMax->n, PS_TYPE_F32);
	psVector *radFill    = psVectorAlloc(radMax->n, PS_TYPE_F32);
	psVectorInit (radFlux,    NAN);
	psVectorInit (radFluxErr, NAN);
	psVectorInit (radFill,    NAN);
	if (!source->radial) goto empty_annuli;
	if (!source->radial->flux) goto empty_annuli;
	if (!source->radial->fill) goto empty_annuli;
	psAssert (source->radial->flux->n <= radFlux->n, "inconsistent vector lengths");
	psAssert (source->radial->fill->n <= radFlux->n, "inconsistent vector lengths");

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

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

        // XXX not sure how to get this : need to load Nimages with weight?
        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
        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->errMag    = 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;

        pmPSF_AxesToModel (PAR, axes);

        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->flux = peakFlux;
        source->peak->dx   = dPAR[PM_PAR_XPOS];
        source->peak->dy   = dPAR[PM_PAR_YPOS];
        source->peak->SN   = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N

        source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
        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->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
        source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
        source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY");

        source->mode = psMetadataLookupU32 (&status, row, "FLAGS");
        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;

    // 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, pmSourceSortBySN);

    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");
    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");

    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");

    // int nRadialBins = 0;
    // if (doAnnuli) {
    // 	// get the max count of radial bins
    // 	for (int i = 0; i < sources->n; i++) {
    // 	    pmSource *source = sources->data[i];
    // 	    if (!source->extpars) continue;
    // 	    if (!source->extpars->radProfile ) continue;
    //         if (!source->extpars->radProfile->binSB) continue;
    // 	    nRadialBins = PS_MAX(nRadialBins, source->extpars->radProfile->binSB->n);
    // 	}
    // }

    // we write out all sources, regardless of quality.  the source flags tell us the state
    for (int i = 0; i < sources->n; i++) {
        // skip source if it is not a ext sourc
        // XXX we have two places that extended source parameters are measured:
        // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
        // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
        // should we require both?

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

        // 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);
            } 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); 
            }
        }

# if (0)
        // Kron measurements
        if (doKron) {
            pmSourceKronValues *kron = source->extpars->kron;
            if (kron) {
                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       kron->mag);
                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
            } else {
                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
            }
        }

        // Isophot measurements
        // XXX insert header data: isophotal level
        if (doIsophotal) {
            pmSourceIsophotalValues *isophot = source->extpars->isophot;
            if (isophot) {
                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
            } else {
                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
            }
        }
# endif

        // 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;
}

// XXX this layout is still the same as PS1_DEV_1
bool pmSourcesWrite_CMF_PS1_SV1_XFIT (psFits *fits, pmReadout *readout, psArray *sources, 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, pmSourceSortBySN);

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

        // 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);

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

            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) {
                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
                } else {
                    psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
                }
            }

            // XXX other parameters which may be set.
            // XXX flag / value to define the model
            // XXX write out the model type, fit status flags

            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;
}
