Index: /trunk/psphot/Makefile
===================================================================
--- /trunk/psphot/Makefile	(revision 5827)
+++ /trunk/psphot/Makefile	(revision 5828)
@@ -21,43 +21,42 @@
 
 PSPHOT = \
-$(SRC)/psphot.$(ARCH).o            \
-$(SRC)/psphotArguments.$(ARCH).o   \
-$(SRC)/psphotSetup.$(ARCH).o	   \
-$(SRC)/psphotImageStats.$(ARCH).o  \
-$(SRC)/psphotSourceStats.$(ARCH).o \
-$(SRC)/psphotChoosePSF.$(ARCH).o   \
-$(SRC)/psphotApplyPSF.$(ARCH).o	   \
-$(SRC)/psphotFixedPSF.$(ARCH).o	   \
+$(SRC)/psphot.$(ARCH).o             \
+$(SRC)/psphotArguments.$(ARCH).o    \
+$(SRC)/psphotSetup.$(ARCH).o	    \
+$(SRC)/psphotImageStats.$(ARCH).o   \
+$(SRC)/psphotSourceStats.$(ARCH).o  \
+$(SRC)/psphotChoosePSF.$(ARCH).o    \
+$(SRC)/psphotApplyPSF.$(ARCH).o	    \
+$(SRC)/psphotFixedPSF.$(ARCH).o	    \
 $(SRC)/psphotEnsemblePSF.$(ARCH).o  \
-$(SRC)/psphotReapplyPSF.$(ARCH).o  \
-$(SRC)/psphotFitGalaxies.$(ARCH).o \
-$(SRC)/psphotOutput.$(ARCH).o      \
-$(SRC)/psphotMarkPSF.$(ARCH).o     \
-$(SRC)/psphotSortBySN.$(ARCH).o    \
-$(SRC)/psphotDefinePixels.$(ARCH).o\
-$(SRC)/psphotMagnitudes.$(ARCH).o  \
-$(SRC)/psphotImageBackground.$(ARCH).o  \
-$(SRC)/psphotBasicDeblend.$(ARCH).o  \
+$(SRC)/psphotReapplyPSF.$(ARCH).o   \
+$(SRC)/psphotFitGalaxies.$(ARCH).o  \
+$(SRC)/psphotOutput.$(ARCH).o       \
+$(SRC)/psphotMarkPSF.$(ARCH).o      \
+$(SRC)/psphotSortBySN.$(ARCH).o     \
+$(SRC)/psphotDefinePixels.$(ARCH).o \
+$(SRC)/psphotMagnitudes.$(ARCH).o   \
+$(SRC)/psphotRadiusChecks.$(ARCH).o \
+$(SRC)/psphotReplaceUnfit.$(ARCH).o \
+$(SRC)/psphotEvalPSF.$(ARCH).o      \
+$(SRC)/psphotEvalFLT.$(ARCH).o      \
+$(SRC)/psphotFullFit.$(ARCH).o      \
+$(SRC)/pmSourceContour.$(ARCH).o    \
+$(SRC)/psLine.$(ARCH).o		    \
+$(SRC)/psModulesUtils.$(ARCH).o	    \
+$(SRC)/pmPeaksSigmaLimit.$(ARCH).o  \
+$(SRC)/pmSourceFitFixed.$(ARCH).o   \
+$(SRC)/psSparse.$(ARCH).o           \
+$(SRC)/psImageData.$(ARCH).o        \
+$(SRC)/psphotModelTest.$(ARCH).o    \
+$(SRC)/psphotImageBackground.$(ARCH).o \
+$(SRC)/psphotBasicDeblend.$(ARCH).o    \
 $(SRC)/psphotModelGroupInit.$(ARCH).o  \
-$(SRC)/psphotRadiusChecks.$(ARCH).o  \
-$(SRC)/psphotReplaceUnfit.$(ARCH).o  \
-$(SRC)/psphotEvalPSF.$(ARCH).o  \
-$(SRC)/psphotEvalFLT.$(ARCH).o  \
-$(SRC)/psphotFullFit.$(ARCH).o  \
-$(SRC)/pmSourceContour.$(ARCH).o \
-$(SRC)/psLine.$(ARCH).o		   \
-$(SRC)/psModulesUtils.$(ARCH).o	   \
-$(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
-$(SRC)/pmSourceFitFixed.$(ARCH).o  \
-$(SRC)/psSparse.$(ARCH).o            \
-$(SRC)/psImageData.$(ARCH).o        \
-$(SRC)/psphotApResid.$(ARCH).o
-
-PSMODULES = \
-$(SRC)/psEllipse.$(ARCH).o         \
-$(SRC)/pmPSF.$(ARCH).o             \
-$(SRC)/pmPSFtry.$(ARCH).o          \
-$(SRC)/pmModelGroup.$(ARCH).o      \
-$(SRC)/pmObjects_EAM.$(ARCH).o
+$(SRC)/pmModelFitSet.$(ARCH).o  \
+$(SRC)/pmSourceFitSet.$(ARCH).o        \
+$(SRC)/psphotSourceFits.$(ARCH).o        \
+$(SRC)/psphotBlendFit.$(ARCH).o        \
+$(SRC)/psphotApResid.$(ARCH).o \
+$(SRC)/psBicube.$(ARCH).o 
 
 MODELS = \
@@ -91,5 +90,5 @@
 psphot: $(BIN)/psphot.$(ARCH)
 $(BIN)/psphot.$(ARCH) : $(PSPHOT)
-$(PSPHOT) $(PSMODULES): $(SRC)/psphot.h
+$(PSPHOT) : $(SRC)/psphot.h
 
 modeltest.install: psphotModelTest.install
Index: unk/psphot/src/modelTestFitSource.c
===================================================================
--- /trunk/psphot/src/modelTestFitSource.c	(revision 5827)
+++ 	(revision )
@@ -1,132 +1,0 @@
-# include "psphot.h"
-
-int DumpImage (psImage *image, char *filename);
-
-bool modelTestFitSource (eamReadout *imdata, psMetadata *config) {
-
-    bool status;
-    int modelType;
-    float obsMag, fitMag, value;
-    char name[64];
-
-    psMetadataItem *item  = NULL;
-
-    // find the model: supplied by user or first in the PSF_MODEL list
-    char *modelName  = psMetadataLookupPtr (&status, config, "TEST_FIT_MODEL");
-    if (modelName == NULL) {
-	// get the list pointers for the PSF_MODEL entries
-	psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL");
-	if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");
-
-	// take the first list element
-	psList *list = (psList *) mdi->data.list;
-	item = psListGet (list, PS_LIST_HEAD);
-	modelName = item->data.V;
-    }
-    modelType = pmModelSetType (modelName);
-    if (modelType < 0) psAbort ("fitsource", "unknown model %s", modelName);
-
-    // find the fitting parameters (try test values first)
-    float INNER = psMetadataLookupF32 (&status, config, "TEST_FIT_INNER_RADIUS");
-    if (!status) {
-	INNER = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS");
-    }
-
-    float OUTER = psMetadataLookupF32 (&status, config, "TEST_FIT_OUTER_RADIUS");
-    if (!status) {
-	OUTER = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
-    }
-
-    float RADIUS = psMetadataLookupF32 (&status, config, "TEST_FIT_RADIUS");
-    if (!status) {
-	RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
-    }
-
-    float mRADIUS = psMetadataLookupF32 (&status, config, "TEST_MOMENTS_RADIUS");
-    if (!status) {
-	mRADIUS = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS");
-    }
-
-    // define the source of interest
-    float xObj     = psMetadataLookupF32 (&status, config, "TEST_FIT_X");
-    float yObj     = psMetadataLookupF32 (&status, config, "TEST_FIT_Y");
-
-    // construct the source structures
-    pmSource *source = pmSourceAlloc();
-    source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
-    psphotDefinePixels (source, imdata, xObj, yObj, OUTER);
-    DumpImage (source->pixels, "object.fits");
-
-    // find the local sky
-    status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
-    if (!status) psAbort ("fitsource", "pmSourceLocalSky error");
-
-    // get the source moments
-    status = pmSourceMoments (source, mRADIUS);
-    if (!status) psAbort ("fitsource", "pmSourceLocalSky error");
-    source->peak->counts = source->moments->Peak;
-
-    fprintf (stderr, "sum: %f\n", source->moments->Sum);
-
-    // get the initial model parameter guess
-    pmModel *model = pmSourceModelGuess (source, modelType);
-    // if any parameters are defined, use those values
-    int nParams = pmModelParameterCount (modelType);
-    psF32 *params = model->params->data.F32;
-    for (int i = 0; i < nParams; i++) {
-	if (i == 2) {
-	    params[i] = xObj;
-	    continue;
-	}
-	if (i == 3) {
-	    params[i] = yObj;
-	    continue;
-	}
-	sprintf (name, "TEST_FIT_PAR%d", i);
-	value = psMetadataLookupF32 (&status, config, name);
-	if (status) {
-	    fprintf (stderr, "using supplied value %f for PAR %d\n", value, i);
-	    params[i] = value;
-	} else {
-	    fprintf (stderr, "using guessed  value %f for PAR %d\n", params[i], i);
-	}
-    }
-
-    float area = params[4]*params[5];
-    fprintf (stderr, "peak: %f\n", source->moments->Sum*area);
-
-    // what fitting mode to use?
-    char *fitModeWord = psMetadataLookupPtr (&status, config, "TEST_FIT_MODE");
-    if (!status) {
-	fitModeWord = psStringCopy ("FLT");
-    }
-    bool fitMode = !strcasecmp (fitModeWord, "PSF");
-
-    // define the pixels used for the fit
-    psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PSPHOT_MASK_MARKED);
-    status = pmSourceFitModel (source, model, fitMode);
-
-    // measure the source mags
-    pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);
-    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag);
-
-    // subtract object, leave local sky
-    pmSourceSubModel (source->pixels, source->mask, model, false, false);
-    
-    for (int i = 0; i < nParams; i++) {
-	fprintf (stderr, "result value %f for PAR %d\n", params[i], i);
-    }
-
-    // write out 
-    DumpImage (source->pixels, "resid.fits");
-    DumpImage (source->mask, "mask.fits");
-    return true;
-}
-
-int DumpImage (psImage *image, char *filename) {
-
-    psFits *fits = psFitsOpen (filename, "w");
-    psFitsWriteImage (fits, NULL, image, 0);
-    psFitsClose (fits);
-    return (TRUE);
-}
Index: /trunk/psphot/src/pmModelFitSet.c
===================================================================
--- /trunk/psphot/src/pmModelFitSet.c	(revision 5828)
+++ /trunk/psphot/src/pmModelFitSet.c	(revision 5828)
@@ -0,0 +1,56 @@
+# include "psphot.h"
+
+// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2, 
+// nPar = nSrc*(nOnePar - 1) + 1
+
+static pmModelFunc mFunc;
+static int nPar;
+static psVector *onePar;
+static psVector *oneDeriv;
+
+bool pmModelFitSetInit (pmModelType type) {
+
+    mFunc = pmModelFunc_GetFunction (type);
+    nPar  = pmModelParameterCount (type);
+
+    onePar = psVectorAlloc (nPar, PS_DATA_F32);
+    oneDeriv = psVectorAlloc (nPar, PS_DATA_F32);
+
+    return true;
+}
+
+psF32 pmModelFitSet(psVector *deriv,
+		    const psVector *params,
+		    const psVector *x)
+{
+
+    psF32 value;
+    psF32 model;
+
+    psF32 *PAR = onePar->data.F32;
+    psF32 *dPAR = oneDeriv->data.F32;
+
+    int nSrc = (params->n - 1) / (nPar - 1);
+
+    model = params->data.F32[0];
+
+    PAR[0] = 0;
+    for (int i = 0; i < nSrc; i++) {
+	for (int n = 1; n < nPar; n++) {
+	    PAR[n] = params->data.F32[i*nPar - i + n];
+	}
+	value = mFunc (oneDeriv, onePar, x);
+	model += value;
+	for (int n = 1; n < nPar; n++) {
+	    deriv->data.F32[i*nPar - i + n] = dPAR[n];
+	}
+    }
+    return (model);
+}
+
+/*
+i:         0           1               2 
+n:         1  2  3  4  5  6  1  2  3  4  5  6  1  2  3  4  5  6
+i*6 + n: 0 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18
+*/
+
Index: /trunk/psphot/src/pmPeaksSigmaLimit.c
===================================================================
--- /trunk/psphot/src/pmPeaksSigmaLimit.c	(revision 5827)
+++ /trunk/psphot/src/pmPeaksSigmaLimit.c	(revision 5828)
@@ -26,5 +26,5 @@
     threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; 
     // threshold = NSIGMA*sky->sampleStdev;
