Index: /trunk/psphot/src/psSparse.c
===================================================================
--- /trunk/psphot/src/psSparse.c	(revision 5653)
+++ /trunk/psphot/src/psSparse.c	(revision 5654)
@@ -1,13 +1,3 @@
 # include "psphot.h"
-
-typedef struct {
-    psVector *Aij;
-    psVector *Bfj;
-    psVector *Qii;
-    psVector *Si;
-    psVector *Sj;
-    int Nelem;
-    int Nrows;
-} psSparse;
 
 void psSparseMatrixTest () {
@@ -16,12 +6,12 @@
     psSparse *sparse = psSparseAlloc (3, 9);
 
-    psSparseMatrixElement (sparse, 0, 0, 1.0);
-    psSparseMatrixElement (sparse, 1, 1, 1.0);
+    psSparseMatrixElement (sparse, 0, 0, 3.0);
+    psSparseMatrixElement (sparse, 1, 1, 2.0);
     psSparseMatrixElement (sparse, 2, 2, 1.0);
 
-    psSparseMatrixElement (sparse, 0, 1, 0.1);
-    psSparseMatrixElement (sparse, 0, 2, 0.1);
-
-    psSparseMatrixResort (sparse);
+    psSparseMatrixElement (sparse, 1, 0, 0.1);
+    psSparseMatrixElement (sparse, 2, 0, -0.1);
+
+    psSparseResort (sparse);
     for (int i = 0; i < sparse->Nelem; i++) {
 	fprintf (stderr, "%d %d %f\n", 
@@ -33,8 +23,9 @@
     psVector *x = psVectorAlloc (3, PS_DATA_F32);
     x->data.F32[0] = 3;
-    x->data.F32[0] = 5;
-    x->data.F32[0] = 7;
-
-    psVector B = psSparseMatrixTimesVector (NULL, sparse, x);
+    x->data.F32[1] = 5;
+    x->data.F32[2] = 7;
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    psVector *B = psSparseMatrixTimesVector (NULL, sparse, x);
     fprintf (stderr, "B: %f %f %f\n", B->data.F32[0], B->data.F32[1], B->data.F32[2]);
 
@@ -43,5 +34,17 @@
     sparse->Bfj->data.F32[2] = B->data.F32[2];
 
-    x = psSparseSolve (x, sparse);
+    x = psSparseSolve (x, sparse, 0);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, sparse, 1);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, sparse, 2);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, sparse, 3);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, sparse, 4);
     fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     return;
@@ -50,10 +53,10 @@
 void psSparseResort (psSparse *sparse) {
 
-    Nelem = sparse->Nelem;
+    int Nelem = sparse->Nelem;
 
     psVector *index = psVectorSortIndex (NULL, sparse->Sj);
-    psVector Aij = sparse->Aij;
-    psVector Si = sparse->Si;
-    psVector Sj = sparse->Sj;
+    psVector *Aij = sparse->Aij;
+    psVector *Si = sparse->Si;
+    psVector *Sj = sparse->Sj;
 
     // allocate new temporary vectors
@@ -61,6 +64,6 @@
     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];
