# include "psphotInternal.h"

bool psphotMakeFluxScale (psImage *image, psMetadata *recipe, pmPSF *psf) {

    bool status = false;

    psTimerStart ("residuals");

    int Nx = psMetadataLookupS32(&status, recipe, "PSF.FLUXSCALE.NX");
    int Ny = psMetadataLookupS32(&status, recipe, "PSF.FLUXSCALE.NY");

    // stats structure for use by ApTrend : XXX make parameters user setable
    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);

    pmTrend2D *trend = pmTrend2DAlloc (PM_TREND_MAP, image, Nx, Ny, stats);

    psVector *xPts = psVectorAlloc (Nx*Ny, PS_TYPE_F32);
    psVector *yPts = psVectorAlloc (Nx*Ny, PS_TYPE_F32);
    psVector *fPts = psVectorAlloc (Nx*Ny, PS_TYPE_F32);

    int Npts = 0;
    bool success = true;                // Function succeeded?

    // generate a set of test normalized PSF fluxes filling the grid
    for (int ix = 0; ix < Nx; ix++) {
        for (int iy = 0; iy < Ny; iy++) {

            float x = psImageBinningGetFineX (trend->map->binning, ix + 0.5);
            float y = psImageBinningGetFineY (trend->map->binning, iy + 0.5);

            // create normalized model object at xc,yc
            pmModel *model = pmModelFromPSFforXY (psf, x, y, 1.0);
            if (!model) {
                psError(PS_ERR_UNKNOWN, false, "Unable to generate model for grid point %d,%d", ix, iy);
                success = false;
                goto DONE;
            }

            // measure the fitMag for this model
            float fitSum = model->modelFlux (model->params);
            assert (fitSum > 0);
            assert (isfinite(fitSum));

            xPts->data.F32[Npts] = x;
            yPts->data.F32[Npts] = y;
            fPts->data.F32[Npts] = fitSum;
            Npts ++;
            assert (Npts <= xPts->nalloc);
            psFree (model);
        }
    }
    xPts->n = Npts;
    yPts->n = Npts;
    fPts->n = Npts;

    if (!pmTrend2DFit (trend, NULL, 0xff, xPts, yPts, fPts, NULL)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to fit trend");
        success = false;
        goto DONE;
    }

    // XXX do something useful to measure residual statistics

    psf->FluxScale = psMemIncrRefCounter(trend);

 DONE:
    psFree(xPts);
    psFree(yPts);
    psFree(fPts);
    psFree(stats);
    psFree(trend);

    return success;
}