-    psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold);
+    psLogMsg ("psphot", 4, "threshold: %f DN\n", threshold);
 
     // find the peaks in the smoothed image 
Index: /trunk/psphot/src/pmSourceFitSet.c
===================================================================
--- /trunk/psphot/src/pmSourceFitSet.c	(revision 5828)
+++ /trunk/psphot/src/pmSourceFitSet.c	(revision 5828)
@@ -0,0 +1,224 @@
+# include "psphot.h"
+
+bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, bool PSF) {
+
+    double x, y, Io;
+
+    FILE *f = fopen (fitset, "r");
+    if (f == NULL) return false;
+
+    psArray *modelSet = psArrayAlloc (16);
+    modelSet->n = 0;
+
+    while (fscanf (f, "%lf %lf %lf", &x, &y, &Io) == 3) {
+	pmModel *model = pmModelAlloc (oneModel->type);
+	
+	for (psS32 i = 0; i < model->params->n; i++) {
+	    model->params->data.F32[i] = oneModel->params->data.F32[i];
+	    model->dparams->data.F32[i] = oneModel->dparams->data.F32[i];
+	}
+	model->params->data.F32[1] = Io;
+	model->params->data.F32[2] = x;
+	model->params->data.F32[3] = y;
+	psArrayAdd (modelSet, 16, model);
+    }
+
+    pmSourceFitSet (source, modelSet, PSF);
+
+    // write out positive object
+    psphotSaveImage (NULL, source->pixels, "object.fits");
+
+    // subtract object, leave local sky
+    for (int i = 0; i < modelSet->n; i++) {
+	pmModel *model = modelSet->data[i];
+	pmSourceSubModel (source->pixels, source->mask, model, false, false);
+    
+	fprintf (stderr, "output parameters (obj %d):\n", i);
+	for (int n = 0; n < model->params->n; n++) {
+	    fprintf (stderr, "%d : %f\n", n, model->params->data.F32[n]);
+	}
+    }
+
+    // write out 
+    psphotSaveImage (NULL, source->pixels, "resid.fits");
+    psphotSaveImage (NULL, source->mask, "mask.fits");
+    return true;
+}
+
+# define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
+# define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
+
+bool pmSourceFitSet (pmSource *source,
+		     psArray *modelSet,
+		     const bool PSF)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->moments, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->weight, false);
+
+    psBool fitStatus = true;
+    psBool onPic     = true;
+    psBool rc        = true;
+
+    // base values on first model
+    pmModel *model = modelSet->data[0];
+
+    // set the static variables
+    pmModelFitSetInit (model->type);
+
+    int nSrc = modelSet->n;
+    int nPar = model->params->n - 1;  // number of object parameters (excluding sky)
+
+    psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
+    psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
+    psVector *paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL;
+
+    // all but the sky are allowed to vary independently (subject to PSF)
+    params->data.F32[0] = model->params->data.F32[0];
+    for (int i = 0; i < nSrc; i++) {
+      model = modelSet->data[i];
+      for (int n = 1; n < nPar + 1; n++) {
+	params->data.F32[i*nPar + n] = model->params->data.F32[n];
+	dparams->data.F32[i*nPar + n] = model->dparams->data.F32[n];
+	if (PSF) {
+	    paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1;
+	}
+      }
+    }
+
+    // PSF model only fits first 4 parameters, FLT model fits all
+    int nParams = PSF ? nSrc*3 + 1 : nSrc*nPar + 1;
+
+    // maximum number of valid pixels
+    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
+
+    // construct the coordinate and value entries
+    psArray *x = psArrayAlloc(nPix);
+    psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);
+    psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32);
+
+    nPix = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            if (source->mask->data.U8[i][j]) {
+		continue;
+	    }
+	    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+	    // Convert i/j to image space:
+	    coord->data.F32[0] = (psF32) (j + source->pixels->col0);
+	    coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+	    x->data[nPix] = (psPtr *) coord;
+	    y->data.F32[nPix] = source->pixels->data.F32[i][j];
+
+	    // psMinimizeLMChi2 takes wt = 1/dY^2
+	    if (source->weight->data.F32[i][j] == 0) {
+		continue;
+	    } 
+	    yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
+	    nPix++;
+        }
+    }
+    if (nPix <  nParams + 1) {
+        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        model->status = PM_MODEL_BADARGS;
+	psFree (x);
+	psFree (y);
+	psFree (yErr);
+        return(false);
+    }
+    x->n = nPix;
+    y->n = nPix;
+    yErr->n = nPix;
+
+    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
+                            PM_SOURCE_FIT_MODEL_TOLERANCE);
+
+    // Set the parameter range checks
+    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
+    psVector *beta_lim = NULL;
+    psVector *params_min = NULL;
+    psVector *params_max = NULL;
+
+    // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
+    //           in the SDRS, I define a new psMinimization which will take these in
+    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
+    modelLimits (&beta_lim, &params_min, &params_max);
+    for (int i = 0; i < nSrc; i++) {
+	for (int n = 1; n < nPar; n++) {
+	    covar->data.F64[0][i*nPar + n] = beta_lim->data.F32[n];
+	    covar->data.F64[1][i*nPar + n] = params_min->data.F32[n];
+	    covar->data.F64[2][i*nPar + n] = params_max->data.F32[n];
+	}
+    }
+    covar->data.F64[0][0] = beta_lim->data.F32[0];
+    covar->data.F64[1][0] = params_min->data.F32[0];
+    covar->data.F64[2][0] = params_max->data.F32[0];
+
+    psTrace (".pmObjects.pmSourceFitSet", 5, "fitting function\n");
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, pmModelFitSet);
+
+    // parameter errors from the covariance matrix
+    for (int i = 0; i < dparams->n; i++) {
+        if ((paramMask != NULL) && paramMask->data.U8[i])
+            continue;
+        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
+    }
+
+    // get the Gauss-Newton distance for fixed model parameters
+    if (paramMask != NULL) {
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
+        psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, pmModelFitSet);
+        for (int i = 0; i < dparams->n; i++) {
+            if (!paramMask->data.U8[i])
+                continue;
+            dparams->data.F32[i] = delta->data.F64[i];
+        }
+    }
+
+    // assign back the parameters to the models
+    for (int i = 0; i < nSrc; i++) {
+      model = modelSet->data[i];
+      model->params->data.F32[0] = params->data.F32[0];
+      for (int n = 1; n < nPar + 1; n++) {
+	  model->params->data.F32[n] = params->data.F32[i*nPar + n];
+	  model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n];
+      }
+      // save the resulting chisq, nDOF, nIter
+      // these are not unique for any one source
+      model->chisq = myMin->value;
+      model->nIter = myMin->iter;
+      model->nDOF  = y->n - nParams;
+
+      // set the model success or failure status
+      model->status = fitStatus ? PM_MODEL_SUCCESS : PM_MODEL_NONCONVERGE;
+
+      // models can go insane: reject these
+      onPic &= (model->params->data.F32[2] >= source->pixels->col0);
+      onPic &= (model->params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
+      onPic &= (model->params->data.F32[3] >= source->pixels->row0);
+      onPic &= (model->params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
+      if (!onPic) {
+	  model->status = PM_MODEL_OFFIMAGE;
+      }
+    }
+
+    source->mode |= PM_SOURCE_FITTED;
+
+    psFree(x);
+    psFree(y);
+    psFree(yErr);
+    psFree(myMin);
+    psFree(covar);
+    psFree(paramMask);
+
+    rc = (onPic && fitStatus);
+    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
+    return(rc);
+}
+
Index: /trunk/psphot/src/psBicube.c
===================================================================
--- /trunk/psphot/src/psBicube.c	(revision 5828)
+++ /trunk/psphot/src/psBicube.c	(revision 5828)
@@ -0,0 +1,143 @@
+# include "psphot.h"
+
+float poly3DTestValue (float x, float y, float z) {
+    float value = 5 + 3*x - 5*x*x - 4*x*y - 2*y*y + 2.5*z - 3*x*z - 6*z*z;
+    return (value);
+}
+
+
+float poly2DTestValue (float x, float y) {
+    float value = 5 + 3*x - 5*x*x - 4*x*y - 2*y*y;
+    return (value);
+}
+
+float bicubeTestValue (int x, int y) {
+    float value = 5 - 5*x*x - 2*y*y + 3*x - 4*x*y;
+    return (value);
+}
+
+bool bicubeTest () {
+
+    psImage *test = psImageAlloc (3, 3, PS_TYPE_F32);
+    for (int ny = 0; ny < 3; ny++) {
+	for (int nx = 0; nx < 3; nx++) {
+	    test->data.F32[ny][nx] = bicubeTestValue (nx - 1, ny - 1);
+	}
+    }
+    psPolynomial2D *poly = psImageBicubeFit (test, 1, 1);
+    psPlane peak = psImageBicubeMin (poly);
+    
+    fprintf (stderr, "done test: %f, %f\n", peak.x, peak.y);
+    return true;
+}
+
+bool poly3DTest_F64 () {
+
+    int Ncell = 50;
+    psVector *x = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F64);
+    psVector *y = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F64);
+    psVector *z = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F64);
+    psVector *f = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F64);
+
+    fprintf (stderr, "starting %s\n", __func__);
+    int cell = 0;
+    for (int nx = 0; nx < Ncell; nx++) {
+	for (int ny = 0; ny < Ncell; ny++) {
+	    for (int nz = 0; nz < Ncell; nz++) {
+		x->data.F64[cell] = nx;
+		y->data.F64[cell] = ny;
+		z->data.F64[cell] = nz;
+		f->data.F64[cell] = poly3DTestValue (nx, ny, nz);
+		cell++;
+	    }
+	}
+    }
+    fprintf (stderr, "built input data\n");
+
+    psPolynomial3D *poly = psPolynomial3DAlloc (2, 2, 2, PS_POLYNOMIAL_ORD);
+
+    psTimerStart ("fit");
+    poly = psVectorFitPolynomial3D (poly, NULL, 0, f, NULL, x, y, z);
+    fprintf (stderr, "fit 3d : %6.1f msec\n", 1000*psTimerMark ("fit"));
+
+    for (int nz = 0; nz <= poly->nZ; nz++) {
+	for (int ny = 0; ny <= poly->nY; ny++) {
+	    for (int nx = 0; nx <= poly->nX; nx++) {
+		fprintf (stderr, "%d %d %d : %f\n", nx, ny, nz, poly->coeff[nx][ny][nz]);
+	    }
+	}
+    }
+    return true;
+}
+
+bool poly3DTest_F32 () {
+
+    int Ncell = 50;
+    psVector *x = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F32);
+    psVector *y = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F32);
+    psVector *z = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F32);
+    psVector *f = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F32);
+
+    fprintf (stderr, "starting %s\n", __func__);
+    int cell = 0;
+    for (int nx = 0; nx < Ncell; nx++) {
+	for (int ny = 0; ny < Ncell; ny++) {
+	    for (int nz = 0; nz < Ncell; nz++) {
+		x->data.F32[cell] = nx;
+		y->data.F32[cell] = ny;
+		z->data.F32[cell] = nz;
+		f->data.F32[cell] = poly3DTestValue (nx, ny, nz);
+		cell++;
+	    }
+	}
+    }
+    fprintf (stderr, "built input data\n");
+
+    psPolynomial3D *poly = psPolynomial3DAlloc (2, 2, 2, PS_POLYNOMIAL_ORD);
+
+    psTimerStart ("fit");
+    poly = psVectorFitPolynomial3D (poly, NULL, 0, f, NULL, x, y, z);
+    fprintf (stderr, "fit 3d : %6.1f msec\n", 1000*psTimerMark ("fit"));
+
+    for (int nz = 0; nz <= poly->nZ; nz++) {
+	for (int ny = 0; ny <= poly->nY; ny++) {
+	    for (int nx = 0; nx <= poly->nX; nx++) {
+		fprintf (stderr, "%d %d %d : %f\n", nx, ny, nz, poly->coeff[nx][ny][nz]);
+	    }
+	}
+    }
+    return true;
+}
+
+bool poly2DTest () {
+
+    int Ncell = 50;
+    psVector *x = psVectorAlloc (Ncell*Ncell, PS_TYPE_F64);
+    psVector *y = psVectorAlloc (Ncell*Ncell, PS_TYPE_F64);
+    psVector *f = psVectorAlloc (Ncell*Ncell, PS_TYPE_F64);
+
+    fprintf (stderr, "starting __func__\n");
+    int cell = 0;
+    for (int nx = 0; nx < Ncell; nx++) {
+	for (int ny = 0; ny < Ncell; ny++) {
+	    x->data.F64[cell] = nx;
+	    y->data.F64[cell] = ny;
+	    f->data.F64[cell] = poly2DTestValue (nx, ny);
+	    cell++;
+	}
+    }
+    fprintf (stderr, "built input data\n");
+
+    psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD);
+
+    psTimerStart ("fit");
+    poly = psVectorFitPolynomial2D (poly, NULL, 0, f, NULL, x, y);
+    fprintf (stderr, "fit 2d : %6.1f msec\n", 1000*psTimerMark ("fit"));
+
+    for (int ny = 0; ny <= poly->nY; ny++) {
+	for (int nx = 0; nx <= poly->nX; nx++) {
+	    fprintf (stderr, "%d %d : %f\n", nx, ny, poly->coeff[nx][ny]);
+	}
+    }
+    return true;
+}
Index: /trunk/psphot/src/psModulesUtils.c
===================================================================
--- /trunk/psphot/src/psModulesUtils.c	(revision 5827)
+++ /trunk/psphot/src/psModulesUtils.c	(revision 5828)
@@ -110,6 +110,9 @@
     pmModel *new = pmModelAlloc (model->type);
     
