Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 10095)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 10096)
@@ -26,4 +26,10 @@
     }
 
+    // supply the measured sky variance for optional constant errors (non-poissonian)
+    float SKY_STDEV = psMetadataLookupF32 (&status, recipe, "SKY_STDEV");
+    if (!status) {
+	SKY_STDEV = 1.0;
+        psWarning("SKY_STDEV is not set --- defaulting to %f\n", SKY_STDEV);
+    }
     // use poissonian errors or local-sky errors
     bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
@@ -32,4 +38,5 @@
         psWarning("POISSON_ERRORS is not set in the recipe --- defaulting to true.\n");
     }
+    pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_STDEV), POISSON_ERRORS);
 
     // how to model the PSF variations across the field
@@ -48,6 +55,4 @@
     }
 
-    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);
-
     stars = psArrayAllocEmpty (sources->n);
 
@@ -94,15 +99,4 @@
         modelName = item->data.V;
         models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask);
-    }
-
-    // XXX test dump of psf stars and model
-    if (1) { 
-	psphotSaveImage (NULL, readout->image,  "testsub.fits");
-	pmSourcesWritePSFs (stars, "psfstars.dat");
-	try = models->data[0];
-	psf = try->psf;
-	psMetadata *psfData = pmPSFtoMetadata (NULL, psf);
-        psMetadataConfigWrite (psfData, "psfmodel.dat");
-	psFree (psfData);
     }
 
@@ -138,4 +132,30 @@
     try = models->data[bestN];
 
+    // XXX test dump of psf star data and psf-subtracted image
+    if (0) { 
+	for (int i = 0; i < try->sources->n; i++) {
+	    // masked for: bad model fit, outlier in parameters
+	    if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
+		continue;
+
+	    pmSource *source = try->sources->data[i];
+	    float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+	    float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+
+	    // set the mask and subtract the PSF model
+	    psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
+	    pmModelSub (source->pixels, source->mask, source->modelPSF, false, false);
+	    psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
+	}
+
+	psphotSaveImage (NULL, readout->image,  "psfstars.fits");
+	pmSourcesWritePSFs (try->sources, "psfstars.dat");
+	psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
+	psMetadataConfigWrite (psfData, "psfmodel.dat");
+	psFree (psfData);
+	psLogMsg ("psphot.choosePSF", 3, "wrote out psf-subtracted image, psf data, exiting\n");
+	exit (0);
+    }
+
     // unset the PSFSTAR flag for stars not used for PSF model
     for (int i = 0; i < try->sources->n; i++) {
Index: /trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- /trunk/psphot/src/psphotEnsemblePSF.c	(revision 10095)
+++ /trunk/psphot/src/psphotEnsemblePSF.c	(revision 10096)
@@ -10,5 +10,5 @@
 // the analysis is performed wrt the simulated pixel values
 
-static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope);
+static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope, bool constant_weights);
 
 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, fitSources, index, false);
+    SetBorderT (border, readout, refSources, fitSources, index, false, CONSTANT_PHOTOMETRIC_WEIGHTS);
 
     psSparseConstraint constraint;
@@ -170,8 +170,9 @@
 
     // solve for normalization terms (need include local sky?)
-    # if (0)
+    # if (1)
     psVector *norm = NULL;
     psVector *skyfit = NULL;
     psSparseBorderSolve (&norm, &skyfit, constraint, border, 3);
+    fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]);
     # else
     psVector *norm = psSparseSolve (NULL, constraint, sparse, 3);
@@ -226,5 +227,5 @@
 
 // calculate the weight terms for the sky fit component of the matrix
-static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope) {
+static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool fitSlope, bool constant_weights) {
 
     // generate the image-wide weight terms
@@ -251,6 +252,6 @@
     psU8 **mask = readout->mask->data.U8;
 
-    double w, x, y, x2, xy, y2, xc, yc;
-    w = x = y = x2 = xy = y2 = 0;
+    double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy;
+    w = x = y = x2 = xy = y2 = fo = fx = fy = 0;
     
     int col0 = readout->image->col0;
@@ -260,13 +261,23 @@
 	for (int i = 0; i < readout->image->numCols; i++) {
 	    if (mask[j][i]) continue;
-	    w  +=       image[j][i]/weight[j][i];
+	    if (constant_weights) {
+		wt = 1.0;
+	    } else {
+		wt = weight[j][i];
+	    }
+	    f = image[j][i];
+	    w   += 1/wt;
+	    fo  += f/wt;
+
 	    if (!fitSlope) continue;
 	    xc = i + col0;
 	    yc = j + row0;
-	    x  +=    xc*image[j][i]/weight[j][i];
-	    y  +=    yc*image[j][i]/weight[j][i];
-	    x2 += xc*xc*image[j][i]/weight[j][i];
-	    xy += xc*yc*image[j][i]/weight[j][i];
-	    y2 += yc*yc*image[j][i]/weight[j][i];
+	    x  +=    xc/wt;
+	    y  +=    yc/wt;
+	    x2 += xc*xc/wt;
+	    xy += xc*yc/wt;
+	    y2 += yc*yc/wt;
+	    fx +=  f*xc/wt;
+	    fy +=  f*yc/wt;
 	}
     }
