/** @file  pmPSF_IO.c
 *
 * This file contains functions to read and write PSF models using the psMetadata Config file
 * format.
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2007-05-04 00:53:15 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

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

/*****************************************************************************/
/* INCLUDE FILES                                                             */
/*****************************************************************************/

#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 "pmModel.h"
#include "pmPSF.h"
#include "pmPSF_IO.h"
#include "pmSource.h"
#include "pmModelGroup.h"
#include "pmSourceIO.h"

bool pmFPAviewCheckDataStatusForPSFmodel (const pmFPAview *view, const pmFPAfile *file)
{
    pmFPA *fpa = file->fpa;

    if (view->chip == -1) {
	bool exists = pmFPACheckDataStatusForPSFmodel (fpa);
        return exists;
    }
    if (view->chip >= fpa->chips->n) {
	psError(PS_ERR_IO, true, "Requested chip == %d >= fpa->chips->n == %ld", view->chip, fpa->chips->n);
        return false;
    }
    pmChip *chip = fpa->chips->data[view->chip];

    if (view->cell == -1) {
        bool exists = pmChipCheckDataStatusForPSFmodel (chip);
        return exists;
    }
    if (view->cell >= chip->cells->n) {
	psError(PS_ERR_IO, true, "Requested cell == %d >= chip->cells->n == %ld", view->cell, chip->cells->n);
        return false;
    }
    pmCell *cell = chip->cells->data[view->cell];

    if (view->readout == -1) {
        bool exists = pmCellCheckDataStatusForPSFmodel (cell);
        return exists;
    }

    if (view->readout >= cell->readouts->n) {
	psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouds->n == %ld", view->readout, cell->readouts->n);
        return false;
    }
    pmReadout *readout = cell->readouts->data[view->readout];

    bool exists = pmReadoutCheckDataStatusForPSFmodel (readout);
    return exists;
}

bool pmFPACheckDataStatusForPSFmodel (const pmFPA *fpa) {

    for (int i = 0; i < fpa->chips->n; i++) {
	pmChip *chip = fpa->chips->data[i];
	if (!chip) continue;
	if (pmChipCheckDataStatusForPSFmodel (chip)) return true;
    }
    return false;
}

bool pmChipCheckDataStatusForPSFmodel (const pmChip *chip) {

    for (int i = 0; i < chip->cells->n; i++) {
	pmCell *cell = chip->cells->data[i];
	if (!cell) continue;
	if (pmCellCheckDataStatusForPSFmodel (cell)) return true;
    }
    return false;
}

bool pmCellCheckDataStatusForPSFmodel (const pmCell *cell) {

    for (int i = 0; i < cell->readouts->n; i++) {
	pmReadout *readout = cell->readouts->data[i];
	if (!readout) continue;
	if (pmReadoutCheckDataStatusForPSFmodel (readout)) return true;
    }
    return false;
}

bool pmReadoutCheckDataStatusForPSFmodel (const pmReadout *readout) {

    bool status;

    // select the psf of interest
    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    if (!psf) return false;
    return true;
}

bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{

    pmFPA *fpa = file->fpa;

    if (view->chip == -1) {
        if (!pmFPAWritePSFmodel (fpa, view, file, config)) {
            psError(PS_ERR_IO, false, "Failed to write PSF for fpa");
            return false;
        }
        return true;
    }

    if (view->chip >= fpa->chips->n) {
        return false;
    }
    pmChip *chip = fpa->chips->data[view->chip];

    if (view->cell == -1) {
        if (!pmChipWritePSFmodel (chip, view, file, config)) {
            psError(PS_ERR_IO, false, "Failed to write PSF for chip");
            return false;
        }
        return true;
    }

    if (view->cell >= chip->cells->n) {
        return false;
    }
    pmCell *cell = chip->cells->data[view->cell];

    if (view->readout == -1) {
        if (!pmCellWritePSFmodel (cell, view, file, config)) {
            psError(PS_ERR_IO, false, "Failed to write PSF for cell");
            return false;
        }
        return true;
    }

    if (view->readout >= cell->readouts->n) {
        return false;
    }
    pmReadout *readout = cell->readouts->data[view->readout];

    if (!pmReadoutWritePSFmodel (readout, view, file, config)) {
	psError(PS_ERR_IO, false, "Failed to write PSF for readout");
	return false;
    }
    return true;
}

