# include "psphotInternal.h"

pmReadout *psphotSelectBackground (pmConfig *config,
                                   const pmFPAview *view) {

    bool status;
    pmReadout *background;

    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL");
    if (!file) return NULL;
    if (file->mode == PM_FPA_MODE_INTERNAL) {
        background = file->readout;
    } else {
        background = pmFPAviewThisReadout (view, file->fpa);
    }
    return background;
}

pmReadout *psphotSelectBackgroundStdev (pmConfig *config,
                                        const pmFPAview *view) {

    bool status;
    pmReadout *background;

    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL.STDEV");
    if (!file) return NULL;
    if (file->mode == PM_FPA_MODE_INTERNAL) {
        background = file->readout;
    } else {
        background = pmFPAviewThisReadout (view, file->fpa);
    }
    return background;
}

// dump source stats for psf stars
bool psphotDumpStats (psArray *sources, char *stage) {

    char filename[64];
    snprintf (filename, 64, "psf.%s.dat", stage);
    FILE *f = fopen (filename, "w");
    for (int i = 0; i < sources->n; i++) {
	pmSource *source = sources->data[i];
	if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;

	pmModel *model = source->modelPSF;
	if (!model) continue;

	// int xc = source->peak->x - source->pixels->col0;
	// int yc = source->peak->y - source->pixels->row0;
	// float mcore = source->modelFlux ? source->modelFlux->data.F32[yc][xc] : NAN;
	// float mpeak = model ? model->params->data.F32[PM_PAR_I0] : NAN;
	// bool subtracted = source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED;
	// fprintf (stderr, "%d %d : %d : %f %f : %f %f\n", source->peak->x, source->peak->y, subtracted, source->peak->flux, source->pixels->data.F32[yc][xc], mcore, mpeak); 

	fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
		 source->peak->xf, source->peak->yf, 
		 model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], 
		 source->psfMag, source->apMag, source->errMag,
		 model->params->data.F32[PM_PAR_I0], 
		 model->params->data.F32[PM_PAR_SXX], model->params->data.F32[PM_PAR_SXY], model->params->data.F32[PM_PAR_SYY], 
		 model->params->data.F32[PM_PAR_7]);
    }
    fclose (f);
    return true;
}

// XXX replace this with a call to a pmConfig function (pmConfigDump...)
bool psphotDumpConfig (pmConfig *config) {

  psMetadataConfigWrite (config->user, "user.md");
  psMetadataConfigWrite (config->camera, "camera.md");
  psMetadataConfigWrite (config->recipes, "recipes.md");
  psMetadataConfigWrite (config->arguments, "arguments.md");
  psMetadataConfigWrite (config->files, "files.md");
  return true;
}

int psphotSaveImage (psMetadata *header, psImage *image, char *filename) {

    psFits *fits = psFitsOpen (filename, "w");
    psFitsWriteImage (fits, NULL, image, 0, NULL);
    psFitsClose (fits);
    return (TRUE);
}

bool psphotDumpMoments (psMetadata *recipe, psArray *sources) {

    bool status;

    // optional dump of all rough source data
    char *output = psMetadataLookupStr (&status, recipe, "MOMENTS_OUTPUT_FILE");
    if (!output) return false;
    if (output[0] == 0) return false;
    if (!strcasecmp (output, "NONE")) return false;

    pmMomentsWriteText (sources, output);
    return true;
}

bool psphotDumpSource (pmSource *source, char *name) {

    FILE *f = fopen (name, "w");
    if (f == NULL) psAbort("can't open file");

    for (int i = 0; i < source->pixels->numRows; i++) {
        for (int j = 0; j < source->pixels->numCols; j++) {
            // skip masked points
            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) {
                continue;
            }
            // skip zero-variance points
            if (source->variance->data.F32[i][j] == 0) {
                continue;
            }

            fprintf (f, "%d %d %f %f %d\n",
                     (j + source->pixels->col0),
                     (i + source->pixels->row0),
                     source->pixels->data.F32[i][j],
                     1.0 / source->variance->data.F32[i][j],
                     source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]);
        }
    }
    fclose (f);
    return true;
}

bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view) {

    bool status = false;

    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");

    // loop over the available readouts
    for (int i = 0; i < num; i++) {
	if (!psphotAddPhotcodeReadout (config, view, "PSPHOT.INPUT", i)) {
            psError (PSPHOT_ERR_CONFIG, false, "failed to add photcode to PSPHOT.INPUT entry %d", i);
	    return false;
	}
    }
    return true;
}

bool psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) {

    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    PS_ASSERT (file, false);

    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);

    // determine PHOTCODE from fpa & view, overwrite in readout->analysis
    char *photcode = pmConceptsPhotcodeForView (file, view);
    PS_ASSERT (photcode, false);

    psMetadataAddStr (readout->analysis, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "photcode from FPA concepts", photcode);
    psLogMsg ("psphot", 3, "PHOTCODE is %s", photcode);

    psFree (photcode);
    return true;
}

bool psphotSetHeaderNstars (psMetadata *header, psArray *sources) {

    int nSrc = 0;
    int nCR  = 0;
    int nEXT = 0;

    // count the number of sources which will be written and other sub-types
    for (int i = 0; (sources != NULL) && (i < sources->n); i++) {
        pmSource *source = (pmSource *) sources->data[i];
        pmModel *model = pmSourceGetModel (NULL, source);
        if (model == NULL)
            continue;
        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
            nEXT ++;
        }
        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) {
            nCR ++;
        }
        nSrc ++;
    }

    for (int i = 0; (sources != NULL) && (i < sources->n); i++) {
        pmSource *source = (pmSource *) sources->data[i];
        pmModel *model = pmSourceGetModel (NULL, source);
        if (model == NULL)
            continue;
        nSrc ++;
    }

    psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS",   PS_META_REPLACE, "Number of sources", nSrc);
    psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_EXT", PS_META_REPLACE, "Number of extended sources", nEXT);
    psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_CR",  PS_META_REPLACE, "Number of cosmic rays", nCR);
    return true;
}

