Index: trunk/psphot/src/psphotSparseMatrix.c
===================================================================
--- trunk/psphot/src/psphotSparseMatrix.c	(revision 5653)
+++ trunk/psphot/src/psphotSparseMatrix.c	(revision 5654)
@@ -1,23 +1,70 @@
 # include "psphot.h"
 
-psSparse *psphotStarOverlaps (psArray *sources) {
+bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
 
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
-    psVector *params = model->params;
+    bool  status;
+    float x;
+    float y;
+    float Sky;
+    int   Nfit = 0;
 
-    // source analysis is done in S/N order (brightest first)
+    psTimerStart ("psphot");
+
+    // source analysis is done in spatial order
     sources = psArraySort (sources, psphotSortByY);
 
-    // fill out the models for each object
-    // XXX need to assign the model parameters using the PSF function, with unit peak
+    float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
+    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
+    float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
+
+    // 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);
+
+    // vector for model input coordinates
+    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
+
+    // this function specifies the radius at this the model hits the given flux
+    pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type);
+
+    // this is the function to get the model function
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+
+    // pre-calculate all model pixels
     psArray  *pixlist = psArrayAlloc (sources->n);
     for (int k = 0; k < sources->n; k++) {
-	src = sources->data[k];
+	source = sources->data[k];
+
+	// 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);
+	// model->params->data.F32[2] = source->peak->x;
+	// model->params->data.F32[3] = source->peak->y;
+	model->params->data.F32[1] = 1.0;
+	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");
+	}
+
+	// build the model image 
 	nx = src->pixels->numCols;
 	ny = src->pixels->numRows;
 	im = psImageAlloc (nx, ny, PS_DATA_F32);
+	mk = psImageAlloc (nx, ny, PS_DATA_U8);
+
+	// fill in the model pixel values
+	psImageKeepCircle (mk, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
+	psVector *params = model->params;
 	for (int i = 0; i < nx; i++) {
 	    for (int j = 0; j < ny; j++) {
+		if (mk->data.U8[j][i]) continue;
 		x->data.F32[0] = (float) i + pix->col0;
 		x->data.F32[1] = (float) j + pix->row0;
@@ -28,8 +75,10 @@
     }
     
-    psVector *Aij = psVectorAlloc (100, PS_DATA_F32);
-    psVector *Bfj = psVectorAlloc (100, PS_DATA_F32);
-    psVector *Si  = psVectorAlloc (100, PS_DATA_S32);
-    psVector *Sj  = psVectorAlloc (100, PS_DATA_S32);
+    // fill out the sparse matrix
+    psSparse *sparse = psSparseAlloc (sources->n, 100);
+    psVector *Aij = sparse->Aij;
+    psVector *Bfj = sparse->Bfj;
+    psVector *Si  = sparse->Si;
+    psVector *Sj  = sparse->Sj;
 
     int Ri, Rj;
@@ -42,14 +91,9 @@
 	// diagonal element (auto-cross-product)
 	f = psphotCrossProduct (pi, pi);
-	Aij->data.F32[Nelem] = f;
-	Si->data.S32[Nelem] = i;
-	Sj->data.S32[Nelem] = j;
-	Nelem ++;
-
-	Qii->data.F32[i] = f;
+	psSparseMatrixElement (sparse, i, i, f);
 
 	// find the image x model value
 	f = psphotCrossProduct (pi, pf);
-	Bfj->data.F32[i] = f;
+	psSparseVectorElement (sparse, i, f);
 
 	// loop over all other stars following this one
@@ -64,40 +108,9 @@
 	    // got an overlap; calculate cross-product and add to output array
 	    f = psphotCrossProduct (pi, pj);
-	    
-	    Aij->data.F32[Nelem] = f;
-	    Si->data.S32[Nelem] = i;
-	    Sj->data.S32[Nelem] = j;
-	    Nelem ++;
-
-	    Aij->data.F32[Nelem] = f;
-	    Si->data.S32[Nelem] = j;
-	    Sj->data.S32[Nelem] = i;
-	    Nelem ++;
+	    psSparseMatrixElement (sparse, j, i, f); 
 	}
     }
-
-    index = psVectorSortIndex (NULL, Sj);
-
-    psVector *tAij = psVectorAlloc (Nelem, PS_DATA_F32);
-    psVector *tSi  = psVectorAlloc (Nelem, PS_DATA_S32);
-    psVector *tSj  = psVectorAlloc (Nelem, PS_DATA_S32);
-    for (i = 0; i < Nelem; i++) {
-	j = index->data.U32[i];
-	tAij->data.F32[i] = Aij->data.F32[j];
-	tSi->data.S32[i]  = Si->data.S32[j];
-	tSj->data.S32[i]  = Sj->data.S32[j];
-    }
-    psFree (Aij);
-    psFree (Si);
-    psFree (Sj);
-
-    sparse->Aij = Aij;
-    sparse->Bfj = Bfj;
-    sparse->Qii = Qii;
-    sparse->Si = Si;
-    sparse->Sj = Sj;
-    sparse->Nelem = Nelem;
-    sparse->Nrows = sources->n;
-
+    
+    psSparseResort (sparse);
     return (match);
 }