Index: /trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- /trunk/psphot/src/psphotEnsemblePSF.c	(revision 10031)
+++ /trunk/psphot/src/psphotEnsemblePSF.c	(revision 10032)
@@ -1,6 +1,4 @@
 # include "psphot.h"
-bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight);
-
-// 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) {
@@ -29,5 +27,4 @@
     psArray  *models = psArrayAllocEmpty (sources->n);
     psVector *index  = psVectorAllocEmpty (sources->n, PS_TYPE_U32);
-    // DROP models->n = index->n = 0;
 
     // option to limit analysis to a specific region
@@ -66,4 +63,5 @@
         // really saturated stars should be re-measured for a better centroid
         // XXX EAM : move this to a 'clear satstar function'
+	// XXX : place this in pmSourceMoments or wherever it is called?
         if (inSource->mode &  PM_SOURCE_MODE_SATSTAR) {
             status = pmSourceMoments (inSource, INNER_RADIUS);
@@ -109,4 +107,5 @@
 
         // make temporary copies of the image pixels and mask
+	// we need to have a copy which will not be modified by changes to its neighbor
         otSource->mask   = psImageCopy (NULL, inSource->mask,   PS_TYPE_U8);
         otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32);
@@ -142,4 +141,8 @@
     psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32);
     psVector *errors = psVectorAlloc (models->n, PS_TYPE_F32);
+
+    // create the border matrix (includes the sparse matrix)
+    // for just sky: 1 row; for x,y terms: 3 rows
+    psSparseBorder *border = psSparseBorderAlloc (sparse, 1);
 
     for (int i = 0; i < models->n; i++) {
@@ -166,4 +169,12 @@
         psSparseVectorElement (sparse, i, f / r);
 
+	// add the per-source weights (border region)
+	p = pmSourceWeight (Mi, CONST, CONSTANT_PHOTOMETRIC_WEIGHTS);
+	// px = pmSourceWeight (Mi, XTERM, CONSTANT_PHOTOMETRIC_WEIGHTS);
+	// py = pmSourceWeight (Mi, YTERM, CONSTANT_PHOTOMETRIC_WEIGHTS);
+	psSparseBorderElementB (border, i, 0, p);
+	// psSparseBorderElementB (border, i, 1, px);
+	// psSparseBorderElementB (border, i, 2, py);
+
         // loop over all other stars following this one
         for (int j = i + 1; j < models->n; j++) {
@@ -183,4 +194,41 @@
     psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
 
+    // generate the image-wide weight terms
+    // first, mask all object pixels
+    for (int i = 0; i < models->n; i++) {
+        int N = index->data.U32[i];
+        pmSource *Fi = sources->data[N];
+        pmSource *Mi = models->data[i];
+
+	// XXX need to select the correct mode (check this against the selections below
+	// - which model is appropriate for the given source?
+	// - which radius is appropriate for the given source?
+	// - where do x,y come from?
+	psImageKeepCircle (Mi->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
+    }	
+
+    // accumulate the image statistics from the masked regions
+    for (int j = 0; j < numRows; j++) {
+	for (int i = 0; i < numCols; i++) {
+	    w  +=     image->data.F32[j][i]/weight->data.F32[j][i];
+	    x  +=   i*image->data.F32[j][i]/weight->data.F32[j][i];
+	    y  +=   j*image->data.F32[j][i]/weight->data.F32[j][i];
+	    x2 += i*i*image->data.F32[j][i]/weight->data.F32[j][i];
+	    xy += i*j*image->data.F32[j][i]/weight->data.F32[j][i];
+	    y2 += j*j*image->data.F32[j][i]/weight->data.F32[j][i];
+	}
+    }
+    psSparseBorderElementT (border, 0, 0, w);
+    if (0) {
+	psSparseBorderElementT (border, 1, 0, x);
+	psSparseBorderElementT (border, 2, 0, y);
+	psSparseBorderElementT (border, 0, 1, x);
+	psSparseBorderElementT (border, 1, 1, x2);
+	psSparseBorderElementT (border, 2, 1, xy);
+	psSparseBorderElementT (border, 0, 2, y);
+	psSparseBorderElementT (border, 1, 2, xy);
+	psSparseBorderElementT (border, 2, 2, y2);
+    }    
+
     psSparseConstraint constraint;
     constraint.paramMin   = 0.0;
@@ -190,5 +238,7 @@
     // solve for normalization terms (need include local sky?)
     psSparseResort (sparse);
-    psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
+    psSparseBorderSolve (&xFit, &yFit, constraint, border, 3);
+
+    // psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
 
     // adjust models, set sources and subtract
@@ -244,20 +294,2 @@
     return true;
 }
-
-bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight) {
-
-    double dC = 0.0;
-    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 ++;
-        }
-    }
-    model->nDOF = Npix - 1;
-    model->chisq = dC;
-
-    return (true);
-}