// read in all chip-level PSFmodel files for this FPA
bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{

    for (int i = 0; i < fpa->chips->n; i++) {

        pmChip *chip = fpa->chips->data[i];
        if (!pmChipWritePSFmodel (chip, view, file, config)) {
            psError(PS_ERR_IO, false, "Failed to write PSF for %dth chip", i);
            return false;
        }
    }
    return true;
}

// read in all cell-level PSFmodel files for this chip
bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{

    for (int i = 0; i < chip->cells->n; i++) {

        pmCell *cell = chip->cells->data[i];
        if (!pmCellWritePSFmodel (cell, view, file, config)) {
            psError(PS_ERR_IO, false, "Failed to write PSF for %dth cell", i);
            return false;
        }
    }
    return true;
}

// read in all readout-level PSFmodel files for this cell
bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{

    for (int i = 0; i < cell->readouts->n; i++) {

        pmReadout *readout = cell->readouts->data[i];
        if (!pmReadoutWritePSFmodel (readout, view, file, config)) {
            psError(PS_ERR_IO, false, "Failed to write PSF for %dth readout", i);
            return false;
        }
    }
    return true;
}

// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, 
// and header + image for the PSF residual images
bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{
    bool status;

    // select the psf of interest
    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    if (!psf) {
	psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout");
	return false;
    }

    // lookup the EXTNAME values used for table data and residual image segments
    char *rule = NULL;

    // Menu of EXTNAME rules
    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
    if (!menu) {
	psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
	return false;
    }
    // EXTNAME for table data
    rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
    if (!rule) {
	psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
	return false;
    }
    char *tableName = pmFPAfileNameFromRule (rule, file, view);

    // write the PSF model parameters in a FITS table

    // we need to write a header for the table, 
    psMetadata *header = psMetadataAlloc();

    char *modelName = pmModelGetType (psf->type);
    psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName);

    psMetadataAddBool (header, PS_LIST_TAIL, "POISSON",  0, "Use Poisson errors in fits?", psf->poissonErrors);

    int nPar = pmModelParameterCount (psf->type)    ;
    psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);

    // save the dimensions of each parameter
    for (int i = 0; i < nPar; i++) {
	char name[9];
        psPolynomial2D *poly = psf->params_NEW->data[i];
        if (poly == NULL) continue;
	snprintf (name, 9, "PAR%02d_NX", i);
        psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX);
	snprintf (name, 9, "PAR%02d_NY", i);
        psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY);
    }

    // other required information describing the PSF
    psMetadataAddF32  (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
    psMetadataAddF32  (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
    psMetadataAddF32  (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
    psMetadataAddS32  (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);

    // XXX can we drop this now?
    psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);

    // build a FITS table of the PSF parameters
    psArray *psfTable = psArrayAllocEmpty (100);
    for (int i = 0; i < nPar; i++) {
        psPolynomial2D *poly = psf->params_NEW->data[i];
        if (poly == NULL) continue; // skip unset parameters (eg, XPOS)
	for (int ix = 0; ix < poly->nX; ix++) {
	    for (int iy = 0; iy < poly->nY; iy++) {

		psMetadata *row = psMetadataAlloc ();
		psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
		psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
		psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
		psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
		psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
		psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);

		psArrayAdd (psfTable, 100, row);
		psFree (row);
	    }
	}
    }
    
    // write an empty FITS segment if we have no PSF information
    if (psfTable->n == 0) {
	// XXX this is probably an error (if we have a PSF, how do we have no data?)
        psFitsWriteBlank (file->fits, header, tableName);
    } else {
	psTrace ("pmFPAfile", 5, "writing psf data %s\n", tableName);
	if (!psFitsWriteTable (file->fits, header, psfTable, tableName)) {
	    psError(PS_ERR_IO, false, "writing psf table data %s\n", tableName);
	    psFree (tableName);
	    psFree (psfTable);
	    psFree (header);
	    return false;
	}
    }
    psFree (tableName);
    psFree (psfTable);
    psFree (header);

    // EXTNAME for residual images
    rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
    if (!rule) {
	psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
	return false;
    }
    char *imageName = pmFPAfileNameFromRule (rule, file, view);

    // write the residual images (3D) 
    header = psMetadataAlloc ();
    if (psf->residuals == NULL) {
	// set some header keywords to make it clear there are no residuals?
        psFitsWriteBlank (file->fits, header, imageName);
	psFree (imageName);
        psFree (header);
        return true;
    }

    psMetadataAddS32 (header, PS_LIST_TAIL, "XBIN",    0, "", psf->residuals->xBin);
    psMetadataAddS32 (header, PS_LIST_TAIL, "YBIN",    0, "", psf->residuals->yBin);
    psMetadataAddS32 (header, PS_LIST_TAIL, "XCENTER", 0, "", psf->residuals->xCenter);
    psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter);

    // write the residuals as three planes of the image
    // this call creates an extension with NAXIS3 = 3
    if (psf->residuals->Rx) {
	// this call creates an extension with NAXIS3 = 3
	psArray *images = psArrayAllocEmpty (3);
	psArrayAdd (images, 1, psf->residuals->Ro);
	psArrayAdd (images, 1, psf->residuals->Rx);
	psArrayAdd (images, 1, psf->residuals->Ry);

	psFitsWriteImageCube (file->fits, header, images, imageName);
	psFree (images);
    } else {
	// this call creates an extension with NAXIS3 = 1
	psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, imageName);
    }
  
    psFree (imageName);
    psFree (header);
    return true;

    // XXX save the growth curve
    // XXX save the ApResid fit

