Index: /branches/eam_branches/ipp-20120601/Ohana/src/libkapa/src/DrawRotString.c
===================================================================
--- /branches/eam_branches/ipp-20120601/Ohana/src/libkapa/src/DrawRotString.c	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/Ohana/src/libkapa/src/DrawRotString.c	(revision 34053)
@@ -77,4 +77,5 @@
     /* check for special characters */
     if (!code && !protect) {
+      /* subscript, starts with _ */
       if (N == 94) {
 	SetRotFont (currentname, (int)(0.8*currentsize));
@@ -83,4 +84,5 @@
 	continue;
       }
+      /* superscript, starts with ^ */
       if (N == 95) { 
 	SetRotFont (currentname, (int)(0.8*currentsize));
@@ -89,4 +91,5 @@
 	continue;
       }
+      /* normal script, starts with | */
       if (N == 124) {
 	SetRotFont (currentname, basesize);
Index: /branches/eam_branches/ipp-20120601/psLib/src/astro/psTime.c
===================================================================
--- /branches/eam_branches/ipp-20120601/psLib/src/astro/psTime.c	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psLib/src/astro/psTime.c	(revision 34053)
@@ -389,5 +389,5 @@
         if (i < numTables) {
             table = psLookupTableAlloc(fullTableName, (const char*)tableFormat, tablesIndex->data.S32[i]);
-            // XXX unused numLines = psLookupTableRead(table);
+            psLookupTableRead(table);
         } else {
             psError(PS_ERR_BAD_PARAMETER_VALUE, no_problem,
Index: /branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelFuncs.h
===================================================================
--- /branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelFuncs.h	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelFuncs.h	(revision 34053)
@@ -29,6 +29,4 @@
 # define PM_MODEL_FUNCS_H
 
-# define HAVE_MODEL_VAR 1
-
 /// @addtogroup Objects Object Detection / Analysis Functions
 /// @{
@@ -57,7 +55,5 @@
     PM_MODEL_OP_NORM     = 0x20,
     PM_MODEL_OP_NOISE    = 0x40,
-# if (HAVE_MODEL_VAR)
     PM_MODEL_OP_MODELVAR = 0x80,
-# endif
 } pmModelOpMode;
 
Index: /branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelUtils.c
===================================================================
--- /branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelUtils.c	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelUtils.c	(revision 34053)
@@ -43,4 +43,13 @@
 
 #include "pmErrorCodes.h"
+
+// XX static bool useModelVar = false;
+// XX 
+// XX void pmModelSetModelVarOption (bool option) {
+// XX   useModelVar = option;
+// XX }
+// XX bool pmModelGetModelVarOption (void) {
+// XX   return useModelVar;
+// XX }
 
 /*****************************************************************************
Index: /branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelUtils.h
===================================================================
--- /branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelUtils.h	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psModules/src/objects/pmModelUtils.h	(revision 34053)
@@ -46,4 +46,7 @@
 bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments);
 
+// XXX void pmModelSetModelVarOption (bool option);
+// XXX bool pmModelGetModelVarOption (void);
+
 /// @}
 # endif /* PM_MODEL_UTILS_H */
Index: /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSource.c
===================================================================
--- /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSource.c	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSource.c	(revision 34053)
@@ -51,7 +51,5 @@
     psFree(tmp->pixels);
     psFree(tmp->variance);
-# if (HAVE_MODEL_VAR)
     psFree(tmp->modelVar);
-# endif
     psFree(tmp->maskObj);
     psFree(tmp->maskView);
@@ -80,7 +78,5 @@
     psFree (source->pixels);
     psFree (source->variance);
-# if (HAVE_MODEL_VAR)
     psFree (source->modelVar);
-# endif
     psFree (source->maskObj);
     psFree (source->maskView);
@@ -90,7 +86,5 @@
     source->pixels = NULL;
     source->variance = NULL;
-# if (HAVE_MODEL_VAR)
     source->modelVar = NULL;
-# endif
     source->maskObj = NULL;
     source->maskView = NULL;
@@ -122,7 +116,5 @@
     source->pixels = NULL;
     source->variance = NULL;
-# if (HAVE_MODEL_VAR)
     source->modelVar = NULL;
-# endif
     source->maskObj = NULL;
     source->maskView = NULL;
@@ -212,7 +204,5 @@
     source->pixels   = in->pixels   ? psImageCopyView(NULL, in->pixels)   : NULL;
     source->variance = in->variance ? psImageCopyView(NULL, in->variance) : NULL;
-# if (HAVE_MODEL_VAR)
     source->modelVar = NULL;
-# endif
     source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL;
 
@@ -1044,8 +1034,6 @@
     bool addNoise = mode & PM_MODEL_OP_NOISE;
 
-# if (HAVE_MODEL_VAR)
     bool addModelVar = mode & PM_MODEL_OP_MODELVAR;
     if (addModelVar) psAssert (source->modelVar, "programming error");
-# endif
 
     // require the use of pmModelAddWithOffset if we are adding noise (because the model size and norm are rescaled)
@@ -1070,9 +1058,6 @@
         }
 
-# if (HAVE_MODEL_VAR)
-        psF32 **target = addModelVar ? source->modelVar->data.F32 : source->pixels->data.F32;
-# else
-        psF32 **target = source->pixels->data.F32;
-# endif
+	psF32 **target = NULL;
+	target = addModelVar ? source->modelVar->data.F32 : source->pixels->data.F32;
 
         for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
@@ -1089,19 +1074,12 @@
             }
         }
-# if (HAVE_MODEL_VAR)
+	// do not change the flag here if we are adding/subtracting from modelVar
 	if (!addModelVar) {
-	    if (add) {
-		source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
-	    } else {
-		source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
-	    }
+	  if (add) {
+	    source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+	  } else {
+	    source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+	  }
 	}
-# else
-	if (add) {
-	    source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
-	} else {
-	    source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
-	}
-# endif
         return true;
     }
@@ -1111,24 +1089,14 @@
         target = source->variance;
     }