-    new->chisq = model->chisq;
-    new->nIter = model->nIter;
+    new->chisq  = model->chisq;
+    new->nDOF   = model->nDOF;
+    new->nIter  = model->nIter;
+    new->status = model->status;
+    new->radius = model->radius;
 
     for (int i = 0; i < new->params->n; i++) {
Index: /trunk/psphot/src/psphot.c
===================================================================
--- /trunk/psphot/src/psphot.c	(revision 5827)
+++ /trunk/psphot/src/psphot.c	(revision 5828)
@@ -12,15 +12,24 @@
     bool         status;
 
+    psTimerStart ("complete");
+
     psphotModelGroupInit ();
 
     config = psphotArguments (&argc, argv);
 
-    // load input data (image and config)
-    // create or load mask and weight images
-    // we have memory leaks here -- may be from psMetadata
+    // load input data (config and images (signal, noise, mask)
+    // XXX we have memory leaks here -- may be from psMetadata
     imdata = psphotSetup (config);
+
     char *breakPt = psMetadataLookupPtr (&status, config, "BREAK_POINT");
-    if (!status) {
-      breakPt = psStringCopy ("NONE");
+    if (!status) breakPt = psStringCopy ("NONE");
+
+    int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE");
+    if (!status) FITMODE = 2;
+
+    // run model fitting tests on a single source
+    if (psMetadataLookupBool (&status, config, "TEST_FIT")) {
+	psphotModelTest (imdata, config);
+	exit (0);
     }
 
@@ -28,5 +37,6 @@
     sky = psphotImageStats (imdata, config);
 
-    // psPolynomial2D *skyModel = psphotImageBackground (imdata, config, sky);
+    // generate a background model (currently, 2D polynomial)
+    // XXX this should be available to be re-added to the original image
     psphotImageBackground (imdata, config, sky);
 
@@ -52,5 +62,4 @@
     // source analysis is done in S/N order (brightest first)
     sources = psArraySort (sources, psphotSortBySN);
-
     psphotDumpMoments (config, sources);
     if (!strcasecmp (breakPt, "CLASS")) exit (0);
@@ -59,10 +68,4 @@
     psf = psphotChoosePSF (config, sources, sky);
     if (!strcasecmp (breakPt, "PSFFIT")) exit (0);
-
-    // XXX add this as conditional output
-    psphotSamplePSFs (psf, imdata->image);
-
-    int FITMODE = psMetadataLookupS32 (&status, config, "FIT_MODE");
-    if (!status) FITMODE = 2;
 
     switch (FITMODE) {
@@ -79,7 +82,8 @@
 
       case 2:
-	// fit extended objects with galaxy models
-	psphotApplyPSF (imdata, config, sources, psf, sky);
-	psphotFitGalaxies (imdata, config, sources, sky);
+	psphotEnsemblePSF (imdata, config, sources, psf, sky);
+	psphotBlendFit (imdata, config, sources, psf, sky);
+	psphotReplaceUnfit (sources);
+	psphotApResid (sources, config, psf);
 	break;
 
@@ -91,5 +95,5 @@
       case 4:
 	// fit extended objects with galaxy models
-	psphotFixedPSF (imdata, config, sources, psf, sky);
+	psphotApplyPSF (imdata, config, sources, psf, sky);
 	psphotFitGalaxies (imdata, config, sources, sky);
 	break;
@@ -97,7 +101,8 @@
 
     // write out data in appropriate format
-    // measure aperture correction
-    // psphotApertureCorrection (sources, config, sky);
+    psphotSamplePSFs (config, psf, imdata->image);
     psphotOutput (imdata, config, sources, psf, sky);
+    psLogMsg ("psphot", 3, "wrote output: %f sec\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete"));
     exit (0);
 }
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 5827)
+++ /trunk/psphot/src/psphot.h	(revision 5828)
@@ -63,6 +63,4 @@
 
 // psphotModelTest functions
-psMetadata     *modelTestArguments (int *argc, char **argv);
-bool 	        modelTestFitSource (eamReadout *imdata, psMetadata *config);
 bool 	        psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
 float           psphotCrossProduct (pmSource *Mi, pmSource *Mj);
@@ -72,17 +70,44 @@
 pmModel        *pmModelCopy (pmModel *model);
 psArray        *pmSourceContour_EAM (psImage *image, int x, int y, float threshold);
-psMetadata     *psphotTestArguments (int *argc, char **argv);
 bool            psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky);
 
 bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
 bool psphotInitLimitsPSF (psMetadata *config);
-bool psphotEvalPSF (pmSource *source);
-bool psphotEvalFLT (pmSource *source);
+bool psphotEvalPSF (pmSource *source, pmModel *model);
+bool psphotEvalDBL (pmSource *source, pmModel *model);
+bool psphotEvalFLT (pmSource *source, pmModel *model);
 bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type);
-bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source);
+bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source, pmModel *model);
 bool psphotInitRadiusFLT (psMetadata *config, psStats *sky, pmModelType type);
-bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source);
-bool psphotSamplePSFs (pmPSF *psf, psImage *image);
+bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source, pmModel *model);
+bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image);
 bool psphotReplaceUnfit (psArray *sources);
 bool psphotDumpMoments (psMetadata *config, psArray *sources);
 bool psphotApResid (psArray *sources, psMetadata *config, pmPSF *psf);
+bool psphotWritePSF (pmPSF *psf, char *filename);
+pmPSF *psphotReadPSF (char *filename);
+
+bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...);
+bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *poly, char *format, ...);
+psPolynomial2D *psPolynomial2DfromMD (psMetadata *folder);
+psPolynomial3D *psPolynomial3DfromMD (psMetadata *folder);
+bool psphotModelTest (eamReadout *imdata, psMetadata *config);
+
+bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF);
+bool pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF);
+psF32 pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x);
+bool pmModelFitSetInit (pmModelType type);
+
+bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+
+bool psphotInitLimitsFLT (psMetadata *config, psStats *sky);
+
+bool psphotFitPSF (eamReadout *imdata, pmSource *source);
+bool psphotFitBlend (eamReadout *imdata, pmSource *source);
+bool psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources);
+
+pmModel *psphotFitFLT (eamReadout *imdata, pmSource *source);
+psArray *psphotFitDBL (eamReadout *imdata, pmSource *source);
+
+psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y);
+psPlane psImageBicubeMin (psPolynomial2D *poly);
Index: /trunk/psphot/src/psphotApResid.c
===================================================================
--- /trunk/psphot/src/psphotApResid.c	(revision 5827)
+++ /trunk/psphot/src/psphotApResid.c	(revision 5828)
@@ -21,5 +21,5 @@
     psVector *rflux   = psVectorAlloc (300, PS_TYPE_F64);
     psVector *apResid = psVectorAlloc (300, PS_TYPE_F64);
-    mask->n = rflux->n = apResid->n = 0;
+    mask->n = xPos->n = yPos->n = rflux->n = apResid->n = 0;
     Npsf = 0;
 
@@ -43,7 +43,8 @@
 
 	psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
-	// XXX EAM : add in source flux
 	status = pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);
 	psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
+
+	pmSourceSubModel (source->pixels, source->mask, model, false, false);
 	if (!status) continue;
 
@@ -55,11 +56,11 @@
 
 	psVectorExtend (mask, 100, 1);
+	psVectorExtend (xPos, 100, 1);
+	psVectorExtend (yPos, 100, 1);
 	psVectorExtend (rflux, 100, 1);
 	psVectorExtend (apResid, 100, 1);
 	Npsf ++;
-
-	pmSourceSubModel (source->pixels, source->mask, model, false, false);
-	psMemCheckCorruption (true);
     }
+    psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec\n", psTimerMark ("psphot"));
 
     // 3hi/1lo sigma clipping on the rflux vs metric fit
@@ -70,10 +71,29 @@
 
     // linear clipped fit of apResid to rflux, xPos, yPos
-    psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, rflux, xPos, yPos);
+    psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux);
     psf->skyBias  = psf->ApTrend->coeff[0][0][1] / (M_PI * PS_SQR(RADIUS));
     psf->ApResid  = psf->ApTrend->coeff[0][0][0];
     psf->dApResid = stats->sampleStdev;
+    psf->ApTrend->coeff[0][0][1] = 0;
 
-    psLogMsg ("ApResid", 4, "aperture residual: %f +/- %f : %f bias\n", psf->ApResid, psf->dApResid, psf->skyBias);
+    /*
+      (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
+      (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
+      (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
+    */
+
+    # if (0)
+    psPolynomial3D *poly = psf->ApTrend;
+    for (int nz = 0; nz <= poly->nZ; nz++) {
+	for (int ny = 0; ny <= poly->nY; ny++) {
+	    for (int nx = 0; nx <= poly->nX; nx++) {
+		fprintf (stderr, "%d %d %d : %22.15g\n", nx, ny, nz, poly->coeff[nx][ny][nz]);
+	    }
+	}
+    }
+    # endif
+
+    psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias\n", psf->ApResid, psf->dApResid, psf->skyBias);
 
     psFree (stats);
Index: /trunk/psphot/src/psphotApplyPSF.c
===================================================================
--- /trunk/psphot/src/psphotApplyPSF.c	(revision 5827)
+++ /trunk/psphot/src/psphotApplyPSF.c	(revision 5828)
@@ -34,5 +34,5 @@
 
 	// sets the model radius (via source->model) and adjusts pixels as needed
-	psphotCheckRadiusPSF (imdata, source);
+	psphotCheckRadiusPSF (imdata, source, model);
 
 	x = model->params->data.F32[2];
@@ -48,7 +48,8 @@
 
 	// check if model fit is acceptable
