Index: /trunk/psphot/src/pmSourceFitFixed.c
===================================================================
--- /trunk/psphot/src/pmSourceFitFixed.c	(revision 5629)
+++ /trunk/psphot/src/pmSourceFitFixed.c	(revision 5629)
@@ -0,0 +1,88 @@
+# include "psphot.h"
+
+bool pmSourceFitFixed (pmSource *source,
+		       pmModel *model)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->moments, 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);
+
+    // XXX EAM : is it necessary for the mask & weight to exist?  the
+    //           tests below could be conditions (!NULL)
+
+    psBool fitStatus = true;
+    psBool onPic     = true;
+    psBool rc        = true;
+    psF32  Ro, ymodel;
+
+    psVector *params = model->params;
+    psVector *dparams = model->dparams;
+    psVector *paramMask = NULL;
+
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+
+    psF32 Sky = params->data.F32[0];
+
+    // construct the coordinate and value entries
+    psVector *x     = psVectorAlloc(100, PS_TYPE_F32);
+    psVector *y     = psVectorAlloc(100, PS_TYPE_F32);
+    psVector *yErr  = psVectorAlloc(100, PS_TYPE_F32);
+    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+    x->n = y->n = yErr->n = 0;
+    psS32 count = 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;
+	    }
+
+	    coord->data.F32[0] = (psF32) (j + source->pixels->col0);
+	    coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+	    ymodel = modelFunc (NULL, params, coord);
+
+	    x->data.F32[count] = source->pixels->data.F32[i][j];
+	    y->data.F32[count] = ymodel;
+
+	    // this test modifies the weight based on deviation from the model flux
+	    Ro = 1.0 + fabs (x->data.F32[count] - ymodel) / sqrt(PS_SQR(ymodel - Sky) + PS_SQR(Sky));
+	    yErr->data.F32[count] = sqrt (source->weight->data.F32[i][j] * Ro);
+
+	    psVectorExtend (x, 1, 100);
+	    psVectorExtend (y, 1, 100);
+	    psVectorExtend (yErr, 1, 100);
+	    count++;
+        }
+    }
+
+    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+
+    // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    // stats->clipIter = 3;
+    // poly = psVectorClipFitPolynomial1D (poly, stats, NULL, 0, y, yErr, x);
+
+    poly = psVectorFitPolynomial1D (poly, stats, NULL, 0, y, yErr, x);
+
+    params->data.F32[0] = poly->coeff[0];
+    params->data.F32[1] = poly->coeff[1];
+
+    // XXX need to get errors back from psVectorClipFit...
+    dparams->data.F32[0] = 0;
+    dparams->data.F32[1] = 0;
+
+    // XXX need to get chisq back from psVectorClipFit...
+    model->chisq = 0;
+    model->nIter = 0;
+    model->nDOF  = 0;
+
+    psFree(x);
+    psFree(y);
+    psFree(yErr);
+    psFree(stats);
+
+    return true;
+}
+
Index: /trunk/psphot/src/psphotFixedPSF.c
===================================================================
--- /trunk/psphot/src/psphotFixedPSF.c	(revision 5629)
+++ /trunk/psphot/src/psphotFixedPSF.c	(revision 5629)
@@ -0,0 +1,91 @@
+# include "psphot.h"
+
+// fit psf model to all objects with centroid fixed 
+
+bool psphotFitFixedPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
+{ 
+    bool  status;
+    float x;
+    float y;
+    int   Nfit = 0;
+    int   Nsub = 0;
+    int   Niter = 0;
+
+    psTimerStart ("psphot");
+
+    // we may set this differently here from the value used to mark likely saturated stars
+    float SATURATION   = psMetadataLookupF32 (&status, config, "SATURATION");
+    float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
+
+    float PSF_MIN_SN   	   = psMetadataLookupF32 (&status, config, "PSF_MIN_SN");
+    float PSF_MAX_CHI  	   = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI");
+    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
+    float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
+    float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
+
+    // set the object surface-brightness limit for fitted pixels
+    float FLUX_LIMIT  = PSF_FIT_NSIGMA * sky->sampleStdev;
+    psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
+
+    // this function specifies the radius at this the model hits the given flux
+    pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type);
+
+    for (int i = 0; i < sources->n; i++) {
+
+	pmSource *source = sources->data[i];
+
+	// skip non-astronomical objects (very likely defects)
+	// XXX EAM : should we try these anyway?
+	if (source->type == PM_SOURCE_DEFECT) continue; 
+	if (source->type == PM_SOURCE_SATURATED) continue;
+
+	// use the source moments, etc to guess basic model parameters
+	pmModel *modelFLT = pmSourceModelGuess (source, psf->type); 
+
+	// set PSF parameters for this model
+	pmModel *model = pmModelFromPSF (modelFLT, psf);
+	x = model->params->data.F32[2];
+	y = model->params->data.F32[3];
+	psFree (modelFLT);
+
+	// set the fit radius based on the object flux limit and the model
+	// XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median)
+	model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING;
+	if (isnan(model->radius)) {
+	  psAbort ("apply_psf_model", "error in radius");
+	}
+	
+	// check if we need to redefine the pixels
+	// XXX EAM : a better test would examine the source pixels
+	if (model->radius > OUTER_RADIUS) {
+	  // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER)
+	    psphotDefinePixels (source, imdata, x, y, model->radius);
+	}
+
+	// fit PSF model (set/unset the pixel mask)
+	psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
+	status = pmSourceFitFixed (source, model);
+	psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED);
+
+	if (!status || (model->params->data.F32[1] < 0)) {
+	  psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
+	  source->type = PM_SOURCE_FAIL_FIT_PSF;  // better choice?
+	  psFree (model);
+	  continue;
+	}
+
+	// model succeeded : tentatively keep it
+	source->modelPSF = model;
+	Niter += model[0].nIter;
+	Nfit ++;
+
+	// subtract object, leave local sky
+	sky = model->params->data.F32[0];
+	model->params->data.F32[0] = 0;
+	pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false);
+	model->params->data.F32[0] = sky;
+    }
+
+    psLogMsg ("psphot", 3, "fit fixed PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
+    return (true);
+}