-# if (HAVE_MODEL_VAR)
     if (addModelVar) {
-        target = source->modelVar;
-    }
-# endif
+      target = source->modelVar;
+    }
 
     if (add) {
 	status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
-# if (HAVE_MODEL_VAR)
 	if (!addNoise && !addModelVar) source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
-# else
-	source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
-# endif
     } else {
 	status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
-# if (HAVE_MODEL_VAR)
 	if (!addNoise && !addModelVar) source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
-# else
-	source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
-# endif
     }
     if (!status) {
Index: /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSource.h
===================================================================
--- /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSource.h	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSource.h	(revision 34053)
@@ -72,7 +72,5 @@
     psImage *pixels;                    ///< Rectangular region including object pixels.
     psImage *variance;			///< Image variance.
-# if (HAVE_MODEL_VAR)
     psImage *modelVar;			///< variance based on current models
-# endif
     psImage *maskObj;                   ///< unique mask for this object which marks included pixels associated with objects.
     psImage *maskView;                  ///< view into global image mask for this object region
Index: /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSourcePhotometry.c
===================================================================
--- /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSourcePhotometry.c	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSourcePhotometry.c	(revision 34053)
@@ -899,21 +899,14 @@
 }
 
-# if (HAVE_MODEL_VAR) 
 double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
-# else
-double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
-# endif
 {
     PS_ASSERT_PTR_NON_NULL(Mi, NAN);
-# if (HAVE_MODEL_VAR) 
-    double flux = 0, wt = 1.0, factor = 0;
-# else
-    double flux = 0, wt = 0, factor = 0;
-# endif
+    double flux = 0;
+    double wt = 1.0;
+    double factor = 0;
 
     const psImage *Pi = Mi->modelFlux;
     assert (Pi != NULL);
 
-# if (HAVE_MODEL_VAR)
     const psImage *Wi = NULL;
     switch (fitVarMode) {
@@ -921,4 +914,152 @@
 	break;
       case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+      case PM_SOURCE_PHOTFIT_MODEL_SKY:
+	Wi = Mi->variance;
+	psAssert (Wi, "programming error");
+	break;
+      case PM_SOURCE_PHOTFIT_MODEL_VAR:
+	Wi = Mi->modelVar;
+	psAssert (Wi, "programming error");
+	break;
+      case PM_SOURCE_PHOTFIT_NONE:
+	psAbort("programming error");
+    }	
+    const psImage *Ti = Mi->maskObj;
+    assert (Ti != NULL);
+
+    for (int yi = 0; yi < Pi->numRows; yi++) {
+	for (int xi = 0; xi < Pi->numCols; xi++) {
+	    if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
+		continue;
+	    if (fitVarMode != PM_SOURCE_PHOTFIT_CONST) {
+		wt = covarFactor * Wi->data.F32[yi][xi];
+		if (wt == 0) continue;
+	    }
+	    switch (term) {
+	      case 0:
+		factor = 1;
+		break;
+	      case 1:
+		factor = xi + Pi->col0;
+		break;
+	      case 2:
+		factor = yi + Pi->row0;
+		break;
+	      default:
+		psAbort("invalid term for pmSourceWeight");
+	    }
+
+	    // wt is 1.0 for CONST
+	    flux += (factor * Pi->data.F32[yi][xi]) / wt;
+	}
+    }
+    return flux;
+}
+
+double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
+    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    double flux;
+    double wt = 1.0;
+
+    const psImage *Pi = Mi->modelFlux;
+    assert (Pi != NULL);
+    const psImage *Pj = Mj->modelFlux;
+    assert (Pj != NULL);
+
+    const psImage *Wi = NULL;
+    switch (fitVarMode) {
+      case PM_SOURCE_PHOTFIT_CONST:
+	break;
+      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+      case PM_SOURCE_PHOTFIT_MODEL_SKY:
+	Wi = Mi->variance;
+	psAssert (Wi, "programming error");
+	break;
+      case PM_SOURCE_PHOTFIT_MODEL_VAR:
+	Wi = Mi->modelVar;
+	psAssert (Wi, "programming error");
+	break;
+      case PM_SOURCE_PHOTFIT_NONE:
+	psAbort("programming error");
+    }	
+
+    const psImage *Ti = Mi->maskObj;
+    assert (Ti != NULL);
+    const psImage *Tj = Mj->maskObj;
+    assert (Tj != NULL);
+
+    Xs = PS_MAX (Pi->col0, Pj->col0);
+    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
+
+    Ys = PS_MAX (Pi->row0, Pj->row0);
+    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
+
+    xIs = Xs - Pi->col0;
+    xJs = Xs - Pj->col0;
+    yIs = Ys - Pi->row0;
+    yJs = Ys - Pj->row0;
+
+    xIe = Xe - Pi->col0;
+    yIe = Ye - Pi->row0;
+
+    // note that weight is addressing the same image pixels
+    flux = 0;
+    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
+        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
+            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
+                continue;
+            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
+                continue;
+
+	    float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
+	    switch (fitVarMode) {
+	      case PM_SOURCE_PHOTFIT_CONST:
+		wt = 1.0;
+		break;
+	      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+	      case PM_SOURCE_PHOTFIT_MODEL_SKY:
+	      case PM_SOURCE_PHOTFIT_MODEL_VAR:
+		wt = covarFactor * Wi->data.F32[yi][xi];
+		break;
+	      case PM_SOURCE_PHOTFIT_NONE:
+		psAbort("programming error");
+	    }
+	    // skip pixels with nonsense weight values
+	    if (wt <= 0) continue;
+
+	    flux += value / wt;
+        }
+    }
+    return flux;
+}
+
+double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
+    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    double flux;
+    double wt = 1.0;
+
+    const psImage *Pi = Mi->pixels;
+    assert (Pi != NULL);
+    const psImage *Pj = Mj->modelFlux;
+    assert (Pj != NULL);
+
+    const psImage *Wi = NULL;
+    switch (fitVarMode) {
+      case PM_SOURCE_PHOTFIT_CONST:
+	break;
+      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+      case PM_SOURCE_PHOTFIT_MODEL_SKY:
 	Wi = Mi->variance;
         psAssert (Wi, "programming error");
@@ -931,105 +1072,4 @@
 	psAbort("programming error");
     }	