-	if (psphotEvalPSF (source)) {
+	if (psphotEvalPSF (source, source->modelPSF)) {
 	    pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
-	    source->mode |= PM_SOURCE_SUBTRACTED;
+	    source->mode |=  PM_SOURCE_SUBTRACTED;
+	    source->mode &= ~PM_SOURCE_TEMPSUB;
 	    Nsub ++;
 	}
Index: /trunk/psphot/src/psphotArguments.c
===================================================================
--- /trunk/psphot/src/psphotArguments.c	(revision 5827)
+++ /trunk/psphot/src/psphotArguments.c	(revision 5828)
@@ -5,72 +5,143 @@
 psMetadata *psphotArguments (int *argc, char **argv) {
 
-  int N;
-  unsigned int Nfail;
-  int mode = PS_DATA_STRING | PS_META_REPLACE;
+    int N;
+    unsigned int Nfail;
+    int mode = PS_DATA_STRING | PS_META_REPLACE;
 
-  // basic pslib options
-  fprintf (stderr, "starting... %s\n", psLibVersion());
-  psLogSetFormat ("M");
-  psArgumentVerbosity (argc, argv);
+    // basic pslib options
+    psLogSetFormat ("M");
+    psArgumentVerbosity (argc, argv);
 
-  // optional mask image - add to config
-  char *mask = NULL;
-  if ((N = psArgumentGet (*argc, argv, "-mask"))) {
-    psArgumentRemove (N, argc, argv);
-    mask = psStringCopy (argv[N]);
-    psArgumentRemove (N, argc, argv);
-  }
+    // optional mask image - add to config
+    char *mask = NULL;
+    if ((N = psArgumentGet (*argc, argv, "-mask"))) {
+	psArgumentRemove (N, argc, argv);
+	mask = psStringCopy (argv[N]);
+	psArgumentRemove (N, argc, argv);
+    }
 
-  // optional weight image - add to config
-  char *weight = NULL;
-  if ((N = psArgumentGet (*argc, argv, "-weight"))) {
-    psArgumentRemove (N, argc, argv);
-    weight = psStringCopy (argv[N]);
-    psArgumentRemove (N, argc, argv);
-  }
+    // optional weight image - add to config
+    char *weight = NULL;
+    if ((N = psArgumentGet (*argc, argv, "-weight"))) {
+	psArgumentRemove (N, argc, argv);
+	weight = psStringCopy (argv[N]);
+	psArgumentRemove (N, argc, argv);
+    }
 
-  // optional output residual image - add to config
-  char *resid = NULL;
-  if ((N = psArgumentGet (*argc, argv, "-resid"))) {
-    psArgumentRemove (N, argc, argv);
-    resid = psStringCopy (argv[N]);
-    psArgumentRemove (N, argc, argv);
-  }
+    // optional output residual image - add to config
+    char *resid = NULL;
+    if ((N = psArgumentGet (*argc, argv, "-resid"))) {
+	psArgumentRemove (N, argc, argv);
+	resid = psStringCopy (argv[N]);
+	psArgumentRemove (N, argc, argv);
+    }
 
-  // optional output residual image - add to config
-  char *photcode = NULL;
-  if ((N = psArgumentGet (*argc, argv, "-photcode"))) {
-    psArgumentRemove (N, argc, argv);
-    photcode = psStringCopy (argv[N]);
-    psArgumentRemove (N, argc, argv);
-  }
+    // optional analysis region - add to config
+    char *region = NULL;
+    if ((N = psArgumentGet (*argc, argv, "-region"))) {
+	psArgumentRemove (N, argc, argv);
+	region = psStringCopy (argv[N]);
+	psArgumentRemove (N, argc, argv);
+    }
 
-  if (*argc != 4) usage ();
+    // optional output residual image - add to config
+    char *photcode = NULL;
+    if ((N = psArgumentGet (*argc, argv, "-photcode"))) {
+	psArgumentRemove (N, argc, argv);
+	photcode = psStringCopy (argv[N]);
+	psArgumentRemove (N, argc, argv);
+    }
 
-  // load config information
-  psMetadata *config = psMetadataAlloc ();
-  psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL);
-  config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE);
+    char *psffile = NULL;
+    if ((N = psArgumentGet (*argc, argv, "-psf"))) {
+	psArgumentRemove (N, argc, argv);
+	psffile = psStringCopy (argv[N]);
+	psArgumentRemove (N, argc, argv);
+    }
 
-  // identify input image & optional weight & mask images
-  // command-line entries override config-file entries
-  psMetadataAdd (config, PS_LIST_HEAD, "IMAGE",       mode, "", argv[1]);
-  psMetadataAdd (config, PS_LIST_HEAD, "OUTPUT_FILE", mode, "", argv[2]);
+    bool ModelTest = false;
+    float ModelTest_X, ModelTest_Y;
+    char *model = NULL;
+    char *fitset = NULL;
+    char *fitmode = NULL;
+    if ((N = psArgumentGet (*argc, argv, "-modeltest"))) {
+	ModelTest = true;
+	psArgumentRemove (N, argc, argv);
+	ModelTest_X = atof (argv[N]);
+	psArgumentRemove (N, argc, argv);
+	ModelTest_Y = atof (argv[N]);
+	psArgumentRemove (N, argc, argv);
 
-  if (mask != NULL) {
-    psMetadataAdd (config, PS_LIST_HEAD, "MASK_IMAGE", mode, "", mask);
-  }
-  if (weight != NULL) {
-    psMetadataAdd (config, PS_LIST_HEAD, "WEIGHT_IMAGE", mode, "", weight);
-  }
-  if (resid != NULL) {
-    psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid);
-  }
-  if (photcode != NULL) {
-    psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", photcode);
-  } else {
-    psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", "NONE");
-  }    
-  return (config);
+	if ((N = psArgumentGet (*argc, argv, "-model"))) {
+	    psArgumentRemove (N, argc, argv);
+	    model = psStringCopy (argv[N]);
+	    psArgumentRemove (N, argc, argv);
+	}
+
+	if ((N = psArgumentGet (*argc, argv, "-fitmode"))) {
+	    psArgumentRemove (N, argc, argv);
+	    fitmode = psStringCopy (argv[N]);
+	    psArgumentRemove (N, argc, argv);
+	}
+	if ((N = psArgumentGet (*argc, argv, "-fitset"))) {
+	    psArgumentRemove (N, argc, argv);
+	    fitset = psStringCopy (argv[N]);
+	    psArgumentRemove (N, argc, argv);
+	}
+    }
+
+    if (*argc != 4) usage ();
+
+    // load config information
+    psMetadata *config = psMetadataAlloc ();
+    psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL);
+    config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE);
+
+    // identify input image & optional weight & mask images
+    // command-line entries override config-file entries
+    psMetadataAdd (config, PS_LIST_HEAD, "IMAGE",       mode, "", argv[1]);
+    psMetadataAdd (config, PS_LIST_HEAD, "OUTPUT_FILE", mode, "", argv[2]);
+
+    if (mask != NULL) {
+	psMetadataAdd (config, PS_LIST_HEAD, "MASK_IMAGE", mode, "", mask);
+    }
+    if (weight != NULL) {
+	psMetadataAdd (config, PS_LIST_HEAD, "WEIGHT_IMAGE", mode, "", weight);
+    }
+    if (resid != NULL) {
+	psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid);
+    }
+    if (region != NULL) {
+	psMetadataAdd (config, PS_LIST_HEAD, "ANALYSIS_REGION", mode, "", region);
+    }
+    if (photcode != NULL) {
+	psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", photcode);
+    } else {
+	psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", "NONE");
+    }    
+    if (psffile != NULL) {
+	psMetadataAdd (config, PS_LIST_HEAD, "PSF_INPUT_FILE", mode, "", psffile);
+    }
+
+    // model related options
+    if (ModelTest) {
+	psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT", PS_DATA_BOOL, "", true);
+
+	int fmode = PS_DATA_F32 | PS_META_REPLACE;
+	psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_X", fmode, "", ModelTest_X);
+	psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_Y", fmode, "", ModelTest_Y);
+
+	if (model != NULL) {
+	    psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODEL", mode, "", model);
+	}
+
+	if (fitmode != NULL) {
+	    psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODE", mode, "", fitmode);
+	}
+	if (fitset != NULL) {
+	    psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_SET", mode, "", fitset);
+	}
+    }
+    return (config);
 }
 
Index: /trunk/psphot/src/psphotBasicDeblend.c
===================================================================
--- /trunk/psphot/src/psphotBasicDeblend.c	(revision 5827)
+++ /trunk/psphot/src/psphotBasicDeblend.c	(revision 5828)
@@ -8,6 +8,4 @@
     pmSource *source, *testSource;
 
-    FILE *f = fopen ("deblend.dat", "w");
-    if (f == NULL) psAbort ("psphot", "can't open deblend.dat output file");
     int Nblend = 0;
 
@@ -82,8 +80,5 @@
 	    // XXX EAM : should the contour input coordinate be in parent or subimage coords? parent, for now
 	    psArray *contour = pmSourceContour_EAM (source->pixels, source->peak->x, source->peak->y, threshold);
-	    if (contour == NULL) {
-		fprintf (stderr, "below threshold? invalid peak?\n");
-		continue;
-	    }
+	    if (contour == NULL) continue;
 
 	    // the source contour consists of two vectors, xv and yv.  the contour is 
@@ -102,4 +97,5 @@
 		    if (xv->data.F32[j+1] < testSource->peak->x) break;
 
+		    # if (0)
 		    int xp0 = source->moments->x - source->pixels->col0;
 		    int xp1 = source->peak->x - source->pixels->col0;
@@ -111,13 +107,22 @@
 		    int yp2 = testSource->moments->y - testSource->pixels->row0;
 		    int yp3 = testSource->peak->y - testSource->pixels->row0;
+		    
+		    fprintf (f, "%d %d (%f, %f) :  %d %d (%f, %f)  vs %f\n",
+			     source->peak->x, source->peak->y,
+			     source->pixels->data.F32[yp0][xp0], source->pixels->data.F32[yp1][xp1],
+			     testSource->peak->x, testSource->peak->y,
+			     testSource->pixels->data.F32[yp2][xp2], testSource->pixels->data.F32[yp3][xp3], threshold
+		       );
+		    # endif
+		    
+		    testSource->mode |= PM_SOURCE_BLEND;
 
-		    fprintf (f, "%d %d (%f, %f) :  %d %d (%f, %f)  vs %f\n", 
-			     source->peak->x, source->peak->y, 
-			     source->pixels->data.F32[yp0][xp0], source->pixels->data.F32[yp1][xp1], 
-			     testSource->peak->x, testSource->peak->y, 
-			     testSource->pixels->data.F32[yp2][xp2], testSource->pixels->data.F32[yp3][xp3], threshold
-			     );
+		    // add this to the list of source->blends
+		    if (source->blends == NULL) {
+			source->blends = psArrayAlloc (16);
+			source->blends->n = 0;
+		    }
+		    psArrayAdd (source->blends, 16, testSource);
 
-		    testSource->mode |= PM_SOURCE_BLEND;
 		    Nblend ++;
 		    j = xv->n;
@@ -126,9 +131,5 @@
 	}
     }
-    fclose (f);
-    psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f)\n", Nblend, psTimerMark ("psphot"));
+    psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f sec)\n", Nblend, psTimerMark ("psphot"));
     return true;
 }
-
-
-
Index: /trunk/psphot/src/psphotBlendFit.c
===================================================================
--- /trunk/psphot/src/psphotBlendFit.c	(revision 5828)
+++ /trunk/psphot/src/psphotBlendFit.c	(revision 5828)
@@ -0,0 +1,45 @@
+# include "psphot.h"
+
+bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
+{ 
+
+    psTimerStart ("psphot");
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, psphotSortBySN);
+    
+    psphotInitLimitsPSF (config);
+    psphotInitLimitsFLT (config, sky);
+    psphotInitRadiusPSF (config, sky, psf->type);
+
+    for (int i = 0; i < sources->n; i++) {
+
+	pmSource *source = sources->data[i];
+
+	// skip non-astronomical objects (very likely defects)
+	if (source->mode &  PM_SOURCE_BLEND) continue;
+	if (source->type == PM_SOURCE_DEFECT) continue; 
+	if (source->type == PM_SOURCE_SATURATED) continue;
+
+	// if model is NULL, we don't have a starting guess
+	if (source->modelPSF == NULL) continue;
+
+	// replace object in image
+	pmSourceAddModel (source->pixels, source->mask, source->modelPSF, false, false);
+
+	psTrace ("psphot.blend", 5, "trying source at %f, %f\n", source->moments->x, source->moments->y);
+
+	if (psphotFitBlend (imdata, source)) continue;
+	if (psphotFitBlob (imdata, source, sources)) continue;
+
+	psTrace ("psphot.blend", 5, "failed fits\n");
+
+	// re-subtract PSF for object, leave local sky
+	pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
+	source->mode |= PM_SOURCE_SUBTRACTED;
+	source->mode |= PM_SOURCE_TEMPSUB;
+    }
+
+    psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects\n", psTimerMark ("psphot"), sources->n);
+    return (true);
+}
Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 5827)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 5828)
@@ -10,4 +10,6 @@
     psArray    	   *stars = NULL;
     psMetadataItem *item  = NULL;
+
+    psTimerStart ("psphot");
 
     // array to store candidate PSF stars
@@ -23,5 +25,5 @@
 	if (source->mode & PM_SOURCE_PSFSTAR) psArrayAdd (stars, 200, source);
     }
