/** @file  pmSourceIO.c
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2008-02-03 22:04:47 $
 *
 *  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 "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

// this output format is valid for psphot analysis of an image, and does not include calibrated
// values derived in the DVO database.
// XXX how do I generate the source tables which I need to send to PSPS?

bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader,
                               psMetadata *tableHeader, char *extname, char *xsrcname)
{
    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;

    // if we request XSRC output, add the XSRC name to this header
    if (xsrcname) {
      psMetadataAddStr (tableHeader, PS_LIST_TAIL, "XSRCNAME", PS_META_REPLACE, "name of XSRC table extension", xsrcname);
    }

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

    table = psArrayAllocEmpty (sources->n);

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

        // no difference between PSF and non-PSF model
        pmModel *model = pmSourceGetModel (NULL, source);

        if (model != NULL) {
            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);
        } else {
            xPos = source->peak->xf;
            yPos = source->peak->yf;
            xErr = source->peak->dx;
            yErr = source->peak->dy;
            axes.major = 0.0;
            axes.minor = 0.0;
            axes.theta = 0.0;
        }

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

        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_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);
        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             PS_MIN (99.0, source->psfMag));
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        PS_MIN (99.0, source->errMag));
        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           PS_MIN (99.0, peakMag));
        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
        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 to top 9 pixels",          source->psfChisq);
        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);

        // XXX these should be major and minor, not 'x' and 'y'
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_WIDTH_X",      PS_DATA_F32, "PSF width in x coordinate",                  axes.major);
        psMetadataAdd (row, PS_LIST_TAIL, "PSF_WIDTH_Y",      PS_DATA_F32, "PSF width in y coordinate",                  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);

        // 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, "FLAGS",            PS_DATA_U16, "psphot analysis flags",                      source->mode);

        // XXX these calibrated values are not supplied by psphot analysis
        // psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "RA from PSF fit",                         RA);
        // psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "DEC from PSF fit",                        DEC);
        // psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF_SIG",       PS_DATA_F32, "Sigma of PSF fit RA",                     dRA);
        // psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF_SIG",      PS_DATA_F32, "Sigma of PSF fit DEC",                    dDEC);
        // psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "Calibrated magnitude",                    calMag);
        // psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_sIG",  PS_DATA_F32, "Sigma of calibrated magnitude",           calMagErr);

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

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

    if (xsrcname) {
      pmSourcesWrite_PS1_DEV_1_XSRC (fits, sources, xsrcname);
    }

    return true;
}

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

    // XXX need to look up the XSRCNAME entries

    // validate a single row of the table (must match SMP)

    // XXX test return values

    // XXX we have a memory problem, which is illustrated here: if I allocate the sources array
    // (line 1) before freeing the table, the data is not really freed (but not a leak?)  if I
    // allocate the sources array *after* freeing the table, the data is actually freed
    // psArray *fooSources = psArrayAllocEmpty (10000);
    // psFree (table);
    // return (fooSources);


    // 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.  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 entries
    // XXX need to chooose PSF vs EXT, based on type?
    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;

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

        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_WIDTH_X");
        axes.minor        = psMetadataLookupF32 (&status, row, "PSF_WIDTH_Y");
        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");

        pmPSF_AxesToModel (PAR, axes);

        float lflux = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
        float flux = (isfinite(lflux)) ? pow(10.0, -0.4*lflux) : NAN;
        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], flux, PM_PEAK_LONE);
        source->peak->flux = flux;

        source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");

	// note that some older versions used PSF_PROBABILITY: this was not well defined.
        source->psfChisq  = psMetadataLookupF32 (&status, row, "PSF_CHISQ");
        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");

        source->mode      = psMetadataLookupU16 (&status, row, "FLAGS");
        assert (status);

        // XXX other values saved but not loaded?
        // psMetadataLookupS64 (&status, row, "IPP_IDET");
        // psMetadataLookupF32 (&status, row, "N_FRAMES");

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

    return sources;
}

bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname)
{

    psArray *table;
    psMetadata *row;
    int i;
    psF32 *PAR, *dPAR;
    psEllipseAxes axes;
    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);

    // we write out all sources, regardless of quality.  the source flags tell us the state
    for (i = 0; i < sources->n; i++) {
	// source->seq is set during the PSF output step, which must be called
	// XXX test here if seq is set?

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

        // choose the convolved EXT model, if available, otherwise the simple one
        pmModel *model = source->modelConv;
        if (model == NULL) {
	    model = source->modelEXT;
	}
        if (model == NULL) continue;

	// XXX do we need to do this?
	if (source->extpars == NULL) 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);

        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);
        psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG",     PS_DATA_F32, "EXT fit instrumental magnitude",             source->extMag);
        psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->errMag);

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

	// XXX at the moment, this will be a fixed size table, but perhaps have multiple output formats depending on what is measured?

	// other values that I need to report:
	// R, Mag Petrosian, .....
	if (source->extpars) {

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

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

	    // Isophot measurements
	    // XXX insert header data: isophotal level
	    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);
	    }

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

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

    if (table->n == 0) {
        psFitsWriteBlank (fits, outhead, extname);
        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(table);
        return false;
    }
    psFree (table);

    return true;
}