-# else
-    const psImage *Wi = Mi->variance;
-    if (!unweighted_sum) {
-        assert (Wi != NULL);
-    }
-# endif
-    const psImage *Ti = Mi->maskObj;
-    assert (Ti != NULL);
-
-    for (int yi = 0; yi < Pi->numRows; yi++) {
-        for (int xi = 0; xi < Pi->numCols; xi++) {
-            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
-                continue;
-# if (HAVE_MODEL_VAR)
-            if (fitVarMode != PM_SOURCE_PHOTFIT_CONST) {
-                wt = covarFactor * Wi->data.F32[yi][xi];
-                if (wt == 0) continue;
-            }
-# else
-            if (!unweighted_sum) {
-                wt = covarFactor * Wi->data.F32[yi][xi];
-                if (wt == 0)
-                    continue;
-            }
-# endif
-
-            switch (term) {
-	      case 0:
-                factor = 1;
-                break;
-	      case 1:
-                factor = xi + Pi->col0;
-                break;
-	      case 2:
-                factor = yi + Pi->row0;
-                break;
-	      default:
-                psAbort("invalid term for pmSourceWeight");
-            }
-
-# if (HAVE_MODEL_VAR)
-	    // wt is 1.0 for CONST
-	    flux += (factor * Pi->data.F32[yi][xi]) / wt;
-# else
-            if (unweighted_sum) {
-                flux += (factor * Pi->data.F32[yi][xi]);
-            } else {
-                flux += (factor * Pi->data.F32[yi][xi]) / wt;
-            }
-# endif
-        }
-    }
-    return flux;
-}
-
-# if (HAVE_MODEL_VAR)
-double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
-# else
-double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
-# endif
-{
-    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
-    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
-    int Xs, Xe, Ys, Ye;
-    int xi, xj, yi, yj;
-    int xIs, xJs, yIs, yJs;
-    int xIe, yIe;
-# if (HAVE_MODEL_VAR)
-    double flux;
-    double wt = 1.0;
-# else
-    double flux, wt;
-# endif
-
-    const psImage *Pi = Mi->modelFlux;
-    assert (Pi != NULL);
-    const psImage *Pj = Mj->modelFlux;
-    assert (Pj != NULL);
-
-# if (HAVE_MODEL_VAR)
-    const psImage *Wi = NULL;
-    switch (fitVarMode) {
-      case PM_SOURCE_PHOTFIT_CONST:
-	break;
-      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
-	Wi = Mi->variance;
-        psAssert (Wi, "programming error");
-	break;
-      case PM_SOURCE_PHOTFIT_MODEL_VAR:
-	Wi = Mi->modelVar;
-        psAssert (Wi, "programming error");
-	break;
-      case PM_SOURCE_PHOTFIT_NONE:
-	psAbort("programming error");
-    }	
-# else
-    const psImage *Wi = Mi->variance;
-    if (!unweighted_sum) {
-        assert (Wi != NULL);
-    }
-# endif
 
     const psImage *Ti = Mi->maskObj;
@@ -1052,5 +1092,5 @@
     yIe = Ye - Pi->row0;
 
-    // note that weight is addressing the same image pixels
+    // note that weight is addressing the same image pixels,
     flux = 0;
     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
@@ -1061,5 +1101,4 @@
                 continue;
 
-# if (HAVE_MODEL_VAR)
 	    float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
 	    switch (fitVarMode) {
@@ -1068,4 +1107,5 @@
 		break;
 	      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
+	      case PM_SOURCE_PHOTFIT_MODEL_SKY:
 	      case PM_SOURCE_PHOTFIT_MODEL_VAR:
                 wt = covarFactor * Wi->data.F32[yi][xi];
@@ -1078,126 +1118,6 @@
 
 	    flux += value / wt;
-# else
-            // XXX skip the nonsense weight pixels?
-            if (unweighted_sum) {
-                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
-            } else {
-                wt = covarFactor * Wi->data.F32[yi][xi];
-                if (wt > 0) {
-                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
-                }
-            }
-# endif
         }
     }
     return flux;
 }
-
-# if (HAVE_MODEL_VAR)
-double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal)
-# else
-double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)
-# endif
-{
-    PS_ASSERT_PTR_NON_NULL(Mi, NAN);
-    PS_ASSERT_PTR_NON_NULL(Mj, NAN);
-    int Xs, Xe, Ys, Ye;
-    int xi, xj, yi, yj;
-    int xIs, xJs, yIs, yJs;
-    int xIe, yIe;
-# if (HAVE_MODEL_VAR)
-    double flux;
-    double wt = 1.0;
-# else
-    double flux, wt;
-# endif
-
-    const psImage *Pi = Mi->pixels;
-    assert (Pi != NULL);
-    const psImage *Pj = Mj->modelFlux;
-    assert (Pj != NULL);
-
-# if (HAVE_MODEL_VAR)
-    const psImage *Wi = NULL;
-    switch (fitVarMode) {
-      case PM_SOURCE_PHOTFIT_CONST:
-	break;
-      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
-	Wi = Mi->variance;
-        psAssert (Wi, "programming error");
-	break;
-      case PM_SOURCE_PHOTFIT_MODEL_VAR:
-	Wi = Mi->modelVar;
-        psAssert (Wi, "programming error");
-	break;
-      case PM_SOURCE_PHOTFIT_NONE:
-	psAbort("programming error");
-    }	
-# else
-    const psImage *Wi = Mi->variance;
-    if (!unweighted_sum) {
-        assert (Wi != NULL);
-    }
-# endif
-
-    const psImage *Ti = Mi->maskObj;
-    assert (Ti != NULL);
-    const psImage *Tj = Mj->maskObj;
-    assert (Tj != NULL);
-
-    Xs = PS_MAX (Pi->col0, Pj->col0);
-    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
-
-    Ys = PS_MAX (Pi->row0, Pj->row0);
-    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
-
-    xIs = Xs - Pi->col0;
-    xJs = Xs - Pj->col0;
-    yIs = Ys - Pi->row0;
-    yJs = Ys - Pj->row0;
-
-    xIe = Xe - Pi->col0;
-    yIe = Ye - Pi->row0;
-
-    // note that weight is addressing the same image pixels,
-    flux = 0;
-    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
-        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
-            if (Ti->data.PS_TYPE_IMAGE_MASK_DATA[yi][xi] & maskVal)
-                continue;
-            if (Tj->data.PS_TYPE_IMAGE_MASK_DATA[yj][xj] & maskVal)
-                continue;
-
-# if (HAVE_MODEL_VAR)
-	    float value = (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
-	    switch (fitVarMode) {
-	      case PM_SOURCE_PHOTFIT_CONST:
-		wt = 1.0;
-		break;
-	      case PM_SOURCE_PHOTFIT_IMAGE_VAR:
-	      case PM_SOURCE_PHOTFIT_MODEL_VAR:
-                wt = covarFactor * Wi->data.F32[yi][xi];
-		break;
-	      case PM_SOURCE_PHOTFIT_NONE:
-		psAbort("programming error");
-	    }
-            // skip pixels with nonsense weight values
-	    if (wt <= 0) continue;
-
-	    flux += value / wt;
-
-# else
-            // XXX skip the nonsense weight pixels?
-            if (unweighted_sum) {
-                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
-            } else {
-                wt = covarFactor * Wi->data.F32[yi][xi];
-                if (wt > 0) {
-                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
-                }
-            }
-# endif
-        }
-    }
-    return flux;
-}
Index: /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSourcePhotometry.h
===================================================================
--- /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSourcePhotometry.h	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psModules/src/objects/pmSourcePhotometry.h	(revision 34053)
@@ -38,5 +38,4 @@
 } pmSourcePhotometryMode;
 