-    psLogMsg ("psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n);
+    psLogMsg ("psphot.pspsf", 4, "selected candidate %d PSF objects\n", stars->n);
 
     // get the fixed PSF fit radius
@@ -83,4 +85,5 @@
 
     modelName = pmModelGetType (psf->type);
+    psLogMsg ("psphot.pspsf", 3, "select psf model: %f sec\n", psTimerMark ("psphot"));
     psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
 
Index: /trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- /trunk/psphot/src/psphotEnsemblePSF.c	(revision 5827)
+++ /trunk/psphot/src/psphotEnsemblePSF.c	(revision 5828)
@@ -1,3 +1,49 @@
 # include "psphot.h"
+
+psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) {
+
+    int ix = x - image->col0;
+    int iy = y - image->row0;
+
+    psF32 *Fm = &image->data.F32[iy - 1][ix];
+    psF32 *Fo = &image->data.F32[iy + 0][ix];
+    psF32 *Fp = &image->data.F32[iy + 1][ix];
+
+    double Fxm = Fm[-1] + Fo[-1] + Fp[-1];
+    double Fxp = Fm[+1] + Fo[+1] + Fp[+1];
+    double Fym = Fm[-1] + Fm[+0] + Fm[+1];
+    double Fyp = Fp[-1] + Fp[+0] + Fp[+1];
+    double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1];
+
+    psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD);
+    poly->mask[2][2] = 1;
+    poly->mask[1][2] = 1;
+    poly->mask[2][1] = 1;
+
+    poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ;
+
+    poly->coeff[1][0] = (Fxp - Fxm)/6.0;
+    poly->coeff[0][1] = (Fyp - Fym)/6.0;
+    
+    poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0;
+    poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0;
+    
+    poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0;
+    
+    return (poly);
+}
+
+psPlane psImageBicubeMin (psPolynomial2D *poly) {
+
+    psPlane min;
+
+    min.xErr = min.yErr = 0;
+
+    double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]);
+
+    min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det;
+    min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det;
+    return (min);
+}
 
 bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
@@ -27,4 +73,15 @@
     index->n = 0;
 
+    bool UseAnalysisRegion = false;
+    psRegion AnalysisRegion;
+    char *region = psMetadataLookupPtr (&status, config, "ANALYSIS_REGION");
+    if (status) {
+	UseAnalysisRegion = true;
+	AnalysisRegion = psRegionFromString (region);
+	psLogMsg ("psphotEnsemblePSF", 4, "using region %f,%f - %f,%f\n", 
+		  AnalysisRegion.x0, AnalysisRegion.y0, 
+		  AnalysisRegion.x1, AnalysisRegion.y1);
+    }
+
     for (int i = 0; i < sources->n; i++) {
 	pmSource *inSource = sources->data[i];
@@ -35,4 +92,11 @@
 	if (inSource->type == PM_SOURCE_DEFECT) continue; 
 	if (inSource->type == PM_SOURCE_SATURATED) continue;
+
+	if (UseAnalysisRegion) {
+	    if (inSource->moments->x < AnalysisRegion.x0) continue;
+	    if (inSource->moments->y < AnalysisRegion.y0) continue;
+	    if (inSource->moments->x > AnalysisRegion.x1) continue;
+	    if (inSource->moments->y > AnalysisRegion.y1) continue;
+	}
 
 	pmSource *otSource = pmSourceAlloc ();
@@ -51,6 +115,23 @@
 	    modelFLT->params->data.F32[2] = inSource->moments->x;
 	    modelFLT->params->data.F32[3] = inSource->moments->y;
-	}
-	// XXX EAM : add option to peak-up on peak (for non-sat objects)
+	} else {
+	    // peak-up on peak (for non-sat objects)
+
+	    int ix = inSource->peak->x;
+	    int iy = inSource->peak->y;
+
+	    psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
+	    psPlane min = psImageBicubeMin (bicube);
+
+	    psTrace ("psphotEnsemblePSF", 5, "peak coord: %f %f -> %f %f\n", 
+		     modelFLT->params->data.F32[2], modelFLT->params->data.F32[3], min.x + ix, min.y + iy);
+	    
+	    // if min point is too deviant, keep the old value
+	    if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
+		modelFLT->params->data.F32[2] = min.x + ix;
+		modelFLT->params->data.F32[3] = min.y + iy;
+	    }
+	    psFree (bicube);
+	}
 
 	// set PSF parameters for this model
@@ -149,9 +230,11 @@
 	// subtract object
 	pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
+	Fi->mode |= PM_SOURCE_SUBTRACTED;
+	Fi->mode |= PM_SOURCE_TEMPSUB;
     }
 
     // XXX EAM : need to free up many things here
 
-    psLogMsg ("psphot.emsemble", 4, "apply models: %f\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f\n", psTimerMark ("psphot"));
     return true;
 }
Index: /trunk/psphot/src/psphotEvalFLT.c
===================================================================
--- /trunk/psphot/src/psphotEvalFLT.c	(revision 5827)
+++ /trunk/psphot/src/psphotEvalFLT.c	(revision 5828)
@@ -1,9 +1,7 @@
 # include "psphot.h"
 
-bool psphotEvalFLT (pmSource *source)
+bool psphotEvalFLT (pmSource *source, pmModel *model)
 { 
     int keep;
-
-    pmModel *model = source->modelFLT;
 
     // do we actually have a valid FLT model?
Index: /trunk/psphot/src/psphotEvalPSF.c
===================================================================
--- /trunk/psphot/src/psphotEvalPSF.c	(revision 5827)
+++ /trunk/psphot/src/psphotEvalPSF.c	(revision 5828)
@@ -46,11 +46,9 @@
 // examine the model->status, fit parameters, etc and decide if the model succeeded
 // set the source->type and source->mode appropriately
-bool psphotEvalPSF (pmSource *source) { 
+bool psphotEvalPSF (pmSource *source, pmModel *model) { 
 
     int keep;
     float dSX, dSY, SX, SY, SN;
     float nSx, nSy, Chi;
-
-    pmModel *model = source->modelPSF;
 
     // do we actually have a valid PSF model?
@@ -127,8 +125,8 @@
     // this source is not a star, warn if it was a PSFSTAR
     if (source->mode & PM_SOURCE_PSFSTAR) {
-	psphotSaveImage (NULL, source->pixels, "failpx.fits");
-	psphotSaveImage (NULL, source->mask, "failmk.fits");
-	psphotSaveImage (NULL, source->weight, "failwt.fits");
 	psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality   (%f, %f  :  %f %f %f %f)\n", 
+		  model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi);
+    } else {
+	psLogMsg ("psphot", 5, "fails PSF fit (%f, %f  :  %f %f %f %f)\n", 
 		  model->params->data.F32[2], model->params->data.F32[3], nSx, nSy, SN, Chi);
     }
@@ -150,2 +148,60 @@
     return false;
 }	
+
+// examine the model->status, fit parameters, etc and decide if the model succeeded
+// set the source->type and source->mode appropriately
+bool psphotEvalDBL (pmSource *source, pmModel *model) { 
+
+    // do we actually have a valid PSF model?
+    if (model == NULL) {
+	source->mode &= ~PM_SOURCE_FITTED;
+	return false;
+    }
+
+    // did the model fit fail for one or another reason?
+    switch (model->status) {
+      case PM_MODEL_SUCCESS:
+	break;
+      case PM_MODEL_UNTRIED:
+	source->mode &= ~PM_SOURCE_FITTED; 
+	return false;
+      case PM_MODEL_BADARGS:
+      case PM_MODEL_NONCONVERGE:
+      case PM_MODEL_OFFIMAGE:
+      default:
+	source->mode |= PM_SOURCE_FAIL;
+	return false;
+    }
+
+    // unless we prove otherwise, this object is a star.
+    source->type = PM_SOURCE_STAR;
+
+    // the following source->mode information pertains to modelPSF:
+    source->mode |= PM_SOURCE_PSFMODEL;
+
+    // if the object has fitted peak above saturation, label as SATSTAR
+    // this is a valid PSF object, but ignore the other quality tests
+    // remember: fit does not use saturated pixels (masked)
+    // XXX no extended object can saturate and stay extended...
+    if (model->params->data.F32[1] >= SATURATION) {
+	if (source->mode & PM_SOURCE_PSFSTAR) {
+	    psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n");
+	}
+	source->mode |=  PM_SOURCE_SATSTAR;
+	return true;
+    } 
+
+    // if the object has a fitted peak below 0, the fit did not converge cleanly
+    if (model->params->data.F32[1] <= 0) {
+	source->mode |= PM_SOURCE_FAIL;
+	return false;
+    } 
+
+    // if the source was predicted to be a SATSTAR, but it fitted below saturation, 
+    // make a note to the user
+    if (source->mode & PM_SOURCE_SATSTAR) {
+	psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n");
+	source->mode &= ~PM_SOURCE_SATSTAR;
+    }
+    return true;
+}	
Index: /trunk/psphot/src/psphotFitGalaxies.c
===================================================================
--- /trunk/psphot/src/psphotFitGalaxies.c	(revision 5827)
+++ /trunk/psphot/src/psphotFitGalaxies.c	(revision 5828)
@@ -43,5 +43,5 @@
 
 	// sets the model radius (via source->model) and adjusts pixels as needed
-	psphotCheckRadiusFLT (imdata, source);
+	psphotCheckRadiusFLT (imdata, source, model);
 
 	// fit FLT (not PSF) model (set/unset the pixel mask)
@@ -56,5 +56,6 @@
 	if (pmModelFitStatus (model)) {
 	    pmSourceSubModel (source->pixels, source->mask, model, false, false);
-	    source->mode |= PM_SOURCE_SUBTRACTED;
+	    source->mode |=  PM_SOURCE_SUBTRACTED;
+	    source->mode &= ~PM_SOURCE_TEMPSUB;
 	    Nsub ++;
 	}
Index: /trunk/psphot/src/psphotFixedPSF.c
===================================================================
--- /trunk/psphot/src/psphotFixedPSF.c	(revision 5827)
+++ /trunk/psphot/src/psphotFixedPSF.c	(revision 5828)
@@ -77,4 +77,6 @@
 	// subtract object, leave local sky
 	pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
+	source->mode |=  PM_SOURCE_SUBTRACTED;
+	source->mode &= ~PM_SOURCE_TEMPSUB;
     }
 
Index: /trunk/psphot/src/psphotFullFit.c
===================================================================
--- /trunk/psphot/src/psphotFullFit.c	(revision 5827)
+++ /trunk/psphot/src/psphotFullFit.c	(revision 5828)
@@ -1,9 +1,3 @@
 # include "psphot.h"
-
-# define CLEAR_TRACE \
-  if (TEST_SRC_ON) { \
-    psTraceSetLevel (".pmObjects.pmSourceFitModel", 0); \
-    psTraceSetLevel ("psMinimizeLMChi2", 0); \
-    TEST_SRC_ON = false; } \
 
 bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
@@ -32,10 +26,4 @@
     psphotInitRadiusFLT (config, sky, modelTypeFLT);
 
