Index: trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- trunk/psphot/src/psphotEnsemblePSF.c	(revision 10104)
+++ trunk/psphot/src/psphotEnsemblePSF.c	(revision 10171)
@@ -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, bool constant_weights);
+static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool constant_weights, int SKY_FIT_ORDER);
 
 bool psphotEnsemblePSF (pmReadout *readout, psArray *refSources, psMetadata *recipe, pmPSF *psf, bool final) {
@@ -19,5 +19,5 @@
     float y;
     float f;
-    float r;
+    // float r;
 
     psTimerStart ("psphot");
@@ -36,8 +36,16 @@
     if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
 
-    const bool CONSTANT_PHOTOMETRIC_WEIGHTS =
+    bool CONSTANT_PHOTOMETRIC_WEIGHTS =
         psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
     if (!status) {
         psAbort(PS_FILE_LINE, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
+    }
+    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
+    if (!status) {
+	SKY_FIT_ORDER = 0;
+    }
+    bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR");
+    if (!status) {
+	SKY_FIT_LINEAR = false;
     }
 
@@ -100,5 +108,5 @@
 
     // vectors to store stats for each object
-    psVector *weight = psVectorAlloc (fitSources->n, PS_TYPE_F32);
+    // psVector *weight = psVectorAlloc (fitSources->n, PS_TYPE_F32);
     psVector *errors = psVectorAlloc (fitSources->n, PS_TYPE_F32);
 
@@ -106,5 +114,6 @@
     // for just sky: 1 row; for x,y terms: 3 rows
     psSparse *sparse = psSparseAlloc (fitSources->n, 100);
-    psSparseBorder *border = psSparseBorderAlloc (sparse, 1);
+    int nBorder = (SKY_FIT_ORDER == 0) ? 1 : 3;
+    psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
 
     // fill out the sparse matrix elements and border elements (B)
@@ -116,7 +125,9 @@
         pmSource *Fi = fitSources->data[i];
 
-        // scale by diagonal element (auto-cross-product)
-        r = pmSourceCrossProduct (Fi, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS);
-        weight->data.F32[i] = r;
+        // diagonal elements of the sparse matrix (auto-cross-product)
+        f = pmSourceCrossProduct (Fi, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        psSparseMatrixElement (sparse, i, i, f);
+
+        // XXX being used? weight->data.F32[i] = r;
 
 	// the formal error depends on the weighting scheme
@@ -125,20 +136,29 @@
 	    errors->data.F32[i] = 1.0 / sqrt(var);
 	} else {
-	    errors->data.F32[i] = 1.0 / sqrt(r);
+	    errors->data.F32[i] = 1.0 / sqrt(f);
 	}
 
-        psSparseMatrixElement (sparse, i, i, 1.0);
 
         // find the image x model value
         f = pmSourceCrossProduct (Ri, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS);
-        psSparseVectorElement (sparse, i, f / r);
+        psSparseVectorElement (sparse, i, f);
 
 	// add the per-source weights (border region)
-	float p = pmSourceWeight (Fi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
-	// float px = pmSourceWeight (Fi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
-	// float py = pmSourceWeight (Fi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
-	psSparseBorderElementB (border, i, 0, p);
-	// psSparseBorderElementB (border, i, 1, px);
-	// psSparseBorderElementB (border, i, 2, py);
+	switch (SKY_FIT_ORDER) {
+	  case 1:
+	    f = pmSourceWeight (Fi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
+	    psSparseBorderElementB (border, i, 1, f);
+	    f = pmSourceWeight (Fi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
+	    psSparseBorderElementB (border, i, 2, f);
+
+	  case 0:
+	    f = pmSourceWeight (Fi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
+	    psSparseBorderElementB (border, i, 0, f);
+	    break;
+
+	  default:
+	    psAbort(PS_FILE_LINE, "Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER);
+	    break;
+	}
 
         // loop over all other stars following this one
@@ -154,5 +174,5 @@
             // got an overlap; calculate cross-product and add to output array
             f = pmSourceCrossProduct (Fi, Fj, CONSTANT_PHOTOMETRIC_WEIGHTS);
-            psSparseMatrixElement (sparse, j, i, f / r);
+            psSparseMatrixElement (sparse, j, i, f);
         }
     }
@@ -161,6 +181,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, CONSTANT_PHOTOMETRIC_WEIGHTS);
+    SetBorderT (border, readout, refSources, fitSources, index, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER);
 
     psSparseConstraint constraint;
@@ -170,13 +189,13 @@
 
     // solve for normalization terms (need include local sky?)
-    # 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);
-    psVector *skyfit = NULL;
-    # endif
+    if (SKY_FIT_LINEAR) {
+	psSparseBorderSolve (&norm, &skyfit, constraint, border, 5);
+	fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]);
+    } else {
+	norm = psSparseSolve (NULL, constraint, sparse, 5);
+	skyfit = NULL;
+    }
 
     // adjust fitSources, set refSources and subtract
@@ -220,5 +239,4 @@
     psFree (norm);
     psFree (skyfit);
-    psFree (weight);
     psFree (errors);
     psFree (border);
@@ -229,5 +247,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, bool constant_weights) {
+static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, psArray *fitSources, psVector *index, bool constant_weights, int SKY_FIT_ORDER) {
 
     // generate the image-wide weight terms
@@ -271,6 +289,6 @@
 	    w   += 1/wt;
 	    fo  += f/wt;
-
-	    if (!fitSlope) continue;
+	    if (SKY_FIT_ORDER == 0) continue;
+
 	    xc = i + col0;
 	    yc = j + row0;
@@ -291,5 +309,5 @@
     psSparseBorderElementG (border, 0, fo);
     psSparseBorderElementT (border, 0, 0, w);
-    if (fitSlope) {
+    if (SKY_FIT_ORDER > 0) {
 	psSparseBorderElementG (border, 0, fx);
 	psSparseBorderElementG (border, 0, fy);