// these values are saved in an output header stub - they are added to either the PHU header
// (CMP) or the MEF table header (CMF).  before the header is created, each readout has these
// values stored on readout->analysis
psMetadata *psphotDefineHeader (psMetadata *analysis) {

    bool status = true;

    psMetadata *header = psMetadataAlloc ();

    // write necessary information to output header
    psMetadataItemSupplement (&status, header, analysis, "ZERO_PT");
    psMetadataItemSupplement (&status, header, analysis, "PHOTCODE");

    psMetadataItemSupplement (&status, header, analysis, "APMIFIT");
    psMetadataItemSupplement (&status, header, analysis, "DAPMIFIT");
    psMetadataItemSupplement (&status, header, analysis, "NAPMIFIT");

    // PSF model parameters (shape values for image center)
    psMetadataItemSupplement (&status, header, analysis, "NPSFSTAR");
    psMetadataItemSupplement (&status, header, analysis, "APLOSS");

    psMetadataItemSupplement (&status, header, analysis, "FWHM_MAJ");
    psMetadataItemSupplement (&status, header, analysis, "FW_MJ_SG");
    psMetadataItemSupplement (&status, header, analysis, "FW_MJ_LQ");
    psMetadataItemSupplement (&status, header, analysis, "FW_MJ_UQ");

    psMetadataItemSupplement (&status, header, analysis, "FWHM_MIN");
    psMetadataItemSupplement (&status, header, analysis, "FW_MN_SG");
    psMetadataItemSupplement (&status, header, analysis, "FW_MN_LQ");
    psMetadataItemSupplement (&status, header, analysis, "FW_MN_UQ");

    psMetadataItemSupplement (&status, header, analysis, "ANGLE");

    // Image Quality measurements
    psMetadataItemSupplement (&status, header, analysis, "IQ_NSTAR");

    psMetadataItemSupplement (&status, header, analysis, "IQ_FW1");
    psMetadataItemSupplement (&status, header, analysis, "IQ_FW1_E");
    psMetadataItemSupplement (&status, header, analysis, "IQ_FW2");
    psMetadataItemSupplement (&status, header, analysis, "IQ_FW2_E");

    psMetadataItemSupplement (&status, header, analysis, "IQ_M2");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M2_ER");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M2_LQ");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M2_UQ");

    psMetadataItemSupplement (&status, header, analysis, "IQ_M2C");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M2C_E");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M2C_L");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M2C_U");

    psMetadataItemSupplement (&status, header, analysis, "IQ_M2S");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M2S_E");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M2S_L");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M2S_U");

    psMetadataItemSupplement (&status, header, analysis, "IQ_M3");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M3_ER");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M3_LQ");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M3_UQ");

    psMetadataItemSupplement (&status, header, analysis, "IQ_M4");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M4_ER");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M4_LQ");
    psMetadataItemSupplement (&status, header, analysis, "IQ_M4_UQ");

    // XXX these need to be defined from elsewhere
    psMetadataAdd (header, PS_LIST_TAIL, "FSATUR",   PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG",      0.0);
    psMetadataAdd (header, PS_LIST_TAIL, "FLIMIT",   PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG",    0.0);
    psMetadataItemSupplement (&status, header, analysis, "NSTARS");

    psMetadataItemSupplement (&status, header, analysis, "NDET_EXT");
    psMetadataItemSupplement (&status, header, analysis, "NDET_CR");

    // sky background model statistics
    psMetadataItemSupplement (&status, header, analysis, "MSKY_MN");
    psMetadataItemSupplement (&status, header, analysis, "MSKY_SIG");
    psMetadataItemSupplement (&status, header, analysis, "MSKY_MIN");
    psMetadataItemSupplement (&status, header, analysis, "MSKY_MAX");
    psMetadataItemSupplement (&status, header, analysis, "MSKY_NX");
    psMetadataItemSupplement (&status, header, analysis, "MSKY_NY");

    psMetadataAddF32 (header, PS_LIST_TAIL, "DT_PHOT", PS_META_REPLACE, "elapsed psphot time", psTimerMark ("psphotReadout"));

    return header;
}

// XXX add args as needed
bool psphotDumpPSFStars (pmReadout *readout, pmPSFtry *try, float radius, psImageMaskType maskVal, psImageMaskType markVal) {

    psphotSaveImage (NULL, readout->image,  "rawstars.fits");

    for (int i = 0; i < try->sources->n; i++) {
        // masked for: bad model fit, outlier in parameters
        if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL)
            continue;

        pmSource *source = try->sources->data[i];
        float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
        float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];

        // set the mask and subtract the PSF model
        // XXX should we be using maskObj? should we be unsetting the mask?
        // use pmModelSub because modelFlux has not been generated
        assert (source->maskObj);
        psImageKeepCircle (source->maskObj, x, y, radius, "OR", markVal);
        pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    }

    FILE *f = fopen ("shapes.dat", "w");
    for (int i = 0; i < try->sources->n; i++) {
        psF32 inPar[10];  // must be psF32 to pmPSF_FitToModel

        // masked for: bad model fit, outlier in parameters
        if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;

        pmSource *source = try->sources->data[i];
        psF32 *outPar = source->modelEXT->params->data.F32;

        psEllipseShape shape;

        shape.sx  = outPar[PM_PAR_SXX] / M_SQRT2;
        shape.sy  = outPar[PM_PAR_SYY] / M_SQRT2;
        shape.sxy = outPar[PM_PAR_SXY];

        psEllipsePol pol = pmPSF_ModelToFit (outPar);
        inPar[PM_PAR_E0] = pol.e0;
        inPar[PM_PAR_E1] = pol.e1;
        inPar[PM_PAR_E2] = pol.e2;
        pmPSF_FitToModel (inPar, 0.1);

        psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0);
        psEllipseAxes axes2 = psEllipsePolToAxes(pol, 0.1);

        psEllipsePol pol2 = psEllipseAxesToPol (axes1);

        fprintf (f, "%3d  %7.2f %7.2f  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f  :  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n",
                 i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS],
                 outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY],
                 pol.e0, pol.e1, pol.e2,
                 pol2.e0, pol2.e1, pol2.e2,
                 inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY],
                 axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD,
                 axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD
            );
    }
    fclose (f);

    psphotSaveImage (NULL, readout->image,  "psfstars.fits");
    pmSourcesWritePSFs (try->sources, "psfstars.dat");
    pmSourcesWriteEXTs (try->sources, "extstars.dat", false);
    // XXX need alternative output function
    // psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
    // psMetadataConfigWrite (psfData, "psfmodel.dat");
    psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n");

    return true;
}