-    float TEST_SRC_X = psMetadataLookupF32 (&status, config, "TEST_SRC_X");
-    float TEST_SRC_Y = psMetadataLookupF32 (&status, config, "TEST_SRC_Y");
-    bool TEST_SRC_ON = false;
-    bool TEST_SRC    = status;
-    fprintf (stderr, "test src: %f %f\n", TEST_SRC_X, TEST_SRC_Y);
-
     for (int i = 0; i < sources->n; i++) {
 
@@ -56,15 +44,10 @@
 	// replace object in image
 	pmSourceAddModel (source->pixels, source->mask, modelPSF, false, false);
+	source->mode &= ~PM_SOURCE_SUBTRACTED;
 
-	psphotCheckRadiusPSF (imdata, source);
+	psphotCheckRadiusPSF (imdata, source, modelPSF);
 
 	x = modelPSF->params->data.F32[2];
 	y = modelPSF->params->data.F32[3];
-
-	if (TEST_SRC && (fabs(TEST_SRC_X - x) < 5) && (fabs(TEST_SRC_Y - y) < 5)) {
-	    psTraceSetLevel (".pmObjects.pmSourceFitModel", 6);
-	    psTraceSetLevel ("psMinimizeLMChi2", 6);
-	    TEST_SRC_ON = true;
-	}
 
 	// fit PSF model (set/unset the pixel mask)
@@ -78,12 +61,23 @@
 
 	// does the PSF model succeed?
-	if (psphotEvalPSF (source)) {
+	if (psphotEvalPSF (source, source->modelPSF)) {
 	    pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
-	    source->mode |= PM_SOURCE_SUBTRACTED;
+	    source->mode |=  PM_SOURCE_SUBTRACTED;
+	    source->mode &= ~PM_SOURCE_TEMPSUB;
 	    psFree (LIN);
 	    Nsub ++;
-	    CLEAR_TRACE;
 	    continue;
 	} 
+
+	// use the moments to get the extended source angle
+	// displace on either side of the main peak by 1sigma or so
+	// cut the peak in half?  use the main peak?
+	// also try fit set on (unsaturated?) objects with a blend within 2-3 sigma
+	// allow all possible nearby blends to be fit?
+	// EllipseMoments moments;
+	// moments.x2 = source->moments->Sx;
+	// moments.y2 = source->moments->Sy;
+	// moments.xy = source->moments->Sxy;
+	// EllipseAxes axes = EllipseMomentsToAxes (moments);
 
 	// skip the source if we don't think it is extended
@@ -101,5 +95,5 @@
 	pmModel *modelFLT = source->modelFLT;
 
-	psphotCheckRadiusFLT (imdata, source);
+	psphotCheckRadiusFLT (imdata, source, modelFLT);
 
 	x = modelFLT->params->data.F32[2];
@@ -116,14 +110,14 @@
 
 	// does the FLT model succeed?
-	if (psphotEvalFLT (source)) {
+	if (psphotEvalFLT (source, source->modelFLT)) {
 	    pmSourceSubModel (source->pixels, source->mask, source->modelFLT, false, false);
-	    source->mode |= PM_SOURCE_SUBTRACTED;
+	    source->mode |=  PM_SOURCE_SUBTRACTED;
+	    source->mode &= ~PM_SOURCE_TEMPSUB;
 	    psFree (LIN);
 	    Nsub ++;
-	    CLEAR_TRACE;
 	    continue;
 	}
 
-	// subtract PSF for object, leave local sky
+	// re-subtract PSF for object, leave local sky
     subLINEAR:
 	psFree (source->modelPSF);
@@ -131,7 +125,6 @@
 	pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false);
 	source->mode |= PM_SOURCE_SUBTRACTED;
-	source->mode |= PM_SOURCE_LINEAR;
+	source->mode |= PM_SOURCE_TEMPSUB;
 	Nsub ++;
-	CLEAR_TRACE;
     }
 
Index: /trunk/psphot/src/psphotImageBackground.c
===================================================================
--- /trunk/psphot/src/psphotImageBackground.c	(revision 5827)
+++ /trunk/psphot/src/psphotImageBackground.c	(revision 5828)
@@ -38,5 +38,4 @@
     }
     x->n = y->n = z->n = Nout;
-    psLogMsg ("psphot", 5, "back: %f sec (select %d points)\n", psTimerMark ("psphot"), Nout);
 
     psPolynomial2D *skyModel = psPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
@@ -49,8 +48,4 @@
 
     skyModel = psVectorFitPolynomial2D (skyModel, NULL, 0, z, NULL, x, y);
-    psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot"));
-    fprintf (stderr, "model: %f %f %f %f\n", 
-	     skyModel->coeff[0][0], skyModel->coeff[1][0], 
-	     skyModel->coeff[0][1], skyModel->coeff[1][1]);
 	     
     // this is a very inefficient way to evaluate the function..
@@ -66,7 +61,5 @@
     sky->sampleMean = 0.0;
 
-    psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot"));
-
-    psphotSaveImage (imdata->header, imdata->image, "backsub.fits");
+    psLogMsg ("psphot", 3, "fit background model: %f sec\n", psTimerMark ("psphot"));
     return (skyModel);
 }
Index: /trunk/psphot/src/psphotImageStats.c
===================================================================
--- /trunk/psphot/src/psphotImageStats.c	(revision 5827)
+++ /trunk/psphot/src/psphotImageStats.c	(revision 5828)
@@ -49,6 +49,6 @@
     sky->sampleStdev  = sqrt(sky->sampleMean/GAIN + PS_SQR(NOISE));
 
-    psLogMsg ("psphot", 4, "stats: %f sec\n", psTimerMark ("psphot"));
-    psLogMsg ("psphot", 3, "background: %f +/- %f\n", sky->sampleMean, sky->sampleStdev);
+    psLogMsg ("psphot", 4, "background: %f +/- %f\n", sky->sampleMean, sky->sampleStdev);
+    psLogMsg ("psphot", 3, "image stats: %f sec\n", psTimerMark ("psphot"));
 
     psFree (stats);
Index: /trunk/psphot/src/psphotMagnitudes.c
===================================================================
--- /trunk/psphot/src/psphotMagnitudes.c	(revision 5827)
+++ /trunk/psphot/src/psphotMagnitudes.c	(revision 5828)
@@ -37,5 +37,5 @@
       rflux   = pow (10.0, 0.4*source->fitMag);
       source->apMag  -= rflux * psf->skyBias * (M_PI * PS_SQR(apRadius));
-      source->fitMag += psf->ApResid;
+      source->fitMag += psPolynomial3DEval (psf->ApTrend, x, y, 0.0);
     }
 
@@ -70,2 +70,8 @@
     return (model);
 }
+
+/*
+aprMag' - fitMag = rflux*skyBias + ApTrend(x,y)
+(aprMag - rflux*skyBias) - fitMAg = ApTrend(x,y)
+(aprMag - rflux*skyBias) = fitMAg + ApTrend(x,y)
+*/
Index: /trunk/psphot/src/psphotModelTest.c
===================================================================
--- /trunk/psphot/src/psphotModelTest.c	(revision 5827)
+++ /trunk/psphot/src/psphotModelTest.c	(revision 5828)
@@ -1,13 +1,169 @@
 # include "psphot.h"
+# include "psEllipse.h"
 
-int main (int argc, char **argv) {
+bool psphotModelTest (eamReadout *imdata, psMetadata *config) {
 
-    psphotModelGroupInit ();
+    bool status;
+    int modelType;
+    float obsMag, fitMag, value;
+    char name[64];
+    pmPSF *psf;
 
-    psMetadata *config = modelTestArguments (&argc, argv);
-    eamReadout *imdata = psphotSetup (config);
+    psMetadataItem *item  = NULL;
 
-    modelTestFitSource (imdata, config);
+    // what fitting mode to use?
+    char *psfModeWord = psMetadataLookupPtr (&status, config, "TEST_FIT_MODE");
+    if (!status) {
+	psfModeWord = psStringCopy ("FLT");
+    }
+    bool psfMode = !strcasecmp (psfModeWord, "PSF");
 
-    exit (0);
+    // in psfMode, psf sets the model type
+    if (psfMode) {
+	char *psfFile = psMetadataLookupPtr (&status, config, "PSF_INPUT_FILE");
+	if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied");
+	psf = psphotReadPSF (psfFile);
+	modelType = psf->type;
+    } else {
+	// find the model: supplied by user or first in the PSF_MODEL list
+	char *modelName  = psMetadataLookupPtr (&status, config, "TEST_FIT_MODEL");
+	if (modelName == NULL) {
+	    // get the list pointers for the PSF_MODEL entries
+	    psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL");
+	    if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");
+
+	    // take the first list element
+	    psList *list = (psList *) mdi->data.list;
+	    item = psListGet (list, PS_LIST_HEAD);
+	    modelName = item->data.V;
+	}
+	modelType = pmModelSetType (modelName);
+	if (modelType < 0) psAbort ("fitsource", "unknown model %s", modelName);
+    }
+
+    // find the fitting parameters (try test values first)
+    float INNER = psMetadataLookupF32 (&status, config, "TEST_FIT_INNER_RADIUS");
+    if (!status) {
+	INNER = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS");
+    }
+
+    float OUTER = psMetadataLookupF32 (&status, config, "TEST_FIT_OUTER_RADIUS");
+    if (!status) {
+	OUTER = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
+    }
+
+    float RADIUS = psMetadataLookupF32 (&status, config, "TEST_FIT_RADIUS");
+    if (!status) {
+	RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
+    }
+
+    float mRADIUS = psMetadataLookupF32 (&status, config, "TEST_MOMENTS_RADIUS");
+    if (!status) {
+	mRADIUS = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS");
+    }
+
+    // define the source of interest
+    float xObj     = psMetadataLookupF32 (&status, config, "TEST_FIT_X");
+    float yObj     = psMetadataLookupF32 (&status, config, "TEST_FIT_Y");
+
+    // construct the source structures
+    pmSource *source = pmSourceAlloc();
+    source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
+    psphotDefinePixels (source, imdata, xObj, yObj, OUTER);
+
+    // find the local sky
+    status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
+    if (!status) psAbort ("fitsource", "pmSourceLocalSky error");
+
+    // get the source moments
+    status = pmSourceMoments (source, mRADIUS);
+    if (!status) psAbort ("fitsource", "pmSourceLocalSky error");
+    source->peak->counts = source->moments->Peak;
+
+    fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->x, source->moments->y);
+
+    EllipseMoments moments;
+    moments.x2 = source->moments->Sx;
+    moments.y2 = source->moments->Sy;
+    moments.xy = source->moments->Sxy;
+    EllipseAxes axes = EllipseMomentsToAxes (moments);
+
+    fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
+
+    // get the initial model parameter guess
+    pmModel *model = pmSourceModelGuess (source, modelType);
+    // if any parameters are defined, use those values
+    int nParams = pmModelParameterCount (modelType);
+    psF32 *params = model->params->data.F32;
+    for (int i = 0; i < nParams; i++) {
+	if (i == 2) {
+	    params[i] = xObj;
+	    continue;
+	}
+	if (i == 3) {
+	    params[i] = yObj;
+	    continue;
+	}
+	sprintf (name, "TEST_FIT_PAR%d", i);
+	value = psMetadataLookupF32 (&status, config, name);
+	if (status) {
+	    params[i] = value;
+	}
+    }
+
+    float area = params[4]*params[5];
+    fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
+
+    if (psfMode) {
+	pmModel *modelPSF = pmModelFromPSF (model, psf);
+	psFree (model);
+	model = modelPSF;
+	params = model->params->data.F32;
+    }
+
+    // list model input shape
+    EllipseShape shape;
+    shape.sx  = 1.4 / model->params->data.F32[4];
+    shape.sy  = 1.4 / model->params->data.F32[5];
+    shape.sxy = model->params->data.F32[6];
+    axes = EllipseShapeToAxes (shape);
+
+    fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
+
+
+    fprintf (stderr, "input parameters: \n");
+    for (int i = 0; i < nParams; i++) {
+	fprintf (stderr, "%d : %f\n", i, params[i]);
+    }
+
+    // define the pixels used for the fit
+    psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PSPHOT_MASK_MARKED);
+
+    char *fitset = psMetadataLookupPtr (&status, config, "TEST_FIT_SET");
+    if (status) {
+	status = psphotFitSet (source, model, fitset, psfMode);
+	return status;
+    }
+
+    status = pmSourceFitModel (source, model, psfMode);
+
+    // measure the source mags
+    pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);
+    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag);
+
+    // write out positive object
+    psphotSaveImage (NULL, source->pixels, "object.fits");
+
+    // subtract object, leave local sky
+    pmSourceSubModel (source->pixels, source->mask, model, false, false);
+    
+    fprintf (stderr, "output parameters: \n");
+    for (int i = 0; i < nParams; i++) {
+	fprintf (stderr, "%d : %f\n", i, params[i]);
+    }
+
+    // write out 
+    psphotSaveImage (NULL, source->pixels, "resid.fits");
+    psphotSaveImage (NULL, source->mask, "mask.fits");
+    return true;
 }
Index: /trunk/psphot/src/psphotOutput.c
===================================================================
--- /trunk/psphot/src/psphotOutput.c	(revision 5827)
+++ /trunk/psphot/src/psphotOutput.c	(revision 5828)
@@ -5,10 +5,14 @@
 
     bool status;
+
+    psTimerStart ("psphot");
 
     char *outputMode = psMetadataLookupPtr (&status, config, "OUTPUT_MODE");
     char *outputFile = psMetadataLookupPtr (&status, config, "OUTPUT_FILE");
     char *residImage = psMetadataLookupPtr (&status, config, "RESID_IMAGE");
+    char *psfFile    = psMetadataLookupPtr (&status, config, "PSF_OUTPUT_FILE");
 
     if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage);
+    if (psfFile != NULL) psphotWritePSF (psf, psfFile);
 
     if (outputFile == NULL) {
@@ -44,4 +48,5 @@
 	return;
     }