# if (0)
    // build a FITS table of the fit to the Aperture Residuals
    psArray *apresTable = psArrayAllocEmpty (100);
    psPolynomial4D *poly = psf->ApTrend;
    for (int ix = 0; ix < poly->nX; ix++) {
	for (int iy = 0; iy < poly->nY; iy++) {

	    row = psMetadataAlloc ();
	    psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
	    psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
	    psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
	    psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
	    psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);

	    psArrayAdd (psfTable, 100, row);
	    psFree (row);
	}
    }
# endif
}

// XXX add in error handling
bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{

    pmFPA *fpa = file->fpa;

    if (view->chip == -1) {
        pmFPAReadPSFmodel (fpa, view, file, config);
        return true;
    }

    if (view->chip >= fpa->chips->n) {
        return false;
    }
    pmChip *chip = fpa->chips->data[view->chip];

    if (view->cell == -1) {
        pmChipReadPSFmodel (chip, view, file, config);
        return true;
    }

    if (view->cell >= chip->cells->n) {
        return false;
    }
    pmCell *cell = chip->cells->data[view->cell];

    if (view->readout == -1) {
        pmCellReadPSFmodel (cell, view, file, config);
        return true;
    }

    if (view->readout >= cell->readouts->n) {
        return false;
    }
    pmReadout *readout = cell->readouts->data[view->readout];

    pmReadoutReadPSFmodel (readout, view, file, config);
    return true;
}

// read in all chip-level PSFmodel files for this FPA
bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{

    for (int i = 0; i < fpa->chips->n; i++) {

        pmChip *chip = fpa->chips->data[i];
        pmChipReadPSFmodel (chip, view, file, config);
    }
    return true;
}

// read in all cell-level PSFmodel files for this chip
bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{

    for (int i = 0; i < chip->cells->n; i++) {

        pmCell *cell = chip->cells->data[i];
        pmCellReadPSFmodel (cell, view, file, config);
    }
    return true;
}

// read in all readout-level PSFmodel files for this cell
bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{

    for (int i = 0; i < cell->readouts->n; i++) {

        pmReadout *readout = cell->readouts->data[i];
        pmReadoutReadPSFmodel (readout, view, file, config);
    }
    return true;
}

// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, 
// and header + image for the PSF residual images
bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{
    bool status; 
    char *rule = NULL;
    psMetadata *header = NULL;

    // Menu of EXTNAME rules
    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
    if (!menu) {
	psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
	return false;
    }
    // EXTNAME for table data
    rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
    if (!rule) {
	psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
	return false;
    }
    char *tableName = pmFPAfileNameFromRule (rule, file, view);
    // EXTNAME for residual images
    rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
    if (!rule) {
	psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
	return false;
    }
    char *imageName = pmFPAfileNameFromRule (rule, file, view);

    // move fits pointer to table and read header
    // advance to the table data extension
    // since we have read the IMAGE header, the TABLE header should exist
    if (!psFitsMoveExtName (file->fits, tableName)) {
	psAbort("cannot find data extension %s in %s", tableName, file->filename);
    }

    // load the PSF model table header
    header = psFitsReadHeader (NULL, file->fits);
    if (!header) psAbort("cannot read table header");

    // load the PSF model parameters from the FITS table
    char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME");
    pmModelType type = pmModelSetType (modelName);

    bool poissonErrors = psMetadataLookupBool (&status, header, "POISSON");

    // we determine the PSF parameter polynomials from the MD-defined polynomials
    pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL);

    // check the number of expected parameters
    int nPar = psMetadataLookupS32 (&status, header, "PSF_NPAR");
    if (nPar != pmModelParameterCount (psf->type))
        psAbort("mismatch model par count");

    // load the dimensions of each parameter
    for (int i = 0; i < nPar; i++) {
	char name[9];
	snprintf (name, 9, "PAR%02d_NX", i);
        int nXorder = psMetadataLookupS32 (&status, header, name);
	snprintf (name, 9, "PAR%02d_NY", i);
        int nYorder = psMetadataLookupS32 (&status, header, name);
	psf->params_NEW->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder);
    }

    // other required information describing the PSF
    psf->ApResid   = psMetadataLookupF32 (&status, header, "AP_RESID");
    psf->dApResid  = psMetadataLookupF32 (&status, header, "AP_ERROR");
    psf->chisq     = psMetadataLookupF32 (&status, header, "CHISQ");
    psf->nPSFstars = psMetadataLookupS32 (&status, header, "NSTARS");

    // XXX can we drop this now?
    psf->skyBias   = psMetadataLookupF32 (&status, header, "SKY_BIAS");

    // read the raw table data
    psArray *table = psFitsReadTable (file->fits);

    // fill in the matching psf->params entries
    for (int i = 0; i > table->n; i++) {
	psMetadata *row = table->data[i];
	int iPar = psMetadataLookupS32 (&status, row, "MODEL_TERM");
	int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
	int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
	// XXX sanity check here

        psPolynomial2D *poly = psf->params_NEW->data[iPar];

	poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
	poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
	poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
    }
    psFree (header);
    psFree (table);

    // move fits pointer to residual image and read header
    // advance to the table data extension
    // since we have read the IMAGE header, the TABLE header should exist
    if (!psFitsMoveExtName (file->fits, imageName)) {
	psAbort("cannot find data extension %s in %s", imageName, file->filename);
    }

    header = psFitsReadHeader (NULL, file->fits);
    int Naxis = psMetadataLookupS32 (&status, header, "NAXIS");
    if (Naxis != 0) {

	int Nx = psMetadataLookupS32 (&status, header, "NAXIS1");
	int Ny = psMetadataLookupS32 (&status, header, "NAXIS2");
	int Nz = psMetadataLookupS32 (&status, header, "NAXIS3");

	int xBin  = psMetadataLookupS32 (&status, header, "XBIN");
	int yBin  = psMetadataLookupS32 (&status, header, "YBIN");

	int xSize = Nx / xBin;
	int ySize = Ny / yBin;

	psf->residuals = pmResidualsAlloc (xSize, ySize, xBin, yBin);

	psf->residuals->xCenter = psMetadataLookupS32 (&status, header, "XCENTER");
	psf->residuals->yCenter = psMetadataLookupS32 (&status, header, "YCENTER");

	psRegion fullImage = {0, 0, 0, 0};
	psFitsReadImageBuffer(psf->residuals->Ro, file->fits, fullImage, 0); // Desired pixels
	if (Nz > 1) {
	    assert (Nz == 3);
	    psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1); // Desired pixels
	    psFitsReadImageBuffer(psf->residuals->Ry, file->fits, fullImage, 2); // Desired pixels
	}
    }

    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
    psFree (psf);

    psFree (tableName);
    psFree (imageName);

    return true;
}

