Index: trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- trunk/psphot/src/psphotEnsemblePSF.c	(revision 9576)
+++ trunk/psphot/src/psphotEnsemblePSF.c	(revision 9595)
@@ -4,5 +4,5 @@
 // 2006.02.07 : no leaks!
 // fit all reasonable sources with the linear PSF model
-bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { 
+bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
 
     bool  status;
@@ -39,103 +39,103 @@
 
     const bool CONSTANT_PHOTOMETRIC_WEIGHTS =
-	psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
+        psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
     if (!status) {
-	psAbort(PS_FILE_LINE, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
+        psAbort(PS_FILE_LINE, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
     }
 
     for (int i = 0; i < sources->n; i++) {
-	pmSource *inSource = sources->data[i];
-
-	// skip non-astronomical objects (very likely defects)
-	// XXX EAM : should we try these anyway?
-	if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue; 
-	if (inSource->type == PM_SOURCE_TYPE_SATURATED) continue;
-	if (final) {
-	    if (inSource->mode &  PM_SOURCE_MODE_SUBTRACTED) continue;
-	} else {
-	    if (inSource->mode &  PM_SOURCE_MODE_BLEND) continue;
-	}
-
-	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 ();
-
-	// really saturated stars should be re-measured for a better centroid
-	// XXX EAM : move this to a 'clear satstar function'
-	if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
-	    status = pmSourceMoments (inSource, INNER_RADIUS);
-	}
-
-	// XXX EAM : add option to use EXT or PSF form
-	// use the source moments, etc to guess basic model parameters
-	pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); 
-	if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
-	    modelEXT->params->data.F32[PM_PAR_XPOS] = inSource->moments->x;
-	    modelEXT->params->data.F32[PM_PAR_YPOS] = inSource->moments->y;
-	} else {
-	    // peak-up on peak (for non-sat objects)
-
-	    // ix,iy must land on inSource->pixels 
-	    int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x));
-	    int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y));
-
-	    psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
-	    psPlane min = psImageBicubeMin (bicube);
-
-	    psTrace ("psphot", 5, "peak coord: %f %f -> %f %f\n", 
-		     modelEXT->params->data.F32[PM_PAR_XPOS], modelEXT->params->data.F32[PM_PAR_YPOS], 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)) {
-		modelEXT->params->data.F32[PM_PAR_XPOS] = min.x + ix;
-		modelEXT->params->data.F32[PM_PAR_YPOS] = min.y + iy;
-	    }
-	    psFree (bicube);
-	}
-
-	// set PSF parameters for this model
-	pmModel *model = pmModelFromPSF (modelEXT, psf);
-	psFree (modelEXT);
-
-	// save the original coords
-	x = model->params->data.F32[PM_PAR_XPOS];
-	y = model->params->data.F32[PM_PAR_YPOS];
-
-	// set the fit radius based on the object flux limit and the model
-	psphotCheckRadiusPSF (readout, inSource, model);
-
-	// make temporary copies of the image pixels and mask
-	otSource->mask   = psImageCopy (NULL, inSource->mask,   PS_TYPE_U8);
-	otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32);
-	otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32);
-
-	// build the model image 
-	psImage *flux = otSource->pixels;
-	psImage *mask = otSource->mask;
-
-	// set model to unit peak, zero sky (we assume sky is subtracted)
-	model->params->data.F32[PM_PAR_SKY] = 0.0;
-	model->params->data.F32[PM_PAR_I0] = 1.0;
-
-	// fill in the model pixel values
-	psImageInit (flux, 0.0);
-	psImageKeepCircle (mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK);
-	pmModelAdd (flux, mask, model, false, false);
-
-	// calculate nDOF (nPix - 1)
-	// int Nmaskpix = psImageCountPixelMask (mask, allArray, PM_MASK_SAT);
-	// model->nDOF  = mask->numCols*mask->numRows - Nmaskpix - 1;
-
-	// save source in list
-	otSource->modelPSF = model;
-	index->data.U32[models->n] = i;
-	psArrayAdd (models, 100, otSource);
-	psFree (otSource);
+        pmSource *inSource = sources->data[i];
+
+        // skip non-astronomical objects (very likely defects)
+        // XXX EAM : should we try these anyway?
+        if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue;
+        if (inSource->type == PM_SOURCE_TYPE_SATURATED) continue;
+        if (final) {
+            if (inSource->mode &  PM_SOURCE_MODE_SUBTRACTED) continue;
+        } else {
+            if (inSource->mode &  PM_SOURCE_MODE_BLEND) continue;
+        }
+
+        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 ();
+
+        // really saturated stars should be re-measured for a better centroid
+        // XXX EAM : move this to a 'clear satstar function'
+        if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
+            status = pmSourceMoments (inSource, INNER_RADIUS);
+        }
+
+        // XXX EAM : add option to use EXT or PSF form
+        // use the source moments, etc to guess basic model parameters
+        pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type);
+        if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
+            modelEXT->params->data.F32[PM_PAR_XPOS] = inSource->moments->x;
+            modelEXT->params->data.F32[PM_PAR_YPOS] = inSource->moments->y;
+        } else {
+            // peak-up on peak (for non-sat objects)
+
+            // ix,iy must land on inSource->pixels
+            int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x));
+            int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y));
+
+            psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy);
+            psPlane min = psImageBicubeMin (bicube);
+
+            psTrace ("psphot", 5, "peak coord: %f %f -> %f %f\n",
+                     modelEXT->params->data.F32[PM_PAR_XPOS], modelEXT->params->data.F32[PM_PAR_YPOS], 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)) {
+                modelEXT->params->data.F32[PM_PAR_XPOS] = min.x + ix;
+                modelEXT->params->data.F32[PM_PAR_YPOS] = min.y + iy;
+            }
+            psFree (bicube);
+        }
+
+        // set PSF parameters for this model
+        pmModel *model = pmModelFromPSF (modelEXT, psf);
+        psFree (modelEXT);
+
+        // save the original coords
+        x = model->params->data.F32[PM_PAR_XPOS];
+        y = model->params->data.F32[PM_PAR_YPOS];
+
+        // set the fit radius based on the object flux limit and the model
+        psphotCheckRadiusPSF (readout, inSource, model);
+
+        // make temporary copies of the image pixels and mask
+        otSource->mask   = psImageCopy (NULL, inSource->mask,   PS_TYPE_U8);
+        otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32);
+        otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32);
+
+        // build the model image
+        psImage *flux = otSource->pixels;
+        psImage *mask = otSource->mask;
+
+        // set model to unit peak, zero sky (we assume sky is subtracted)
+        model->params->data.F32[PM_PAR_SKY] = 0.0;
+        model->params->data.F32[PM_PAR_I0] = 1.0;
+
+        // fill in the model pixel values
+        psImageInit (flux, 0.0);
+        psImageKeepCircle (mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK);
+        pmModelAdd (flux, mask, model, false, false);
+
+        // calculate nDOF (nPix - 1)
+        // int Nmaskpix = psImageCountPixelMask (mask, allArray, PM_MASK_SAT);
+        // model->nDOF  = mask->numCols*mask->numRows - Nmaskpix - 1;
+
+        // save source in list
+        otSource->modelPSF = model;
+        index->data.U32[models->n] = i;
+        psArrayAdd (models, 100, otSource);
+        psFree (otSource);
     }
     psLogMsg ("psphot.emsemble", 4, "built models: %f (%ld objects)\n", psTimerMark ("psphot"), sources->n);