+
     psAbort ("psphot", "unknown output mode %s", outputMode);
 }
@@ -287,5 +292,4 @@
 
     // set NAXIS to 0
-    fprintf (stderr, "setting naxis\n"); 
     mdi = psMetadataLookup (imdata->header, "NAXIS");
     mdi->data.S32 = 0;
@@ -540,5 +544,4 @@
     }
 
-    fprintf (stderr, "writing out moments\n");
     for (i = 0; i < sources->n; i++) {
         source = sources->data[i];
@@ -604,8 +607,16 @@
 
 
-bool psphotSamplePSFs (pmPSF *psf, psImage *image) {
-
-  // make sample PSFs for 4 corners and the center
+bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image) {
+
+    bool status;
+
+    // make sample PSFs for 4 corners and the center
     psImage *sample;
+
+    // optional dump of all rough source data
+    char *output = psMetadataLookupPtr (&status, config, "PSF_SAMPLE_FILE");
+    if (!status) return false;
+    if (output == NULL) return false;
+    if (output[0] == 0) return false;
 
     pmModel *modelFLT = pmModelAlloc (psf->type);
@@ -613,5 +624,5 @@
     modelFLT->params->data.F32[1] = 1;
 
-    psFits *fits = psFitsOpen ("sample.psf.fits", "w");
+    psFits *fits = psFitsOpen (output, "w");
 
     sample = pmModelPSFatXY (image, modelFLT, psf, 25, 25, 25, 25);
@@ -630,4 +641,206 @@
 }
 
+psPolynomial2D *psPolynomial2DfromMD (psMetadata *folder) {
+
+    bool status;
+    char keyword[80];
+
+    // get polynomial orders
+    // XXX add status failures tests
+    int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
+    int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
+
+    psPolynomial2D *poly = psPolynomial2DAlloc (nXorder, nYorder, PS_POLYNOMIAL_ORD);
+
+    for (int nx = 0; nx < poly->nX + 1; nx++) {
+	for (int ny = 0; ny < poly->nY + 1; ny++) {
+	    sprintf (keyword, "VAL_X%02d_Y%02d", nx, ny);
+	    poly->coeff[nx][ny] = psMetadataLookupF32 (&status, folder, keyword);
+	    if (!status) poly->mask[nx][ny] = 1;
+
+	    sprintf (keyword, "ERR_X%02d_Y%02d", nx, ny);
+	    poly->coeffErr[nx][ny] = psMetadataLookupF32 (&status, folder, keyword);
+	}
+    }
+    return (poly);
+}
+
+psPolynomial3D *psPolynomial3DfromMD (psMetadata *folder) {
+
+    bool status;
+    char keyword[80];
+
+    // get polynomial orders
+    // XXX add status failures tests
+    int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
+    int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
+    int nZorder = psMetadataLookupS32 (&status, folder, "NORDER_Z");
+
+    psPolynomial3D *poly = psPolynomial3DAlloc (nXorder, nYorder, nZorder, PS_POLYNOMIAL_ORD);
+
+    for (int nx = 0; nx < poly->nX + 1; nx++) {
+	for (int ny = 0; ny < poly->nY + 1; ny++) {
+	    for (int nz = 0; nz < poly->nZ + 1; nz++) {
+		sprintf (keyword, "VAL_X%02d_Y%02d_Z%02d", nx, ny, nz);
+		poly->coeff[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword);
+		if (!status) poly->mask[nx][ny][nz] = 1;
+
+		sprintf (keyword, "ERR_X%02d_Y%02d_Z%02d", nx, ny, nz);
+		poly->coeffErr[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword);
+	    }
+	}
+    }
+    return (poly);
+}
+
+// XXX : these may need F64, or %g format for output
+bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...) {
+
+    int Nbyte;
+    char tmp;
+    char *root;
+    va_list argp;  
+
+    va_start (argp, format);
+    Nbyte = vsnprintf (&tmp, 0, format, argp);
+    va_end (argp);
+
+    if (!Nbyte) return false;
+
+    va_start (argp, format);
+    root = (char *) psAlloc (Nbyte + 1);
+    memset (root, 0, Nbyte + 1);
+    vsnprintf (root, Nbyte + 1, format, argp);
+    va_end (argp);
+
+    psMetadata *folder = psMetadataAlloc ();
+    psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 2D polynomial", folder);
+    psFree (root);
+
+    // specify the polynomial orders
+    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
+    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
+
+    // place polynomial entries on folder
+    for (int nx = 0; nx < poly->nX + 1; nx++) {
+	for (int ny = 0; ny < poly->nY + 1; ny++) {
+	    if (poly->mask[nx][ny]) continue;
+	    psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny], nx, ny);
+	    psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient error", poly->coeffErr[nx][ny], nx, ny);
+	}
+    }
+    return true;
+}
+
+bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *poly, char *format, ...) {
+
+    int Nbyte;
+    char tmp;
+    char *root;
+    va_list argp;  
+
+    va_start (argp, format);
+    Nbyte = vsnprintf (&tmp, 0, format, argp);
+    va_end (argp);
+
+    if (!Nbyte) return false;
+
+    va_start (argp, format);
+    root = (char *) psAlloc (Nbyte + 1);
+    memset (root, 0, Nbyte + 1);
+    vsnprintf (root, Nbyte + 1, format, argp);
+    va_end (argp);
+
+    psMetadata *folder = psMetadataAlloc ();
+    psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 3D polynomial", folder);
+    psFree (root);
+
+    // specify the polynomial orders
+    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
+    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
+    psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Z", PS_DATA_S32, "number of z orders", poly->nZ);
+
+    // place polynomial entries on folder
+    for (int nx = 0; nx < poly->nX + 1; nx++) {
+	for (int ny = 0; ny < poly->nY + 1; ny++) {
+	    for (int nz = 0; nz < poly->nZ + 1; nz++) {
+		if (poly->mask[nx][ny][nz]) continue;
+		psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny][nz], nx, ny, nz);
+		psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coeffficient error", poly->coeffErr[nx][ny][nz], nx, ny, nz);
+	    }
+	}
+    }
+    return true;
+}
+
+bool psphotWritePSF (pmPSF *psf, char *filename) {
+
+    psMetadata *psfdata = psMetadataAlloc ();
+
+    char *modelName = pmModelGetType (psf->type);
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
+
+    int nPar = pmModelParameterCount (psf->type)    ;
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
+
+    for (int i = 0; i < nPar - 4; i++) {
+	psPolynomial2D *poly = psf->params->data[i];
+	psPolynomial2DtoMD (psfdata, poly, "PSF_PAR%02d", i);
+    }
+    psPolynomial3DtoMD (psfdata, psf->ApTrend, "APTREND");
+
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid);
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid);
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
+
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq);
+    psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars);
+    
+    psMetadataConfigWrite (psfdata, filename);
+    psFree (psfdata);
+
+    return true;
+}
+
+pmPSF *psphotReadPSF (char *filename) {
+
+    bool status;
+    unsigned int Nfail;
+    char keyword[80];
+
+    psMetadata *psfdata = psMetadataConfigParse (NULL, &Nfail, filename, FALSE);
+
+    char *modelName = psMetadataLookupPtr (&status, psfdata, "PSF_MODEL_NAME");
+    pmModelType type = pmModelSetType (modelName);
+
+    pmPSF *psf = pmPSFAlloc (type);
+
+    int nPar = psMetadataLookupS32 (&status, psfdata, "PSF_MODEL_NPAR");
+    if (nPar != pmModelParameterCount (psf->type)) psAbort ("read PSF" , "mismatch model par count");
+
+    for (int i = 0; i < nPar - 4; i++) {
+	sprintf (keyword, "PSF_PAR%02d", i);
+	psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
+	psPolynomial2D *poly = psPolynomial2DfromMD (folder);
+	psFree (psf->params->data[i]);
+	psf->params->data[i] = poly;
+    }
+    sprintf (keyword, "APTREND");
+    psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
+    psPolynomial3D *poly = psPolynomial3DfromMD (folder);
+    psFree (psf->ApTrend);
+    psf->ApTrend = poly;
+
+    psf->ApResid = psMetadataLookupF32 (&status, psfdata, "PSF_AP_RESID");
+    psf->dApResid = psMetadataLookupF32 (&status, psfdata, "PSF_dAP_RESID");
+    psf->skyBias = psMetadataLookupF32 (&status, psfdata, "PSF_SKY_BIAS");
+
+    psf->chisq = psMetadataLookupF32 (&status, psfdata, "PSF_CHISQ");
+    psf->nPSFstars = psMetadataLookupS32 (&status, psfdata, "PSF_NSTARS");
+
+    psFree (psfdata);
+    return (psf);
+}
+
 bool psphotDumpMoments (psMetadata *config, psArray *sources) {
 
Index: /trunk/psphot/src/psphotRadiusChecks.c
===================================================================
--- /trunk/psphot/src/psphotRadiusChecks.c	(revision 5827)
+++ /trunk/psphot/src/psphotRadiusChecks.c	(revision 5828)
@@ -20,7 +20,5 @@
 }
 
-bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source) {
-
-    pmModel *model = source->modelPSF;
+bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source, pmModel *model) {
 
     // set the fit radius based on the object flux limit and the model
@@ -64,7 +62,5 @@
 }
 
-bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source) {
-
-    pmModel *model = source->modelFLT;
+bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source, pmModel *model) {
 
     // set the fit radius based on the object flux limit and the model
Index: /trunk/psphot/src/psphotReplaceUnfit.c
===================================================================
--- /trunk/psphot/src/psphotReplaceUnfit.c	(revision 5827)
+++ /trunk/psphot/src/psphotReplaceUnfit.c	(revision 5828)
@@ -3,5 +3,4 @@
 bool psphotReplaceUnfit (psArray *sources) { 
 
-    pmModel *model;
     pmSource *source;
 
@@ -11,35 +10,35 @@
       source = sources->data[i];
 
-	// skip non-astronomical objects (very likely defects)
-	// these were never fitted and subtracted
-	if (source->mode &  PM_SOURCE_BLEND) continue;
-	if (source->type == PM_SOURCE_SATURATED) continue;
+	if (!(source->mode & PM_SOURCE_TEMPSUB)) continue;
+	if (source->modelPSF == NULL) continue;
 
-	if (source->type == PM_SOURCE_DEFECT) {
-	    if (source->mode &  PM_SOURCE_SUBTRACTED) goto addPSF;
-	    continue; 
-	}
-
-	if (source->type == PM_SOURCE_GALAXY) {
-	    if (source->mode & PM_SOURCE_FAIL) goto addPSF;
-	    if (source->mode & PM_SOURCE_POOR) goto addPSF;
-	    continue;
-	}
-
-	// need to skip the successful fits
-	if (source->type == PM_SOURCE_STAR) {
-	    if (source->mode & PM_SOURCE_FAIL) goto addPSF;
-	    if (source->mode & PM_SOURCE_POOR) goto addPSF;
-	    continue;
-	}
-
-    addPSF:
-	model = source->modelPSF;
-	if (model == NULL) continue;
-
-	pmSourceAddModel (source->pixels, source->mask, model, false, false);
+	pmSourceAddModel (source->pixels, source->mask, source->modelPSF, false, false);
 	source->mode &= ~PM_SOURCE_SUBTRACTED;
+	source->mode &= ~PM_SOURCE_TEMPSUB;
     }
-    psLogMsg ("psphot.replace", 4, "replace models: %f (%d objects)\n", psTimerMark ("psphot"), sources->n);
+    psLogMsg ("psphot.replace", 3, "replace unfitted models: %f sec (%d objects)\n", psTimerMark ("psphot"), sources->n);
     return true;
 }
+
+# if (0)
+if (source->mode &  PM_SOURCE_BLEND) continue;
+if (source->type == PM_SOURCE_SATURATED) continue;
+
+if (source->type == PM_SOURCE_DEFECT) {
+    if (source->mode &  PM_SOURCE_SUBTRACTED) goto addPSF;
+    continue; 
+}
+
+if (source->type == PM_SOURCE_GALAXY) {
+    if (source->mode & PM_SOURCE_FAIL) goto addPSF;
+    if (source->mode & PM_SOURCE_POOR) goto addPSF;
+    continue;
+}
+
+// need to skip the successful fits
+if (source->type == PM_SOURCE_STAR) {
+    if (source->mode & PM_SOURCE_FAIL) goto addPSF;
+    if (source->mode & PM_SOURCE_POOR) goto addPSF;
+    continue;
+}
+# endif
Index: /trunk/psphot/src/psphotSetup.c
===================================================================
--- /trunk/psphot/src/psphotSetup.c	(revision 5827)
+++ /trunk/psphot/src/psphotSetup.c	(revision 5828)
@@ -102,10 +102,7 @@
     }
 