-# if (HAVE_MODEL_VAR)
 typedef enum {
     PM_SOURCE_PHOTFIT_NONE       = 0,
@@ -44,6 +43,6 @@
     PM_SOURCE_PHOTFIT_IMAGE_VAR  = 2,
     PM_SOURCE_PHOTFIT_MODEL_VAR  = 3,
+    PM_SOURCE_PHOTFIT_MODEL_SKY  = 4,   // XXX bad name: set variance floor based on mean variance image (variance of sky)
 } pmSourceFitVarMode;
-# endif
 
 bool pmSourcePhotometryModel(
@@ -84,13 +83,7 @@
 bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);
 
-# if (HAVE_MODEL_VAR)
 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
 double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
 double pmSourceModelWeight(const pmSource *Mi, int term, const pmSourceFitVarMode fitVarMode, const float covarFactor, psImageMaskType maskVal);
-# else
-double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
-double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
-double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal);
-# endif
 
 bool pmSourceNeighborFlags (pmSource *source);
Index: /branches/eam_branches/ipp-20120601/psphot/src/psphot.h
===================================================================
--- /branches/eam_branches/ipp-20120601/psphot/src/psphot.h	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psphot/src/psphot.h	(revision 34053)
@@ -61,4 +61,7 @@
 bool            psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
 
+bool            psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool            psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+
 bool            psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
@@ -99,10 +102,5 @@
 
 bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, const char *filerule, bool final);
-
-# if (HAVE_MODEL_VAR)
 bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode);
-# else
-bool            psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final);
-# endif
 
 bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, const char *filerule, bool getPSFsize);
Index: /branches/eam_branches/ipp-20120601/psphot/src/psphotEfficiency.c
===================================================================
--- /branches/eam_branches/ipp-20120601/psphot/src/psphotEfficiency.c	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psphot/src/psphotEfficiency.c	(revision 34053)
@@ -420,9 +420,5 @@
 
     // psphotFitSourcesLinearReadout subtracts the model fits
-# if (HAVE_MODEL_VAR)
     if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true, PM_SOURCE_PHOTFIT_CONST)) {
-# else
-    if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) {
-# endif
         psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources.");
         psFree(fakeSources);
Index: /branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- /branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c	(revision 34053)
@@ -12,9 +12,7 @@
 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
 
-# if (HAVE_MODEL_VAR)
 bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources);
 pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe);
 bool psphotFreeModelVariance (pmReadout *readout, psArray *sources);
-# endif
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -30,5 +28,4 @@
     assert (recipe);
 
-# if (HAVE_MODEL_VAR)
     pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe);
     if (!fitVarMode) {
@@ -40,5 +37,4 @@
     // do a single pass.
     pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
-# endif
 
     int num = psphotFileruleCount(config, filerule);
@@ -68,15 +64,9 @@
         psAssert (psf, "missing psf?");
 
-# if (HAVE_MODEL_VAR)
-        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1))
-# else
-        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final))
-# endif
-	{
+        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
             return false;
         }
 
-# if (HAVE_MODEL_VAR)
 	// the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance
 	// after we have determined the initial set of fits, then we can generate the variance image and 
@@ -104,5 +94,4 @@
 	    }
 	}
-# endif
 
 	psphotVisualShowResidualImage (readout, (num > 0)); 
@@ -113,5 +102,4 @@
 }
 
-# if (HAVE_MODEL_VAR)
 // look up the fit variance mode from the recipe; older recipes do not have the value
 // 'LINEAR_FIT_VARIANCE_MODE'; in those cases, look for 'CONSTANT_PHOTOMETRIC_WEIGHTS' as a boolean and
@@ -136,4 +124,7 @@
 	return PM_SOURCE_PHOTFIT_IMAGE_VAR;
     }
+    if (!strcasecmp(fitVarModeString, "SKY")   || !strcasecmp(fitVarModeString, "MODEL_SKY")) {
+	return PM_SOURCE_PHOTFIT_MODEL_SKY;
+    }
     if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
 	return PM_SOURCE_PHOTFIT_MODEL_VAR;
@@ -142,12 +133,6 @@
     return PM_SOURCE_PHOTFIT_NONE;
 }