+    for (int i = 0; i < Nelem; i++) {
+	int j = index->data.U32[i];
 	tAij->data.F32[i] = Aij->data.F32[j];
 	tSi->data.S32[i]  = Si->data.S32[j];
@@ -80,4 +83,6 @@
 void psSparseMatrixElement (psSparse *sparse, int i, int j, float value) {
 
+    int k;
+
     if (i < j) {
 	fprintf (stderr, "*** error: subdiagonal element ***\n");
@@ -91,10 +96,7 @@
 	// check vectors lengths and extend if needed
 	if (sparse->Nelem >= sparse->Aij->nalloc) {
-	    sparse->Aij->nalloc += 100;	    
-	    psVectorRealloc (sparse->Aij, sparse->Aij->nalloc);
-	    sparse->Si->nalloc += 100;	    
-	    psVectorRealloc (sparse->Si, sparse->Si->nalloc);
-	    sparse->Sj->nalloc += 100;	    
-	    psVectorRealloc (sparse->Sj, sparse->Sj->nalloc);
+	    psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100);
+	    psVectorRealloc (sparse->Si,  sparse->Si->nalloc + 100);
+	    psVectorRealloc (sparse->Sj,  sparse->Sj->nalloc + 100);
 	}
 
@@ -111,10 +113,7 @@
 	// check vectors lengths and extend if needed
 	if (sparse->Nelem >= sparse->Aij->nalloc - 1) {
-	    sparse->Aij->nalloc += 100;	    
-	    psVectorRealloc (sparse->Aij, sparse->Aij->nalloc);
-	    sparse->Si->nalloc += 100;	    
-	    psVectorRealloc (sparse->Si, sparse->Si->nalloc);
-	    sparse->Sj->nalloc += 100;	    
-	    psVectorRealloc (sparse->Sj, sparse->Sj->nalloc);
+	    psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100);
+	    psVectorRealloc (sparse->Si,  sparse->Si->nalloc + 100);
+	    psVectorRealloc (sparse->Sj,  sparse->Sj->nalloc + 100);
 	}
 
@@ -139,5 +138,5 @@
 void psSparseVectorElement (psSparse *sparse, int i, float value) {
 
-    sparse->Bjf->data.F32[i] = value;
+    sparse->Bfj->data.F32[i] = value;
     return;
 }
@@ -146,6 +145,9 @@
 psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector) {
 
+    int i, Nelem;
+    float F;
+
     if (output == NULL) {
-	psVector *output = psVectorAlloc (vector->n, PS_DATA_F32);
+	output = psVectorAlloc (vector->n, PS_DATA_F32);
     }
 
@@ -154,6 +156,7 @@
 	F = 0;
 	while (matrix->Sj->data.S32[Nelem] == j) {
-	    i = matrix->Sj->data.S32[Nelem];
-	    F += vector->data.F32[i] * matrix->Aij->data.F32[i];
+	    i = matrix->Si->data.S32[Nelem];
+	    F += vector->data.F32[i] * matrix->Aij->data.F32[Nelem];
+	    Nelem++;
 	}
 	output->data.F32[j] = F;
@@ -162,5 +165,5 @@
 }
 
-psVector *psSparseSolve (psVector *guess, psSparse *sparse) {
+psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter) {
 
     psF32 dG;
@@ -169,12 +172,12 @@
     psVector *Bfj = sparse->Bfj;
 
-    guess = psVectorCopy (guess, Bfj);
+    guess = psVectorCopy (guess, Bfj, PS_DATA_F32);
 
     // temporary storage for intermediate results
     psVector *dQ = psVectorAlloc (guess->n, PS_DATA_F32);
 
-    for (int j = 0; j < 2; j++) {
+    for (int j = 0; j < Niter; j++) {
 	dQ = psSparseMatrixTimesVector (dQ, sparse, guess);
-	for (int i = 0; i < dG->n; i++) {
+	for (int i = 0; i < dQ->n; i++) {
 	    dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
 	    guess->data.F32[i] -= dG;
@@ -187,9 +190,9 @@
 static void psSparseFree (psSparse *sparse) {
     if (sparse == NULL) return;
-    psFree (sparse.Aij);
-    psFree (sparse.Bfj);
-    psFree (sparse.Qii);
-    psFree (sparse.Si);
-    psFree (sparse.Sj);
+    psFree (sparse->Aij);
+    psFree (sparse->Bfj);
+    psFree (sparse->Qii);
+    psFree (sparse->Si);
+    psFree (sparse->Sj);
     return;
 }
@@ -212,4 +215,5 @@
 
     sparse->Nrows = Nrows;
+    psMemSetDeallocator(sparse, (psFreeFunc) psSparseFree);
     return (sparse);
 }
Index: /trunk/psphot/src/psSparse.h
===================================================================
--- /trunk/psphot/src/psSparse.h	(revision 5654)
+++ /trunk/psphot/src/psSparse.h	(revision 5654)
@@ -0,0 +1,18 @@
+
+typedef struct {
+    psVector *Aij;
+    psVector *Bfj;
+    psVector *Qii;
+    psVector *Si;
+    psVector *Sj;
+    int Nelem;
+    int Nrows;
+} psSparse;
+
+void psSparseMatrixTest ();
+void psSparseResort (psSparse *sparse);
+void psSparseMatrixElement (psSparse *sparse, int i, int j, float value);
+void psSparseVectorElement (psSparse *sparse, int i, float value);
+psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector);
+psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter);
+psSparse *psSparseAlloc (int Nrows, int Nelem);
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 5653)
+++ /trunk/psphot/src/psphot.h	(revision 5654)
@@ -9,4 +9,5 @@
 # include "pmPSFtry.h"
 # include "pmModelGroup.h"
+# include "psSparse.h"
 
 typedef struct {
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);
 }
Index: /trunk/psphot/src/psphotTest.c
===================================================================
--- /trunk/psphot/src/psphotTest.c	(revision 5653)
+++ /trunk/psphot/src/psphotTest.c	(revision 5654)
@@ -3,5 +3,5 @@
 int main (int argc, char **argv) {
 
-  psphotSparseMatrixTest ();
+  psSparseMatrixTest ();
 
 }