/** @file  pmPSFtryFitPSF.c
 *  @brief Fit the PSF model to the candidate PSF stars (part of PSF model generation)
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.69 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-01-27 06:39:38 $
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */

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

# include <pslib.h>
#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPAMaskWeight.h"
#include "pmSpan.h"
#include "pmFootprint.h"
#include "pmPeaks.h"
#include "pmMoments.h"
#include "pmResiduals.h"
#include "pmGrowthCurve.h"
#include "pmTrend2D.h"
#include "pmPSF.h"
#include "pmModel.h"
#include "pmSource.h"
#include "pmSourceUtils.h"
#include "pmPSFtry.h"
#include "pmModelClass.h"
#include "pmModelUtils.h"
#include "pmSourceFitModel.h"
#include "pmSourcePhotometry.h"
#include "pmSourceVisual.h"

// stage 3: Refit with fixed shape parameters.  This function uses the LMM fitting, but could
// be re-written to use the simultaneous linear fitting (see psphotFitSourcesLinear.c)
bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) {

    bool status;

    psTimerStart ("psf.fit");

    // maskVal is used to test for rejected pixels, and must include markVal
    maskVal |= markVal;

    // DEBUG code: save the PSF model fit data in detail
# ifdef DEBUG
    char filename[64];
    snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy);
    FILE *f = fopen (filename, "w");
    psAssert (f, "failed open");
# endif

    int Npsf = 0;
    for (int i = 0; i < psfTry->sources->n; i++) {

        pmSource *source = psfTry->sources->data[i];
	psAssert (source->moments, "how can a psf source not have moments?");

        // masked for: bad model fit, outlier in parameters
        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) {
            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y);
            continue;
        }

        // set shape for this model based on PSF
	psFree (source->modelPSF);
        source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
        if (source->modelPSF == NULL) {
            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
            continue;
        }
	// PSF fit and aperture mags use different radii
        source->modelPSF->fitRadius = options->fitRadius;
        source->apRadius = options->apRadius;

        // set object mask to define valid pixels for PSF model fit
        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);

        // fit the PSF model to the source
        status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_NORM, maskVal);

        // skip poor fits
        if (!status) {
            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
            continue;
        }

        // set object mask to define valid pixels for APERTURE magnitude
	if (options->fitRadius != options->apRadius) {
            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
	    psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
	}

	// This function calculates the psf and aperture magnitudes
        status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal); // raw PSF mag, AP mag
        if (!status || isnan(source->apMag)) {
            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
            continue;
        }

        // clear object mask to define valid pixels
	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask

        psfTry->fitMag->data.F32[i] = source->psfMag;
        psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
        psfTry->metricErr->data.F32[i] = source->errMag;

	// XXX this did not work: modifies shape of psf too much
        // psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag;

# ifdef DEBUG
	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, 
		 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], 
		 source->psfMag, source->apMag, source->errMag,
		 source->modelPSF->params->data.F32[PM_PAR_I0], 
		 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], 
		 source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_7]);
# endif

        psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
        Npsf ++;
    }
    psfTry->psf->nPSFstars = Npsf;

# ifdef DEBUG
    fclose (f);
# endif

    pmSourceVisualShowModelFits (psfTry->psf, psfTry->sources, maskVal);

    psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, psfTry->sources->n);
    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, psfTry->sources->n);

    if (Npsf == 0) {
        psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built.");
        return false;
    }

    return true;
}