-# endif
-
-# if (HAVE_MODEL_VAR)
-bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode)
-# else
-bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final)
-# endif
-{
+
+bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) {
     bool status;
     float x;
@@ -185,11 +170,4 @@
     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
 
-# if (!HAVE_MODEL_VAR)
-    bool CONSTANT_PHOTOMETRIC_WEIGHTS =
-        psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
-    if (!status) {
-        psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
-    }
-# endif
     int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
     if (!status) {
@@ -319,5 +297,4 @@
     psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
 
-# if (HAVE_MODEL_VAR)
     // if fitVarMode is MODEL_VAR, then we need to generate the model image variance
     // XXX we have two possibilities here: 
@@ -327,5 +304,4 @@
 
     // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?)
-# endif
 
     // fill out the sparse matrix elements and border elements (B)
@@ -336,14 +312,9 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-# if (HAVE_MODEL_VAR)
         f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
-# else
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
         psSparseMatrixElement (sparse, i, i, f);
 
-# if (HAVE_MODEL_VAR)
-        // if we have used CONSTANT errors, then we need to calculate the value of the parameter error
-        if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) {
+        // if we have used CONSTANT errors, then we need to re-calculate the value of the parameter error
+        if (fitVarMode == PM_SOURCE_PHOTFIT_CONST) {
             float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
             errors->data.F32[i] = 1.0 / sqrt(var);
@@ -351,20 +322,7 @@
             errors->data.F32[i] = 1.0 / sqrt(f);
         }
-# else
-        // the formal error depends on the weighting scheme
-        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
-            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
-            errors->data.F32[i] = 1.0 / sqrt(var);
-        } else {
-            errors->data.F32[i] = 1.0 / sqrt(f);
-        }
-# endif
 
         // find the image x model value
-# if (HAVE_MODEL_VAR)
         f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
-# else
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
         psSparseVectorElement (sparse, i, f);
 
@@ -372,22 +330,11 @@
         switch (SKY_FIT_ORDER) {
           case 1:
-# if (HAVE_MODEL_VAR)
             f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 1, f);
             f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 2, f);
-# else
-            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-            psSparseBorderElementB (border, i, 1, f);
-            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-            psSparseBorderElementB (border, i, 2, f);
-# endif
 
           case 0:
-# if (HAVE_MODEL_VAR)
             f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal);
-# else
-            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
             psSparseBorderElementB (border, i, 0, f);
             break;
@@ -409,9 +356,5 @@
 
             // got an overlap; calculate cross-product and add to output array
-# if (HAVE_MODEL_VAR)
             f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal);
-# else
-            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
-# endif
             psSparseMatrixElement (sparse, j, i, f);
         }
@@ -451,9 +394,5 @@
 
     // set the sky, sky_x, sky_y components of border matrix
-# if (HAVE_MODEL_VAR)
     SetBorderMatrixElements (border, readout, fitSources, (fitVarMode == PM_SOURCE_PHOTFIT_CONST), SKY_FIT_ORDER, markVal);
-# else
-    SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
-# endif
 
     psSparseConstraint constraint;
@@ -614,5 +553,4 @@
 }
 
-# if (HAVE_MODEL_VAR)
 bool psphotModelBackgroundReadout(psImage *model,  // Model image
 				  psImage *modelStdev, // Model stdev image
@@ -650,4 +588,5 @@
     psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
 
+    // generate an image of the mean variance image in DN
     if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
 	psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
@@ -668,6 +607,13 @@
     psFree (varModelStdev);
 
-    // XXX for a test:
-    psphotSaveImage (NULL, modelVar, "model.bck.fits");
+    float gain = 1.0;  // accept 1.0 as a default since it is not critical to the analysis
+    pmCell *cell = readout->parent; // The parent cell
+    if (cell) {
+      gain = psMetadataLookupF32(&status, cell->concepts, "CELL.GAIN"); // Cell gain
+      if (!status) {
+	gain = 1.0;	      // set note above
+      }
+    }
+    if (gain > 2.0) { /* warn? */ }
 
     // insert all of the source models
@@ -690,10 +636,9 @@
 
 	// add the source model to the model variance image
+	// XXX note that this should be added with gain applied
+	// var_DN = flux_DN / gain [e/DN]
+	// to do this requires an API upgrade...
 	pmSourceAdd (source, PM_MODEL_OP_MODELVAR, maskVal);
     }
-
-    // XXX for a test:
-    psphotSaveImage (NULL, modelVar, "model.var.fits");
-    psphotSaveImage (NULL, readout->variance, "image.var.fits");
 
     // we save the model variance for future reference
@@ -724,3 +669,2 @@
     return true;
 }
-# endif
Index: /branches/eam_branches/ipp-20120601/psphot/src/psphotMaskReadout.c
===================================================================
--- /branches/eam_branches/ipp-20120601/psphot/src/psphotMaskReadout.c	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psphot/src/psphotMaskReadout.c	(revision 34053)
@@ -94,5 +94,5 @@
 
     // test output of files at this stage
