Index: /trunk/psphot/src/psphotSourceFits.c
===================================================================
--- /trunk/psphot/src/psphotSourceFits.c	(revision 7436)
+++ /trunk/psphot/src/psphotSourceFits.c	(revision 7437)
@@ -2,4 +2,110 @@
 
 // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
+
+bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { 
+
+    float x, y, dR;
+
+    // if this source is not a possible blend, just fit as PSF
+    if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
+	bool status = psphotFitPSF (readout, source, psf);
+	return status;
+    }
+    psTrace ("psphot.blend", 5, "trying blend...\n");
+
+    // save the PSF model from the Ensemble fit
+    pmModel *PSF = pmModelCopy (source->modelPSF);
+    if (isnan(PSF->params->data.F32[1])) psAbort ("psphot", "nan in blend fit primary");
+
+    x = PSF->params->data.F32[2];
+    y = PSF->params->data.F32[3];
+
+    psArray *modelSet = psArrayAlloc (source->blends->n + 1);
+    psArrayAdd (modelSet, 16, PSF);
+    
+    psArray *sourceSet = psArrayAlloc (source->blends->n + 1);
+    psArrayAdd (sourceSet, 16, source);
+
+    // we need to include all blends in the fit (unless primary is saturated?)
+    dR = 0;
+    for (int i = 0; i < source->blends->n; i++) {
+	pmSource *blend = source->blends->data[i];
+
+	// find the blend which is furthest from source
+	dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y));
+
+	// create the model and guess parameters for this blend
+	pmModel *model = pmModelAlloc (PSF->type);
+	for (int j = 0; j < model->params->n; j++) {
+	    model->params->data.F32[j] = PSF->params->data.F32[j];
+	    model->dparams->data.F32[j] = PSF->dparams->data.F32[j];
+	}
+
+	// XXX assume local sky is 0.0?
+	model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky;
+	if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit");
+	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);
+
+	// free to avoid double counting model 
+	psFree (model);
+    }
+
+    // extend source radius as needed
+    psphotCheckRadiusPSFBlend (readout, source, PSF, dR);
+
+    // fit PSF model (set/unset the pixel mask)
+    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK);
+    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
+    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
+
+    // correct model chisq for flux trend
+    double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
+    PSF->chisqNorm = PSF->chisq / chiTrend;
+
+    // evaluate the blend objects, subtract if good, free otherwise
+    for (int i = 1; i < modelSet->n; i++) {
+	pmSource *blend = sourceSet->data[i];
+	pmModel *model  = modelSet->data[i];
+
+	// correct model chisq for flux trend
+	chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]);
+	model->chisqNorm = model->chisq / chiTrend;
+
+	// if this one failed, skip it
+	if (!psphotEvalPSF (blend, model)) {
+	    psTrace ("psphot.blend", 5, "failed on blend %d of %d\n", i, modelSet->n);
+	    continue;
+	}
+
+	// otherwise, supply the resulting model to the corresponding blend
+	blend->modelPSF = psMemIncrRefCounter (model);
+	psTrace ("psphot.blend", 5, "fitted blend as PSF\n");
+	pmModelSub (source->pixels, source->mask, model, false, false);
+	blend->mode |=  PM_SOURCE_MODE_SUBTRACTED;
+	blend->mode &= ~PM_SOURCE_MODE_TEMPSUB;
+    }
+    psFree (modelSet);
+    psFree (sourceSet);
+
+    // evaluate the primary object
+    if (!psphotEvalPSF (source, PSF)) {
+	psTrace ("psphot.blend", 5, "failed on blend 0 of %d\n", modelSet->n);
+	psFree (PSF);
+	return false;
+    }
+
+    psTrace ("psphot.blend", 5, "fitted primary as PSF\n");
+    pmModelSub (source->pixels, source->mask, PSF, false, false);
+    psFree (source->modelPSF);
+    source->modelPSF = PSF;
+    source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
+    source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
+    return true;
+}
 
 bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { 
@@ -26,5 +132,4 @@
     chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
     PSF->chisqNorm = PSF->chisq / chiTrend;
-    // PSF->chisqNorm = PSF->chisq;
 
     // does the PSF model succeed?
@@ -43,7 +148,7 @@
     source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
     source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
-    return true;
-}
-
+
+    return true;
+}
 
 static float EXT_MIN_SN;
