/** @file  pmSourceIO.c
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.3.2.2 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-19 17:59:50 $
 *
 *  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

bool pmSourcesWrite_CMF_PS1_V1 (psFits *fits, pmReadout *readout, psArray *sources, 
				psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
{
    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;
    int i;
    psF32 *PAR, *dPAR;
    psEllipseAxes axes;
    psF32 xPos, yPos;
    psF32 xErr, yErr;
    psF32 errMag, chisq, apRadius;
    psS32 nPix, nDOF;

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

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

    // 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 write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    for (i = 0; i < sources->n; i++) {
        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 = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius?
	    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, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
        psMetadataAdd (row, PS_LIST_TAIL, "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, "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, "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 ??
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      source->moments->Mxx);
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      source->moments->Mxy);
        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      source->moments->Myy);

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

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

    if (table->n == 0) {
        psFitsWriteBlank (fits, tableHeader, extname);
        psFree (table);
        return true;
    }

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

    return true;
}

// read in a readout from the fits file
psArray *pmSourcesRead_CMF_PS1_V1 (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

        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->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");

	// 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");
        model->radiusFit  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");

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

// XXX this layout is still the same as PS1_DEV_1 
bool pmSourcesWrite_CMF_PS1_V1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
{

    bool status;
    psArray *table;
    psMetadata *row;
    psF32 *PAR, *dPAR;
    psF32 xPos, yPos;
    psF32 xErr, yErr;

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

    table = psArrayAllocEmpty (sources->n);

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

    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
    psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    assert (radialBinsLower->n == radialBinsUpper->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);

	// Petrosian measurements
	// XXX insert header data: petrosian ref radius, flux ratio
	if (doPetrosian) {
	    pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
	    if (petrosian) {
		psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
		psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
		psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
		psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
	    } 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);
	    }
	} 

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

	// Flux Annuli
	if (doAnnuli) {
	    pmSourceAnnuli *annuli = source->extpars->annuli;
	    if (annuli) {
		psVector *fluxVal = annuli->flux;
		psVector *fluxErr = annuli->fluxErr;
		psVector *fluxVar = annuli->fluxVar;

		for (int j = 0; j < fluxVal->n; j++) {
		    char name[32];
		    sprintf (name, "FLUX_VAL_R_%02d", j);
		    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
		    sprintf (name, "FLUX_ERR_R_%02d", j);
		    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
		    sprintf (name, "FLUX_VAR_R_%02d", j);
		    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
		} 
	    } else {
		for (int j = 0; j < radialBinsLower->n; j++) {
		    char name[32];
		    sprintf (name, "FLUX_VAL_R_%02d", j);
		    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
		    sprintf (name, "FLUX_ERR_R_%02d", j);
		    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
		    sprintf (name, "FLUX_VAR_R_%02d", j);
		    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
		} 
	    }
	}

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

    if (table->n == 0) {
	psFitsWriteBlank (fits, outhead, extname);
	psFree (outhead);
	psFree (table);
	return true;
    }

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

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

bool pmSourceLocalAstrometry (psSphere *ptSky, float *posAngle, float *pltScale, pmChip *chip, float xPos, float yPos) {

    pmFPA *fpa = chip->parent;

    if (!chip->toFPA) goto escape;
    if (!fpa->toTPA) goto escape;
    if (!fpa->toSky) goto escape;

    // generate RA,DEC
    psPlane ptCH, ptFP, ptTP_o, ptTP_x, ptTP_y;

    // calculate the astrometry for the coordinate of interest
    ptCH.x = xPos;
    ptCH.y = yPos;
    psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    psPlaneTransformApply (&ptTP_o, fpa->toTPA, &ptFP);
    psDeproject (ptSky, &ptTP_o, fpa->toSky);

    // calculate the astrometry for the coordinate + 1pix in X
    ptCH.x = xPos + 1.0;
    ptCH.y = yPos;
    psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    psPlaneTransformApply (&ptTP_x, fpa->toTPA, &ptFP);

    // calculate the astrometry for the coordinate + 1pix in Y
    ptCH.x = xPos;
    ptCH.y = yPos + 1.0;
    psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    psPlaneTransformApply (&ptTP_y, fpa->toTPA, &ptFP);
    
    // the resulting Tangent Plane coordinates are in TP pixels; convert to local Tangent Plane
    // degrees

    float dTPx_dCHx = fpa->toSky->Xs * (ptTP_x.x - ptTP_o.x);
    float dTPy_dCHx = fpa->toSky->Ys * (ptTP_x.y - ptTP_o.y);

    float dTPx_dCHy = fpa->toSky->Xs * (ptTP_y.x - ptTP_o.x);
    float dTPy_dCHy = fpa->toSky->Ys * (ptTP_y.y - ptTP_o.y);

    float pltScale_x = hypot(dTPx_dCHx, dTPy_dCHx);
    float pltScale_y = hypot(dTPx_dCHy, dTPy_dCHy);
    *pltScale = 0.5*(pltScale_x + pltScale_y);

    float posAngle_x = atan2 (+dTPy_dCHx, +dTPx_dCHx);
    float posAngle_y = atan2 (-dTPy_dCHy, +dTPx_dCHy);
    *posAngle = 0.5*(posAngle_x + posAngle_y);

    return true;

escape:
    // no astrometry calibration, give up
    ptSky->r = NAN;
    ptSky->d = NAN;
    *posAngle = NAN;
    *pltScale = NAN;
    
    return false;
}