/************ old support functions, deprecate? **************/

psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf)
{

    if (metadata == NULL) {
        metadata = psMetadataAlloc ();
    }

    char *modelName = pmModelGetType (psf->type);
    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);

    int nPar = pmModelParameterCount (psf->type)    ;
    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);

    for (int i = 0; i < nPar; i++) {
        psPolynomial2D *poly = psf->params_NEW->data[i];
        if (poly == NULL)
            continue;
        psPolynomial2DtoMetadata (metadata, poly, "PSF_PAR%02d", i);
    }
    psPolynomial4DtoMetadata (metadata, psf->ApTrend, "APTREND");

    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid);
    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid);
    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);

    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq);
    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars);
    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_POISSON_ERRORS", PS_DATA_BOOL, "Poisson errors for fits", psf->poissonErrors);

    return metadata;
}

pmPSF *pmPSFfromMetadata (psMetadata *metadata)
{

    bool status;
    char keyword[80];

    char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME");
    pmModelType type = pmModelSetType (modelName);

    bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS");
    if (!status)
        poissonErrors = true;

    // we determine the PSF parameter polynomials from the MD-defined polynomials
    pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL);

    int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR");
    if (nPar != pmModelParameterCount (psf->type))
        psAbort("mismatch model par count");

    // un-fitted terms, not in the Metadata, are left NULL
    // XXX add a double-check of the expected number?
    for (int i = 0; i < nPar; i++) {
        sprintf (keyword, "PSF_PAR%02d", i);
        psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword);
        if (!status)
            continue;
        psPolynomial2D *poly = psPolynomial2DfromMetadata (folder);
        psFree (psf->params_NEW->data[i]);
        psf->params_NEW->data[i] = poly;
    }
    sprintf (keyword, "APTREND");
    psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword);
    psPolynomial4D *poly = psPolynomial4DfromMetadata (folder);
    psFree (psf->ApTrend);
    psf->ApTrend = poly;

    psf->ApResid = psMetadataLookupF32 (&status, metadata, "PSF_AP_RESID");
    psf->dApResid = psMetadataLookupF32 (&status, metadata, "PSF_dAP_RESID");
    psf->skyBias = psMetadataLookupF32 (&status, metadata, "PSF_SKY_BIAS");

    psf->chisq = psMetadataLookupF32 (&status, metadata, "PSF_CHISQ");
    psf->nPSFstars = psMetadataLookupS32 (&status, metadata, "PSF_NSTARS");

    psFree (metadata);
    return (psf);
}

// read in all readout-level Objects files for this cell
bool pmReadoutWritePSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{
    bool status;
    char *filename;
    char *realname;

    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");

    switch (file->type) {
      case PM_FPA_FILE_PSF:
        filename = pmFPAfileNameFromRule (file->filerule, file, view);
        bool create = file->mode == PM_FPA_MODE_WRITE ? true : false;
        realname = pmConfigConvertFilename (filename, config, create);

        psMetadata *psfData = pmPSFtoMetadata (NULL, psf);
        psMetadataConfigWrite (psfData, realname);
        psFree (psfData);
        psFree (realname);
        psFree (filename);
        return true;

      default:
        fprintf (stderr, "warning: type mismatch\n");
        break;
    }
    return false;
}

// read in all readout-level Objects files for this cell
bool pmReadoutReadPSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
{

    unsigned int Nfail;
    char *filename;
    char *realname;

    switch (file->type) {
      case PM_FPA_FILE_PSF:
        filename = pmFPAfileNameFromRule (file->filerule, file, view);
        bool create = file->mode == PM_FPA_MODE_WRITE ? true : false;
        realname = pmConfigConvertFilename (filename, config, create);

        psMetadata *psfData = psMetadataConfigRead(NULL, &Nfail, realname, FALSE);
        pmPSF *psf = pmPSFfromMetadata (psfData);
        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf);

        psFree (psf);
        psFree (psfData);
        psFree (realname);
        psFree (filename);

        return true;

      default:
        fprintf (stderr, "warning: type mismatch\n");
        break;
    }
    return false;
}