-    if (psTraceGetLevel("psphot") >= 5) {
+    if (psTraceGetLevel("psphot.imsave") >= 5) {
         psphotSaveImage (NULL, readout->image,  "image.fits");
         psphotSaveImage (NULL, readout->mask,   "mask.fits");
@@ -105,2 +105,101 @@
     return true;
 }
+
+// XXX this function and support below was created to test the theory that the faint-end
+// bias results from the Poisson variation of the background pixels.  This is NOT the
+// case.  Using the code below maintains the faint-end bias.
+bool psphotUpdateVariance (pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    bool status = false;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+	// Generate the mask and weight images, including the user-defined analysis region of interest
+	if (!psphotUpdateVarianceReadout (config, view, filerule, i, recipe)) {
+	    psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for %s entry %d", filerule, i);
+	    return false;
+	}
+    }
+    return true;
+}
+
+// determine the mean variance image (equivalent to the background model, but for the variance image)
+// set the variance image to the MAX(input, mean)
+bool psphotUpdateVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+
+    bool status;
+
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
+    psAssert (file, "missing file?");
+
+    // find the currently selected readout
+    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmSourceFitVarMode varMode = psphotGetFitVarMode (recipe);
+    if (varMode == PM_SOURCE_PHOTFIT_NONE) {
+      psError (PSPHOT_ERR_CONFIG, false, "need valid LINEAR_FIT_VARIANCE_MODE");
+      return false;
+    }
+
+    // make this an option via the recipe
+    if (varMode != PM_SOURCE_PHOTFIT_MODEL_SKY) return true;
+
+    // create a model variance image (full-scale image to take result of psImageUnbin below)
+    psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
+
+    // find the binning information
+    psImageBinning *backBinning = psphotBackgroundBinning (modelVar, config);
+    assert (backBinning);
+    
+    psImage *varModel = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
+    psImage *varModelStdev = psImageAlloc(backBinning->nXruff, backBinning->nYruff, PS_TYPE_F32); // Background model
+
+    if (!psphotModelBackgroundReadout(varModel, varModelStdev, NULL, readout, backBinning, config, true)) {
+	psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
+	psFree (varModel);
+	psFree (varModelStdev);
+	return false;
+    }
+
+    // linear interpolation to full-scale
+    if (!psImageUnbin (modelVar, varModel, backBinning)) {
+	psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
+	psFree (varModel);
+	psFree (varModelStdev);
+	return false;
+    }
+
+    // XXX save these?
+    psFree (varModel);
+    psFree (varModelStdev);
+
+    psImage *im = readout->image;
+    psImage *wt = readout->variance;
+    for (int j = 0; j < im->numRows; j++) {
+      for (int i = 0; i < im->numCols; i++) {
+	if (!isfinite(im->data.F32[j][i])) continue;
+	if (!isfinite(wt->data.F32[j][i])) continue;
+	// XXX for a test, make variance constant wt->data.F32[j][i] = PS_MAX(wt->data.F32[j][i], modelVar->data.F32[j][i]);
+	wt->data.F32[j][i] = modelVar->data.F32[j][i];
+      }
+    }
+
+    // test output of files at this stage
+    if (psTraceGetLevel("psphot.imsave") >= 5) {
+        psphotSaveImage (NULL, readout->image,  "image.varsky.fits");
+        psphotSaveImage (NULL, readout->mask,   "mask.varsky.fits");
+        psphotSaveImage (NULL, readout->variance, "variance.varsky.fits");
+    }
+
+    psFree (modelVar);
+
+    return true;
+}
Index: /branches/eam_branches/ipp-20120601/psphot/test/tap_psphot_varmodel.pro
===================================================================
--- /branches/eam_branches/ipp-20120601/psphot/test/tap_psphot_varmodel.pro	(revision 34052)
+++ /branches/eam_branches/ipp-20120601/psphot/test/tap_psphot_varmodel.pro	(revision 34053)
@@ -2,85 +2,61 @@
 # -*-sh-*-
 
-# $KAPA = kapa -noX
-
-# PSF.CONVOLVE : if true, we insert delta functions (and optionally
-#                galaxies) and smooth the image with the psf model
-#                (uses a GAUSS regardless of the model). Note that
-#                PSF.CONVOLVE = T is faster than F, but (a) only
-#                allows Gauss models and (b) only yields quantized
-#                locations
-
-# config for ppImage to generate chip, mask, weight
-$ppImageConfig = -recipe PPIMAGE PPIMAGE_N
-$ppImageConfig = $ppImageConfig -Db BACKGROUND T
-$ppImageConfig = $ppImageConfig -Db CHIP.FITS T
-$ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T
-$ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T
-$ppImageConfig = $ppImageConfig -Db BASE.FITS F
-$ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T
-$ppImageConfig = $ppImageConfig -Db PHOTOM F
-
-# basic options for the these images (filter, location, obstype)
-$BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0
-$BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
-
-# create an image with fake sources and insert the resulting cmf file into a dvodb
-$RefConfig = -camera SIMTEST 
-$RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE 
-$RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref 
-$RefConfig = $RefConfig -Db PSF.CONVOLVE F
-
-# options for the reference image
-$RefOptions = $BaseOptions -exptime 100.0 
-$RefOptions = $RefOptions -seeing 1.0 
-$RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS 
-$RefOptions = $RefOptions -Df STARS.DENSITY 10.0 
-$RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
-
-# basic config for ppSim with randomly distributed stars and NO galaxies
-$RealConfig = -camera SIMTEST 
-$RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN 
-$RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
-$RealConfig = $RealConfig -Db STARS.FAKE F
-$RealConfig = $RealConfig -Db STARS.REAL T 
-$RealConfig = $RealConfig -Db MATCH.DENSITY F 
-$RealConfig = $RealConfig -Db PSF.CONVOLVE F
-$RealConfig = $RealConfig -Df STARS.DENSITY 10.0
-$RealConfig = $RealConfig -Df STARS.SIGMA.LIM 2.5
-$RealConfig = $RealConfig -Db GALAXY.FAKE F 
-$RealConfig = $RealConfig -Db GALAXY.GRID F 
-
-# options for the repeated images
-$RealOptions = $BaseOptions -exptime 30.0
-  
-$ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS
-
-# sample alternate options:
-# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
-# $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2
-# $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0
-# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
-
-list fwhm 
- 1.0 
- 1.1 
- 1.2 
- 1.5
-end
+# This script includes a set of tests to demonstrate the dependence of the faint-end bias on the weighting scheme
+# We have 3 weighting schemes :
+# CONSTANT -- per-pixel weight is fixed (disadvantage: lower S/N, especially for higher sky?)
+# IMAGE_VAR - per-pixel weight is Poisson from image (disadvantage: faint-end bias)
+# MODEL_VAR - per-pixel weight is Poisson from model (disadvantage: 2 linear-fit passes)
+
+# Functions:
+# 
+# init  : initialize variables
+# mkref : generate a fake reference catalog in a DVO database
+# mkexp : generate a fake exposure from fake catalog and detrend it (saves basename.in.cmf & basename.fits)
+# runphot : run psphot on an exposure (saves basename.cmf)
+# ckchip.mags : generate a set of summary plots for the psphot analysis
+
+# I would like to produce a grid of tests:
+# sky (bright, middle, dark)
+# fwhm (1.0, 1.3, 1.6)
+# PSF_MODEL input (GAUSS PS1_V1)
+# PSF_MODEL apply (GAUSS PS1_V1)
+# variance mode: CONSTANT, IMAGE_VAR, MODEL_VAR
 
 macro go
