Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 10075)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 10076)
@@ -97,5 +97,6 @@
 
     // XXX test dump of psf stars and model
-    if (0) { 
+    if (1) { 
+	psphotSaveImage (NULL, readout->image,  "testsub.fits");
 	pmSourcesWritePSFs (stars, "psfstars.dat");
 	try = models->data[0];
Index: /trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- /trunk/psphot/src/psphotEnsemblePSF.c	(revision 10075)
+++ /trunk/psphot/src/psphotEnsemblePSF.c	(revision 10076)
@@ -10,5 +10,5 @@
 // the analysis is performed wrt the simulated pixel values
 
-static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, bool fitSlope);
+static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope);
 
 bool psphotEnsemblePSF (pmReadout *readout, psArray *refSources, psMetadata *recipe, pmPSF *psf, bool final) {
@@ -162,5 +162,5 @@
     // set the sky, sky_x, sky_y components of border matrix
     // XXX ignore sky slope terms for now
-    SetBorderT (border, readout, refSources, false);
+    SetBorderT (border, readout, refSources, fitSources, index, false);
 
     psSparseConstraint constraint;
@@ -170,7 +170,11 @@
 
     // solve for normalization terms (need include local sky?)
+    # if (0)
     psVector *norm = NULL;
     psVector *skyfit = NULL;
     psSparseBorderSolve (&norm, &skyfit, constraint, border, 3);
+    # else
+    psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
+    # endif
 
     // adjust fitSources, set refSources and subtract
@@ -215,4 +219,5 @@
     psFree (weight);
     psFree (errors);
+    psFree (border);
 
     psLogMsg (__func__, 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot"));
@@ -221,5 +226,5 @@
 
 // calculate the weight terms for the sky fit component of the matrix
-static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, bool fitSlope) {
+static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope) {
 
     // generate the image-wide weight terms
@@ -231,7 +236,9 @@
     // turn off MARK for all object pixels
     // we must use the refSources here since their pixels point at the science image
-    for (int i = 0; i < refSources->n; i++) {
-        pmSource *refSource = refSources->data[i];
+    for (int i = 0; i < fitSources->n; i++) {
+        int N = index->data.U32[i];
+        pmSource *refSource = refSources->data[N];
 	pmModel *model = pmSourceGetModel (NULL, refSource);
+	if (model == NULL) continue;
         float x = model->params->data.F32[PM_PAR_XPOS];
         float y = model->params->data.F32[PM_PAR_YPOS];