/** @file ppSubReadoutPhotometry.c
 *
 *  @brief
 *
 *  @ingroup ppSub
 *
 *  @author IfA
 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-18 00:31:20 $
 *  Copyright 2009 Institute for Astronomy, University of Hawaii
 */

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

#include <stdio.h>
#include <pslib.h>
#include <psmodules.h>
#include <psphot.h>

#include "ppSub.h"

bool ppSubReadoutPhotometry(const char *name, ppSubData *data)
{
    psAssert(data, "Require processing data");
    pmConfig *config = data->config;    // Configuration
    psAssert(config, "Require configuration");

    if (!data->photometry) {
        return true;
    }

    // Look up recipe values
    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim
    psAssert(recipe, "We checked this earlier, so it should be here.");

    if (!psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
        pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
        pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
        return true;
    }

    // The PSF (measured in ppSubMakePSF) is stored on the chip->analysis of PSPHOT.INPUT
    // In order to use an incoming PSF, it must be stored on the chip->analysis of PSPHOT.PSF.LOAD
    pmFPAview *view = ppSubViewReadout(); // View to readout
    pmChip *psfInputChip = pmFPAfileThisChip(config->files, view, "PSPHOT.INPUT"); // Chip with PSF
    psAssert (psfInputChip, "should have been generated for ppSubMakePSF");
    pmChip *psfLoadChip = pmFPAfileThisChip(config->files, view, "PSPHOT.PSF.LOAD"); // Chip to have PSF
    psAssert (psfLoadChip, "PSPHOT.PSF.LOAD should have been defined in ppSubCamera");
    pmPSF *psf = psMetadataLookupPtr(NULL, psfInputChip->analysis, "PSPHOT.PSF"); // PSF for photometry
    if (!psf) {
        psErrorStackPrint(stderr, "No PSF available");
        psWarning("No PSF available --- suspect bad data quality.");
        ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
        return true;
    }
    psMetadataAddPtr(psfLoadChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE,
                     "PSF from ppSubMakePSF", psf);

    bool mdok = false;

    // psphotReadoutMinimal performs the photometry analysis on PSPHOT.INPUT; we need to move
    // around the pointers so PSPHOT.INPUT corresponds to the output image; previously, it was
    // equivalent to the minuend image.
    pmReadout *inRO = pmFPAfileThisReadout(config->files, view, name); // Readout with image and sources

    pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file
    pmFPAfile *inFile = psMetadataLookupPtr(&mdok, config->files, name); // Input file
    if (!pmFPACopy(photFile->fpa, inFile->fpa)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to copy FPA for photometry");
        psFree(view);
        return false;
    }
    pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer
    if (psMetadataLookup(photRO->analysis, "PSPHOT.SOURCES")) {
        psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
    }

    psMetadataAddPtr(photRO->parent->parent->analysis, PS_LIST_TAIL, "PSPHOT.PSF",
                     PS_META_REPLACE | PS_DATA_UNKNOWN, "Point-spread function", data->psf);

    psFree(photRO->analysis);
    photRO->analysis = psMetadataAlloc();

    if (!psphotReadoutMinimal(config, view)) {
        // This is likely a data quality issue
        // XXX Split into multiple cases using error codes?
        psErrorStackPrint(stderr, "Unable to perform photometry on image");
        psWarning("Unable to perform photometry on image --- suspect bad data quality.");
        ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
    }

    if (data->stats) {
        psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
        bool mdok;
        int numSources = psMetadataLookupS32(&mdok, data->stats, "NUM_SOURCES"); // Number of sources
        numSources += sources ? sources->n : 0;
        psMetadataAddS32(data->stats, PS_LIST_TAIL, "NUM_SOURCES", PS_META_REPLACE,
                         "Total number of sources detected", numSources);
        float newTime = psTimerClear("PPSUB_PHOT"); // Time for photometry
        float oldTime = psMetadataLookupF32(&mdok, data->stats, "TIME_PHOT"); // Previous time for photometry
        psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_PHOT", PS_META_REPLACE, "Time to do photometry",
                         isfinite(oldTime) ? oldTime + newTime : newTime);
    }

    if (!data->quality) {
        if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.SOURCES")) {
            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to copy PSPHOT.SOURCES");
            return false;
        }
        if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.HEADER")) {
            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to copy PSPHOT.HEADER");
            return false;
        }
    }

    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");

    return true;
}





#ifdef TESTING
    // Record data about sources: not everything gets into the output CMF files
    {
        pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
        psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
        FILE *sourceFile = fopen("sources.dat", "w"); // File for sources
        fprintf(sourceFile,
                "# x y mag mag_err psf_chisq cr_nsigma ext_nsigma psf_qf flags m_x m_y m_xx m_xy m_yy\n");
        for (int i = 0; i < sources->n; i++) {
            pmSource *source = sources->data[i];
            if (!source) {
                continue;
            }

            float x, y;             // Position of source
            float chi2;             // chi^2 for source
            if (source->modelPSF) {
                x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
                y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
                chi2 = source->modelPSF->chisq;
            } else if (source->peak) {
                x = source->peak->xf;
                y = source->peak->yf;
                chi2 = NAN;
            } else {
                psWarning("No position available for source.");
                continue;
            }

            float xMoment = NAN, yMoment = NAN, xxMoment = NAN, xyMoment = NAN, yyMoment = NAN;
            if (source->moments) {
                xMoment = source->moments->Mx;
                yMoment = source->moments->My;
                xxMoment = source->moments->Mxx;
                xyMoment = source->moments->Mxy;
                yyMoment = source->moments->Myy;
            }

            fprintf(sourceFile, "%f %f %f %f %f %f %f %f %d %f %f %f %f %f\n",
                    x, y, source->psfMag, source->errMag, chi2, source->crNsigma, source->extNsigma,
                    source->pixWeight, source->mode, xMoment, yMoment, xxMoment, xyMoment, yyMoment);
        }
        fclose(sourceFile);
    }
#endif