-
   mkdir test
 
-  $ExtraOptions = -D PSF.MODEL PS_MODEL_GAUSS
-  mkexp test/image.00 1.0
-  $ExtraOptions = -D PSF.MODEL PS_MODEL_PGAUSS
-  mkexp test/image.01 1.0
-  $ExtraOptions = -D PSF.MODEL PS_MODEL_PS1_V1
-  mkexp test/image.02 1.0
+  local sky fwhm psf_in psf_out
+
+  # mkref creates refimage.* and catdir.ref
+  file catdir.ref found
+  if (not($found))
+    mkref
+  end
+
+  $KAPA = kapa -noX
+  if (not($?RefConfig)) init
+
+  foreach sky 19.0 20.0 21.0
+    foreach fwhm 1.0 1.3 1.6
+      foreach psf_in GAUSS PS1_V1
+
+        sprintf name "test/test.%02d.%02d.%s" $sky {10*$fwhm} $psf_in
+        # mkexp $name $sky $fwhm $psf_in
+        
+        foreach psf_out GAUSS PS1_V1
+          foreach mode CONSTANT IMAGE_VAR MODEL_VAR
+            # runphot $name $name.$psf_out.$mode "-D LINEAR_FIT_VARIANCE_MODE $mode -D PSF_MODEL PS_MODEL_$psf_out"
+            ckchip.mags $name.in.cmf $name.$psf_out.$mode.cmf $name.$psf_out.$mode 0.0
+          end
+        end
+      end
+    end
+  end
 end
 
 # create a reference database of fake stars to be used by ppSim below
 macro mkref
+  if (not($?RefConfig)) init
+
   exec rm -rf catdir.ref
   exec rm -f refimage.fits
@@ -100,4 +76,29 @@
 # create a realistic distribution of fake stars, GAUSS PSF
 macro mkexp
+  if ($0 != 5)
+    echo "USAGE: mkexp basename sky fwhm psf_model"
+    break
+  end
+
+  local fwhm basename psf_model
+  $basename = $1
+  $sky = $2
+  $fwhm = $3
+  $psf_model = $4
+
+  $ExtraOptions = -D PSF.MODEL PS_MODEL_$psf_model
+
+  # create the raw image
+  echo ppSim -seeing $fwhm -skymags $sky -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
+  exec ppSim -seeing $fwhm -skymags $sky -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
+  exec /bin/mv -f $basename.cmf $basename.in.cmf
+
+  # create the chip output
+  echo ppImage $ppImageConfig -file $basename.fits $basename
+  exec ppImage $ppImageConfig -file $basename.fits $basename
+end
+
+# create a realistic distribution of fake stars, GAUSS PSF
+macro mkexp.deep
   if ($0 != 3)
     echo "USAGE: mkexp basename fwhm"
@@ -109,7 +110,35 @@
   $fwhm = $2
 
+  $RealOptionsDeep = $BaseOptions -exptime 100
+
   # create the raw image
-  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
-  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptions $ExtraOptions $basename
+  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
+  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
+  exec /bin/mv -f $basename.cmf $basename.in.cmf
+
+  # create the chip output
+  echo ppImage $ppImageConfig -file $basename.fits $basename
+  exec ppImage $ppImageConfig -file $basename.fits $basename
+end
+
+# create a realistic distribution of fake stars, GAUSS PSF
+macro mkexp.bright
+  if ($0 != 3)
+    echo "USAGE: mkexp basename fwhm"
+    break
+  end
+
+  local fwhm basename
+  $basename = $1
+  $fwhm = $2
+
+  # basic options for the these images (filter, location, obstype)
+  $BaseOptions = -type OBJECT -filter r -skymags 20.0 -ra 270.70 -dec -23.70 -pa 0.0
+  $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
+  $RealOptionsDeep = $BaseOptions -exptime 100
+
+  # create the raw image
+  echo ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
+  exec ppSim -seeing $fwhm -nx 3000 -ny 3000 $RealConfig $RealOptionsDeep $ExtraOptions $basename
   exec /bin/mv -f $basename.cmf $basename.in.cmf
 
@@ -133,4 +162,98 @@
   echo psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options
   exec psphot -threads 4 -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $outname $options