-    psLogMsg ("psphot", 4, "load data: %f sec\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot"));
 
-    psphotSaveImage (NULL, weight, "weight.fits");
-    psphotSaveImage (NULL, mask, "mask.fits");
-
-    // save the image data & return it
+    // return image data
     eamReadout *imdata = eamReadoutAlloc(image, weight, mask, header);
     return (imdata);
Index: /trunk/psphot/src/psphotSourceFits.c
===================================================================
--- /trunk/psphot/src/psphotSourceFits.c	(revision 5828)
+++ /trunk/psphot/src/psphotSourceFits.c	(revision 5828)
@@ -0,0 +1,287 @@
+# include "psphot.h"
+# include "psEllipse.h"
+
+// given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
+
+bool psphotFitPSF (eamReadout *imdata, pmSource *source) { 
+
+    float x, y;
+
+    // save the PSF model from the Ensemble fit
+    pmModel *PSF = pmModelCopy (source->modelPSF);
+
+    // extend source radius as needed
+    psphotCheckRadiusPSF (imdata, source, PSF);
+
+    x = PSF->params->data.F32[2];
+    y = PSF->params->data.F32[3];
+
+    // fit PSF model (set/unset the pixel mask)
+    psImageKeepCircle (source->mask, x, y, PSF->radius, "OR", PSPHOT_MASK_MARKED);
+    pmSourceFitModel (source, PSF, true);
+    psImageKeepCircle (source->mask, x, y, PSF->radius, "AND", ~PSPHOT_MASK_MARKED);
+
+    // does the PSF model succeed?
+    if (!psphotEvalPSF (source, PSF)) {
+	psFree (PSF);
+	return false;
+    }
+
+    psTrace ("psphot.blend", 5, "fitted as PSF\n");
+    pmSourceSubModel (source->pixels, source->mask, PSF, false, false);
+
+    // free old model, save new model 
+    psFree (source->modelPSF);
+    source->modelPSF = PSF;
+
+    source->mode |=  PM_SOURCE_SUBTRACTED;
+    source->mode &= ~PM_SOURCE_TEMPSUB;
+    return true;
+}
+
+
+static float GAL_MIN_SN;
+static float GAL_MOMENTS_RAD;
+static pmModelType modelTypeFLT;
+
+bool psphotInitLimitsFLT (psMetadata *config, psStats *sky) {
+
+    bool status;
+
+    // galaxy model parameters
+    GAL_MIN_SN       = psMetadataLookupF32 (&status, config, "GAL_MIN_SN");
+    GAL_MOMENTS_RAD  = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
+
+    // 'galaxy' model descriptions
+    char *modelNameFLT = psMetadataLookupPtr (&status, config, "GAL_MODEL");
+    modelTypeFLT = pmModelSetType (modelNameFLT);
+    psphotInitRadiusFLT (config, sky, modelTypeFLT);
+
+    return true;
+}
+
+bool psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources) { 
+
+    bool okFLT, okDBL;
+    float chiFLT, chiDBL;
+
+    // skip the source if we don't think it is extended
+    if (source->type == PM_SOURCE_UNKNOWN) return false;
+    if (source->type == PM_SOURCE_DEFECT) return false;
+    if (source->type == PM_SOURCE_SATURATED) return false;
+    if (source->moments->SN < GAL_MIN_SN) return false;
+
+    // recalculate the source moments using the larger galaxy moments radius
+    if (!pmSourceMoments (source, GAL_MOMENTS_RAD)) return false;
+
+    psTrace ("psphot.blend", 5, "trying blob...\n");
+
+    pmSource *tmpSrc = pmSourceAlloc ();
+
+    pmModel *FLT = psphotFitFLT (imdata, source);
+    okFLT = psphotEvalFLT (tmpSrc, FLT);
+    chiFLT = FLT->chisq / FLT->nDOF;
+
+    psArray *DBL = psphotFitDBL (imdata, source);
+    okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
+    okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
+    pmModel *ONE = DBL->data[0];
+    chiDBL = ONE->chisq / ONE->nDOF;
+
+    if (okFLT && okDBL) {
+	psTrace ("psphot.blend", 5, "blob chisq: %f vs %f\n", chiFLT, chiDBL);
+	// XXX EAM : a bogus bias: need to examine this better
+	if (3*chiFLT > chiDBL) goto keepDBL;
+	goto keepFLT;
+    }
+
+    if (okFLT && !okDBL) goto keepFLT;
+    if (!okFLT && okDBL) goto keepDBL;
+    return false;
+
+keepFLT:
+    // sub FLT
+    pmSourceSubModel (source->pixels, source->mask, FLT, false, false);
+    psTrace ("psphot.blend", 5, "blob as FLT: %f %f\n", FLT->params->data.F32[2], FLT->params->data.F32[3]);
+
+    // save new model
+    source->modelFLT = FLT;
+    source->mode |=  PM_SOURCE_SUBTRACTED;
+    source->mode &= ~PM_SOURCE_TEMPSUB;
+    return true;
+
+keepDBL:
+    // sub DLB
+    pmSourceSubModel (source->pixels, source->mask, (pmModel *) DBL->data[0], false, false);
+    pmSourceSubModel (source->pixels, source->mask, (pmModel *) DBL->data[1], false, false);
+    psTrace ("psphot.blend", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[2], ONE->params->data.F32[3]);
+
+    source->modelPSF = (pmModel *) DBL->data[0];
+    source->mode |=  PM_SOURCE_SUBTRACTED;
+    source->mode &= ~PM_SOURCE_TEMPSUB;
+
+    // save new model
+    // tmpSrc->modelPSF = (pmModel *) modelSet->data[1];
+    // psArrayAdd (sources, 100, tmpSrc);
+
+    return true;
+}
+
+// fit a double PSF source to an extended blob
+psArray *psphotFitDBL (eamReadout *imdata, pmSource *source) { 
+
+    float x, y, dx, dy;
+    pmModel *DBL;
+    pmModel *PSF;
+    EllipseAxes axes;
+    EllipseMoments moments;
+    psArray *modelSet;
+    
+    // make a guess at the position of the two sources
+    moments.x2 = source->moments->Sx;
+    moments.y2 = source->moments->Sy;
+    moments.xy = source->moments->Sxy;
+    axes = EllipseMomentsToAxes (moments);
+
+    dx = 2 * cos (axes.theta);
+    dy = 2 * sin (axes.theta);
+
+    // save the PSF model from the Ensemble fit
+    PSF = source->modelPSF;
+    psphotCheckRadiusPSF (imdata, source, PSF);
+
+    modelSet = psArrayAlloc (2);
+
+    DBL = pmModelCopy (PSF);
+    DBL->params->data.F32[1] *= 0.5;
+    DBL->params->data.F32[2] += dx;
+    DBL->params->data.F32[3] += dy;
+    modelSet->data[0] = DBL;
+
+    DBL = pmModelCopy (PSF);
+    DBL->params->data.F32[1] *= 0.5;
+    DBL->params->data.F32[2] -= dx;
+    DBL->params->data.F32[3] -= dy;
+    modelSet->data[1] = DBL;
+
+    x = PSF->params->data.F32[2];
+    y = PSF->params->data.F32[3];
+
+    // fit FLT (not PSF) model (set/unset the pixel mask)
+    psImageKeepCircle (source->mask, x, y, PSF->radius, "OR", PSPHOT_MASK_MARKED);
+    pmSourceFitSet (source, modelSet, true);
+    psImageKeepCircle (source->mask, x, y, PSF->radius, "AND", ~PSPHOT_MASK_MARKED);
+
+    return (modelSet);
+}
+
+pmModel *psphotFitFLT (eamReadout *imdata, pmSource *source) { 
+
+    float x, y;
+
+    // use the source moments, etc to guess basic model parameters
+    pmModel *FLT = pmSourceModelGuess (source, modelTypeFLT); 
+
+    psphotCheckRadiusFLT (imdata, source, FLT);
+
+    x = FLT->params->data.F32[2];
+    y = FLT->params->data.F32[3];
+
+    // fit FLT (not PSF) model (set/unset the pixel mask)
+    psImageKeepCircle (source->mask, x, y, FLT->radius, "OR", PSPHOT_MASK_MARKED);
+    pmSourceFitModel (source, FLT, false);
+    psImageKeepCircle (source->mask, x, y, FLT->radius, "AND", ~PSPHOT_MASK_MARKED);
+
+    return (FLT);
+}
+
+bool psphotFitBlend (eamReadout *imdata, pmSource *source) { 
+
+    float x, y, dR;
+
+    // if this source is not a possible blend, just fit as PSF
+    if ((source->blends == NULL) || (source->mode & PM_SOURCE_SATSTAR)) {
+	bool status = psphotFitPSF (imdata, source);
+	return status;
+    }
+
+    psTrace ("psphot.blend", 5, "trying blend...\n");
+
+    // save the PSF model from the Ensemble fit
+    pmModel *PSF = pmModelCopy (source->modelPSF);
+
+    // extend source radius as needed
+    psphotCheckRadiusPSF (imdata, source, PSF);
+
+    x = PSF->params->data.F32[2];
+    y = PSF->params->data.F32[3];
+
+    psArray *modelSet = psArrayAlloc (source->blends->n + 1);
+    modelSet->n = 0;
+    psArrayAdd (modelSet, 16, PSF);
+    
+    psArray *sourceSet = psArrayAlloc (source->blends->n + 1);
+    sourceSet->n = 0;
+    psArrayAdd (sourceSet, 16, source);
+
+    for (int i = 0; i < source->blends->n; i++) {
+	pmSource *blend = source->blends->data[i];
+
+	// is this object close enough to include in the fit
+	// XXX this test is bogus: need to think about this
+	dR = hypot (blend->peak->x - x, blend->peak->y - y);
+	if (dR > PSF->radius - 3) continue;
+
+	// create the model and guess parameters for this blend
+	blend->modelPSF = pmModelAlloc (PSF->type);
+	pmModel *model = blend->modelPSF;
+
+	for (psS32 i = 0; i < model->params->n; i++) {
+	    model->params->data.F32[i] = PSF->params->data.F32[i];
+	    model->dparams->data.F32[i] = PSF->dparams->data.F32[i];
+	}
+
+	model->params->data.F32[1] = blend->peak->counts;
+	model->params->data.F32[2] = blend->peak->x;
+	model->params->data.F32[3] = blend->peak->y;
+
+	// add this blend to the list
+	psArrayAdd (modelSet, 16, model);
+	psArrayAdd (sourceSet, 16, blend);
+    }
+
+    // fit PSF model (set/unset the pixel mask)
+    psImageKeepCircle (source->mask, x, y, PSF->radius, "OR", PSPHOT_MASK_MARKED);
+    pmSourceFitSet (source, modelSet, true);
+    psImageKeepCircle (source->mask, x, y, PSF->radius, "AND", ~PSPHOT_MASK_MARKED);
+
+    // evaluate the blend objects, subtract if good, free otherwise
+    for (int i = 1; i < modelSet->n; i++) {
+	pmSource *src = sourceSet->data[i];
+	pmModel *model = modelSet->data[i];
+
+	if (!psphotEvalPSF (src, model)) {
+	    psFree (model);
+	    src->modelPSF = NULL;
+	    continue;
+	}
+
+	psTrace ("psphot.blend", 5, "fitted blend as PSF\n");
+	pmSourceSubModel (source->pixels, source->mask, model, false, false);
+	src->mode |=  PM_SOURCE_SUBTRACTED;
+	src->mode &= ~PM_SOURCE_TEMPSUB;
+    }
+
+    // evaluate the primary object
+    if (!psphotEvalPSF (source, PSF)) {
+	psFree (PSF);
+	return false;
+    }
+
+    psTrace ("psphot.blend", 5, "fitted primary as PSF\n");
+    pmSourceSubModel (source->pixels, source->mask, PSF, false, false);
+    source->modelPSF = PSF;
+    source->mode |=  PM_SOURCE_SUBTRACTED;
+    source->mode &= ~PM_SOURCE_TEMPSUB;
+    return true;
+}