@@ -244,107 +349,2 @@
 }
 
-bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { 
-
-    float x, y, dR;
-
-    // if this source is not a possible blend, just fit as PSF
-    if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
-	bool status = psphotFitPSF (readout, source, psf);
-	return status;
-    }
-    psTrace ("psphot.blend", 5, "trying blend...\n");
-
-    // save the PSF model from the Ensemble fit
-    pmModel *PSF = pmModelCopy (source->modelPSF);
-    if (isnan(PSF->params->data.F32[1])) psAbort ("psphot", "nan in blend fit primary");
-
-    x = PSF->params->data.F32[2];
-    y = PSF->params->data.F32[3];
-
-    psArray *modelSet = psArrayAlloc (source->blends->n + 1);
-    psArrayAdd (modelSet, 16, PSF);
-    
-    psArray *sourceSet = psArrayAlloc (source->blends->n + 1);
-    psArrayAdd (sourceSet, 16, source);
-
-    // we need to include all blends in the fit (unless primary is saturated?)
-    dR = 0;
-    for (int i = 0; i < source->blends->n; i++) {
-	pmSource *blend = source->blends->data[i];
-
-	// find the blend which is furthest from source
-	dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y));
-
-	// create the model and guess parameters for this blend
-	pmModel *model = pmModelAlloc (PSF->type);
-	for (int j = 0; j < model->params->n; j++) {
-	    model->params->data.F32[j] = PSF->params->data.F32[j];
-	    model->dparams->data.F32[j] = PSF->dparams->data.F32[j];
-	}
-
-	// XXX assume local sky is 0.0?
-	model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky;
-	if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit");
-	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);
-
-	// free to avoid double counting model 
-	psFree (model);
-    }
-
-    // extend source radius as needed
-    psphotCheckRadiusPSFBlend (readout, source, PSF, dR);
-
-    // fit PSF model (set/unset the pixel mask)
-    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK);
-    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
-    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
-
-    // correct model chisq for flux trend
-    double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
-    PSF->chisqNorm = PSF->chisq / chiTrend;
-
-    // evaluate the blend objects, subtract if good, free otherwise
-    for (int i = 1; i < modelSet->n; i++) {
-	pmSource *blend = sourceSet->data[i];
-	pmModel *model  = modelSet->data[i];
-
-	// correct model chisq for flux trend
-	chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]);
-	model->chisqNorm = model->chisq / chiTrend;
-
-	// if this one failed, skip it
-	if (!psphotEvalPSF (blend, model)) {
-	    psTrace ("psphot.blend", 5, "failed on blend %d of %d\n", i, modelSet->n);
-	    continue;
-	}
-
-	// otherwise, supply the resulting model to the corresponding blend
-	blend->modelPSF = psMemIncrRefCounter (model);
-	psTrace ("psphot.blend", 5, "fitted blend as PSF\n");
-	pmModelSub (source->pixels, source->mask, model, false, false);
-	blend->mode |=  PM_SOURCE_MODE_SUBTRACTED;
-	blend->mode &= ~PM_SOURCE_MODE_TEMPSUB;
-    }
-    psFree (modelSet);
-    psFree (sourceSet);
-
-    // evaluate the primary object
-    if (!psphotEvalPSF (source, PSF)) {
-	psTrace ("psphot.blend", 5, "failed on blend 0 of %d\n", modelSet->n);
-	psFree (PSF);
-	return false;
-    }
-
-    psTrace ("psphot.blend", 5, "fitted primary as PSF\n");
-    pmModelSub (source->pixels, source->mask, PSF, false, false);
-    psFree (source->modelPSF);
-    source->modelPSF = PSF;
-    source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
-    source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
-    return true;
-}