-    
+
     // fill out the sparse matrix
     psSparse *sparse = psSparseAlloc (models->n, 100);
@@ -144,32 +144,32 @@
 
     for (int i = 0; i < models->n; i++) {
-	int N = index->data.U32[i];
-	pmSource *Fi = sources->data[N];
-	pmSource *Mi = models->data[i];
-
-	// scale by diagonal element (auto-cross-product)
-	r = pmSourceCrossProduct (Mi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS);
-	weight->data.F32[i] = r;
-
-	psSparseMatrixElement (sparse, i, i, 1.0);
-
-	// find the image x model value
-	f = pmSourceCrossProduct (Fi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS);
-	psSparseVectorElement (sparse, i, f / r);
-
-	// loop over all other stars following this one
-	for (int j = i + 1; j < models->n; j++) {
-	    pmSource *Mj = models->data[j];
-
-	    // skip over disjoint source images, break after last possible overlap
-	    if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break;
-	    if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue;
-	    if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue;
-	    if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue;
-	    
-	    // got an overlap; calculate cross-product and add to output array
-	    f = pmSourceCrossProduct (Mi, Mj, CONSTANT_PHOTOMETRIC_WEIGHTS);
-	    psSparseMatrixElement (sparse, j, i, f / r); 
-	}
+        int N = index->data.U32[i];
+        pmSource *Fi = sources->data[N];
+        pmSource *Mi = models->data[i];
+
+        // scale by diagonal element (auto-cross-product)
+        r = pmSourceCrossProduct (Mi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        weight->data.F32[i] = r;
+
+        psSparseMatrixElement (sparse, i, i, 1.0);
+
+        // find the image x model value
+        f = pmSourceCrossProduct (Fi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        psSparseVectorElement (sparse, i, f / r);
+
+        // loop over all other stars following this one
+        for (int j = i + 1; j < models->n; j++) {
+            pmSource *Mj = models->data[j];
+
+            // skip over disjoint source images, break after last possible overlap
+            if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break;
+            if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue;
+            if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue;
+            if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue;
+
+            // got an overlap; calculate cross-product and add to output array
+            f = pmSourceCrossProduct (Mi, Mj, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            psSparseMatrixElement (sparse, j, i, f / r);
+        }
     }
     psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
@@ -186,40 +186,40 @@
     // adjust models, set sources and subtract
     for (int i = 0; i < models->n; i++) {
-	int N = index->data.U32[i];
-	pmSource *Fi = sources->data[N];
-	pmSource *Mi = models->data[i];
-
-	// if we already have a PSF model, free it.
-	psFree (Fi->modelPSF);
-
-	// need to increment counter so we can free models here and sources above
-	Fi->modelPSF = psMemIncrRefCounter (Mi->modelPSF);
-
-	// assign linearly-fitted normalization
-	if (isnan(norm->data.F32[i])) {
-	    psAbort ("psphot", "ensemble source is nan");
-	}
-	Fi->modelPSF->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
-	Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
-	// XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it
-
-	// subtract object
-	pmModelSub (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
-	Fi->mode |= PM_SOURCE_MODE_SUBTRACTED;
-	if (!final) Fi->mode |= PM_SOURCE_MODE_TEMPSUB;
+        int N = index->data.U32[i];
+        pmSource *Fi = sources->data[N];
+        pmSource *Mi = models->data[i];
+
+        // if we already have a PSF model, free it.
+        psFree (Fi->modelPSF);
+
+        // need to increment counter so we can free models here and sources above
+        Fi->modelPSF = psMemIncrRefCounter (Mi->modelPSF);
+
+        // assign linearly-fitted normalization
+        if (isnan(norm->data.F32[i])) {
+            psAbort ("psphot", "ensemble source is nan");
+        }
+        Fi->modelPSF->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
+        Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i]));
+        // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it
+
+        // subtract object
+        pmModelSub (Fi->pixels, Fi->mask, Fi->modelPSF, false, false);
+        Fi->mode |= PM_SOURCE_MODE_SUBTRACTED;
+        if (!final) Fi->mode |= PM_SOURCE_MODE_TEMPSUB;
     }
 
     // measure chisq for each source
     for (int i = 0; final && (i < models->n); i++) {
-	int N = index->data.U32[i];
-	pmSource *Fi = sources->data[N];
-	pmModel *model = Fi->modelPSF;
-
-	x = model->params->data.F32[PM_PAR_XPOS];
-	y = model->params->data.F32[PM_PAR_YPOS];
-
-	psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK);
-	pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight);
-	psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
+        int N = index->data.U32[i];
+        pmSource *Fi = sources->data[N];
+        pmModel *model = Fi->modelPSF;
+
+        x = model->params->data.F32[PM_PAR_XPOS];
+        y = model->params->data.F32[PM_PAR_YPOS];
+
+        psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "OR", PM_MASK_MARK);
+        pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight);
+        psImageKeepCircle (Fi->mask, x, y, model->radiusTMP, "AND", PS_NOT_U8(PM_MASK_MARK));
     }
 
@@ -239,10 +239,10 @@
     int Npix = 0;
     for (int j = 0; j < image->numRows; j++) {
-	for (int i = 0; i < image->numCols; i++) {
-	    if (mask->data.U8[j][i]) continue;
-	    if (weight->data.F32[j][i] <= 0) continue;
-	    dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i];
-	    Npix ++;
-	}
+        for (int i = 0; i < image->numCols; i++) {
+            if (mask->data.U8[j][i]) continue;
+            if (weight->data.F32[j][i] <= 0) continue;
+            dC += PS_SQR (image->data.F32[j][i]) / weight->data.F32[j][i];
+            Npix ++;
+        }
     }
     model->nDOF = Npix - 1;