Index: /trunk/psphot/src/Makefile.am
===================================================================
--- /trunk/psphot/src/Makefile.am	(revision 14962)
+++ /trunk/psphot/src/Makefile.am	(revision 14963)
@@ -68,4 +68,5 @@
 	psphotSourceSize.c	 \
 	psphotDiagnosticPlots.c	 \
+	psphotMakeFluxScale.c	 \
 	psphotAddNoise.c
 
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 14962)
+++ /trunk/psphot/src/psphot.h	(revision 14963)
@@ -168,3 +168,5 @@
 bool psphotDiagnosticPlots (pmConfig *config, char *name, ...);
 
+bool psphotMakeFluxScale (psImage *image, psMetadata *recipe, pmPSF *psf);
+
 #endif
Index: /trunk/psphot/src/psphotMakeFluxScale.c
===================================================================
--- /trunk/psphot/src/psphotMakeFluxScale.c	(revision 14963)
+++ /trunk/psphot/src/psphotMakeFluxScale.c	(revision 14963)
@@ -0,0 +1,63 @@
+# 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;
+}