+end
+
+# compare two cmf files with extname Chip.psf 
+# things to compare:
+# * completeness (which sources in (1) are not detected in (2)
+# * positions (X_PSF, Y_PSF) 
+# * instrumental psf mags
+# * position errors (no input errors; use a model?)
+# * measured FWHM?
+# * kron mags (fluxes)
+# * etc, etc
+macro ckchip.mags
+  if ($0 != 5)
+    echo "USAGE: ckchip.mags (raw) (out) (output) (zpt_off)"
+    break
+  end
+
+list pairs
+  PSF_INST_MAG_out      M_raw                 PSF_INST_MAG_raw 2 -0.21 0.21 V
+  AP_MAG_out            M_raw                 PSF_INST_MAG_raw 2 -0.21 0.21 V
+end
+
+  load.cmf $1 Chip.psf raw
+  load.cmf $2 Chip.psf out
+
+  # images generated with convolution will not have the right output positions
+  set X_raw = int(X_PSF_raw) + 0.5
+  set Y_raw = int(Y_PSF_raw) + 0.5
+  set M_raw = PSF_INST_MAG_raw + $4
+  set K_out = -2.5*log(KRON_FLUX_out)
+  match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2
+
+  local i NX NY nx ny N
+
+  resize 1000 1000
+
+  clear
+  section a0 0.0 0.0 1.0 0.5
+  label -fn courier 14; 
+  style -pt 0 -sz 0.4
+  show.pair 0
+  delete -q imag_V dmag_V smag_V Dmag_V dmag_Vraw
+  for imag -11 -6 0.2
+    subset dmsub = delta if (rv > $imag) && (rv < $imag + 0.2)
+    vstat -q dmsub
+    concat $imag imag_V
+    concat $MEAN dmag_V
+    concat $MEDIAN Dmag_V
+    concat $SIGMA smag_V
+  end
+
+  vstat -q dmag_V
+  $offset = $MEDIAN
+  set dmag_V = dmag_V - $offset
+  set Dmag_V = Dmag_V - $offset
+
+  section a1 0.0 0.50 1.0 0.25
+  lim rv -0.025 0.075; label -fn courier 14; box; 
+  plot -c red   -pt 7 -sz 2.0 imag_V dmag_V
+  plot -c darkgreen -pt 3 -sz 2.0 imag_V Dmag_V
+  # plot -c gold -pt 4 -sz 2.0 imag_V dmag_Vraw
+  plot -c blue  -pt 2 -sz 2.0 imag_V smag_V
+  # label -y "mean (red), median (green), unclipped mean (gold), sigma (blue) of delta"
+  label -y "mean (red), median (green), sigma (blue) of delta"
+  sprintf line "OFFSET: %6.3f" $offset
+  textline -frac 0.1 0.8 -fn courier 24 "$line"
+
+  set lChiNorm = log(PSF_CHISQ_out / PSF_NDOF_out)
+  reindex chi = lChiNorm using index1
+  section a2 0.0 0.75 1.0 0.25
+  label -fn courier 14; 
+  lim rv chi; box; 
+  plot -c red -pt 0 -sz 0.5 rv chi
+
+  label -y "chisq"
+  png -name $3.png
+
+  # section a1 0.0 0.5 1.0 0.5
+  # style -pt 0 -sz 0.4
+  # show.pair 1
+end
+
+macro go.phot
+ runphot test.00 test.00.varmode.C "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE CONSTANT"
+ runphot test.00 test.00.varmode.I "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE IMAGE_VAR"
+ runphot test.00 test.00.varmode.M "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE MODEL_VAR"
+ runphot test.00 test.00.varmode.S "-Db SAVE.RESID T -D LINEAR_FIT_VARIANCE_MODE MODEL_SKY"
+end
+
+macro go.vars
+ dev -n varC; ckchip.mags test.00.in.cmf test.00.varmode.C.cmf test.00.varmode.C 0.0
+ dev -n varI; ckchip.mags test.00.in.cmf test.00.varmode.I.cmf test.00.varmode.I 0.0
+ dev -n varM; ckchip.mags test.00.in.cmf test.00.varmode.M.cmf test.00.varmode.M 0.0
+ dev -n varS; ckchip.mags test.00.in.cmf test.00.varmode.S.cmf test.00.varmode.S 0.0
 end
 
@@ -347,5 +470,6 @@
     lim rv $word:4 $word:5; box; plot rv delta
   end
-  label -y '$word:0' -x '$word:2'
+  $line = '$word:0' - '$word:1'
+  label -y "$line" -x '$word:2'
 end
 
@@ -637,4 +761,71 @@
 end
 
+macro init
+  # config for ppImage to generate chip, mask, weight
+  $ppImageConfig = -recipe PPIMAGE PPIMAGE_N
+  $ppImageConfig = $ppImageConfig -Db BACKGROUND T
+  $ppImageConfig = $ppImageConfig -Db CHIP.FITS T
+  $ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T
+  $ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T
+  $ppImageConfig = $ppImageConfig -Db BASE.FITS F
+  $ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T
+  $ppImageConfig = $ppImageConfig -Db PHOTOM F
+  
+  # basic options for the these images (filter, location, obstype)
+  $BaseOptions = -type OBJECT -filter r -ra 270.70 -dec -23.70 -pa 0.0
+  $BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0
+  
+  # PSF.CONVOLVE : if true, we insert delta functions (and optionally
+  #                galaxies) and smooth the image with the psf model
+  #                (uses a GAUSS regardless of the model). Note that
+  #                PSF.CONVOLVE = T is faster than F, but (a) only
+  #                allows Gauss models and (b) only yields quantized
+  #                locations
+
+  # create an image with fake sources and insert the resulting cmf file into a dvodb
+  $RefConfig = -camera SIMTEST 
+  $RefConfig = $RefConfig -recipe PPSIM STACKTEST.MAKE 
+  $RefConfig = $RefConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref 
+  $RefConfig = $RefConfig -Db PSF.CONVOLVE F
+  
+  # options for the reference image
+  $RefOptions = $BaseOptions
+  $RefOptions = $RefOptions -exptime 100.0 
+  $RefOptions = $RefOptions -seeing 1.0 
+  $RefOptions = $RefOptions -skymags 21.0  
+  $RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS 
+  $RefOptions = $RefOptions -Df STARS.DENSITY 10.0 
+  $RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
+
+  # basic config for ppSim with randomly distributed stars and NO galaxies
+  $RealConfig = -camera SIMTEST 
+  $RealConfig = $RealConfig -recipe PPSIM STACKTEST.RUN 
+  $RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
+  $RealConfig = $RealConfig -Db STARS.FAKE F
+  $RealConfig = $RealConfig -Db STARS.REAL T 
+  $RealConfig = $RealConfig -Db MATCH.DENSITY F 
+  $RealConfig = $RealConfig -Db PSF.CONVOLVE F
+  $RealConfig = $RealConfig -Df STARS.DENSITY 10.0
+  $RealConfig = $RealConfig -Df STARS.SIGMA.LIM 2.5
+  $RealConfig = $RealConfig -Db GALAXY.FAKE F 
+  $RealConfig = $RealConfig -Db GALAXY.GRID F 
+  
+  # options for the repeated images
+  $RealOptions = $BaseOptions -exptime 30.0
+    
+  # sample alternate options:
+  # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
+  # $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2
+  # $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0
+  # $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
+  
+  list fwhm 
+   1.0 
+   1.1 
+   1.2 
+   1.5
+  end
+end
+
 if ($SCRIPT)
   fulltest 4