# 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;

    // 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);

	    // 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;
    
    // XXX test for errors here
    pmTrend2DFit (trend, NULL, 0xff, xPts, yPts, fPts, NULL);
    
    // XXX do something useful to measure residual statistics

    psf->FluxScale = trend;

    psFree (xPts);
    psFree (yPts);
    psFree (fPts);
    psFree (stats);

    return true;
}