@@ -276,6 +287,9 @@
 
     // set the Border T elements
+    psSparseBorderElementG (border, 0, fo);
     psSparseBorderElementT (border, 0, 0, w);
     if (fitSlope) {
+	psSparseBorderElementG (border, 0, fx);
+	psSparseBorderElementG (border, 0, fy);
 	psSparseBorderElementT (border, 1, 0, x);
 	psSparseBorderElementT (border, 2, 0, y);
Index: /trunk/psphot/src/psphotImageMedian.c
===================================================================
--- /trunk/psphot/src/psphotImageMedian.c	(revision 10095)
+++ /trunk/psphot/src/psphotImageMedian.c	(revision 10096)
@@ -124,5 +124,5 @@
     stats = psImageStats (stats, model->image, NULL, 0);
     psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MEAN", PS_DATA_F32 | PS_META_REPLACE, "sky model mean",          stats->sampleMean);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_STDV", PS_DATA_F32 | PS_META_REPLACE, "sky model stdev",         stats->sampleStdev);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_STDEV", PS_DATA_F32 | PS_META_REPLACE, "sky model stdev",         stats->sampleStdev);
     psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MAX",  PS_DATA_F32 | PS_META_REPLACE, "sky model maximum value", stats->max);
     psMetadataAdd (recipe, PS_LIST_TAIL, "SKY_MIN",  PS_DATA_F32 | PS_META_REPLACE, "sky model minimum value", stats->min);
Index: /trunk/psphot/src/psphotModelTest.c
===================================================================
--- /trunk/psphot/src/psphotModelTest.c	(revision 10095)
+++ /trunk/psphot/src/psphotModelTest.c	(revision 10096)
@@ -17,5 +17,5 @@
     bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
     if (!status) POISSON_ERRORS = true;
-    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);
+    pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS);
 
     // run model fitting tests on a single source
Index: /trunk/psphot/src/psphotOutput.c
===================================================================
--- /trunk/psphot/src/psphotOutput.c	(revision 10095)
+++ /trunk/psphot/src/psphotOutput.c	(revision 10096)
@@ -80,5 +80,5 @@
     // sky background model statistics
     psMetadataItemTransfer (header, recipe, "SKY_MEAN");
-    psMetadataItemTransfer (header, recipe, "SKY_STDV");
+    psMetadataItemTransfer (header, recipe, "SKY_STDEV");
     psMetadataItemTransfer (header, recipe, "SKY_MIN");
     psMetadataItemTransfer (header, recipe, "SKY_MAX");
Index: /trunk/psphot/src/psphotTestPSF.c
===================================================================
--- /trunk/psphot/src/psphotTestPSF.c	(revision 10095)
+++ /trunk/psphot/src/psphotTestPSF.c	(revision 10096)
@@ -32,4 +32,5 @@
         psWarning("POISSON_ERRORS is not set in the recipe --- defaulting to true.\n");
     }
+    pmSourceFitModelInit (15, 0.1, 1.0, POISSON_ERRORS);
 
     // how to model the PSF variations across the field
@@ -47,6 +48,4 @@
         }
     }
-
-    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);
 
     stars = psArrayAllocEmpty (sources->n);
Index: /trunk/psphot/src/psphotWeightBias.c
===================================================================
--- /trunk/psphot/src/psphotWeightBias.c	(revision 10095)
+++ /trunk/psphot/src/psphotWeightBias.c	(revision 10096)
@@ -18,5 +18,16 @@
 
     // set fitting method to use non-poisson errors (local sky error)
-    pmSourceFitModelInit (15, 0.1, false);
+    float SKY_STDEV = psMetadataLookupF32 (&status, recipe, "SKY_STDEV");
+    if (!status) {
+      SKY_STDEV = 1.0;
+        psWarning("SKY_STDEV is not set --- defaulting to %f\n", SKY_STDEV);
+    }
+    // use poissonian errors or local-sky errors
+    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
+    if (!status) {
+        POISSON_ERRORS = true;
+        psWarning("POISSON_ERRORS is not set in the recipe --- defaulting to true.\n");
+    }
+    pmSourceFitModelInit (15, 0.1, PS_SQR(SKY_STDEV), POISSON_ERRORS);
 
     // option to limit analysis to a specific region