Index: /trunk/psModules/src/objects/Makefile.am
===================================================================
--- /trunk/psModules/src/objects/Makefile.am	(revision 5843)
+++ /trunk/psModules/src/objects/Makefile.am	(revision 5844)
@@ -15,5 +15,5 @@
 	models/pmModel_QGAUSS.c \
 	models/pmModel_SGAUSS.c 
-    
+
 psmoduleincludedir = $(includedir)
 psmoduleinclude_HEADERS = \
Index: /trunk/psModules/src/objects/pmModelGroup.c
===================================================================
--- /trunk/psModules/src/objects/pmModelGroup.c	(revision 5843)
+++ /trunk/psModules/src/objects/pmModelGroup.c	(revision 5844)
@@ -10,14 +10,59 @@
 #include "models/pmModel_SGAUSS.c"
 
-static pmModelGroup models[] = {
-                                   {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
-                                   {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
-                                   {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
-                                   {"PS_MODEL_SGAUSS",       9, pmModelFunc_SGAUSS,  pmModelFlux_SGAUSS,  pmModelRadius_SGAUSS,  pmModelLimits_SGAUSS,  pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS},
-                               };
+static pmModelGroup defaultModels[] = {
+                                          {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
+                                          {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
+                                          {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
+                                          {"PS_MODEL_SGAUSS",       9, pmModelFunc_SGAUSS,  pmModelFlux_SGAUSS,  pmModelRadius_SGAUSS,  pmModelLimits_SGAUSS,  pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS},
+                                      };
+
+static pmModelGroup *models = NULL;
+static int Nmodels = 0;
+
+static void ModelGroupFree (pmModelGroup *modelGroup)
+{
+
+    if (modelGroup == NULL)
+        return;
+    psFree (modelGroup);
+    return;
+}
+
+pmModelGroup *pmModelGroupAlloc (int nModels)
+{
+
+    pmModelGroup *modelGroup = (pmModelGroup *) psAlloc (nModels * sizeof(pmModelGroup));
+    psMemSetDeallocator(modelGroup, (psFreeFunc) ModelGroupFree);
+    return (modelGroup);
+}
+
+void pmModelGroupAdd (pmModelGroup *model)
+{
+
+    if (models == NULL) {
+        pmModelGroupInit ();
+    }
+
+    Nmodels ++;
+    models = (pmModelGroup *) psRealloc (models, Nmodels*sizeof(pmModelGroup));
+    models[Nmodels-1] = model[0];
+    return;
+}
+
+void pmModelGroupInit (void)
+{
+
+    int Nnew = sizeof (defaultModels) / sizeof (pmModelGroup);
+
+    models = pmModelGroupAlloc (Nnew);
+    for (int i = 0; i < Nnew; i++) {
+        models[i] = defaultModels[i];
+    }
+    Nmodels = Nnew;
+    return;
+}
 
 pmModelFunc pmModelFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -29,5 +74,4 @@
 pmModelFlux pmModelFlux_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -39,5 +83,4 @@
 pmModelRadius pmModelRadius_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -49,5 +92,4 @@
 pmModelLimits pmModelLimits_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -59,5 +101,4 @@
 pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -69,5 +110,4 @@
 pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -79,5 +119,4 @@
 pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -89,5 +128,4 @@
 psS32 pmModelParameterCount (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -99,6 +137,4 @@
 psS32 pmModelSetType (char *name)
 {
-
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     for (int i = 0; i < Nmodels; i++) {
         if (!strcmp(models[i].name, name)) {
@@ -111,6 +147,4 @@
 char *pmModelGetType (pmModelType type)
 {
-
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
Index: /trunk/psModules/src/objects/pmModelGroup.h
===================================================================
--- /trunk/psModules/src/objects/pmModelGroup.h	(revision 5843)
+++ /trunk/psModules/src/objects/pmModelGroup.h	(revision 5844)
@@ -9,6 +9,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-10-10 19:53:40 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -209,3 +209,12 @@
 pmModelGroup;
 
+// allocate a pmModelGroup to hold nModels entries
+pmModelGroup *pmModelGroupAlloc (int nModels);
+
+// initialize the internal (static) model group with the default models
+void pmModelGroupInit (void);
+
+// add a new model to the internal (static) model group
+void pmModelGroupAdd (pmModelGroup *model);
+
 # endif
Index: /trunk/psModules/src/objects/pmObjects.c
===================================================================
--- /trunk/psModules/src/objects/pmObjects.c	(revision 5843)
+++ /trunk/psModules/src/objects/pmObjects.c	(revision 5844)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -85,5 +85,5 @@
 psVector containing the specified row of data from the psImage.
  
-XXX: Is there a better way to do this?  
+XXX: Is there a better way to do this?
 XXX EAM: does this really need to alloc a new vector???
 *****************************************************************************/
@@ -271,4 +271,7 @@
     tmp->chisq = 0.0;
     tmp->nIter = 0;
+    tmp->radius = 0;
+    tmp->status = PM_MODEL_UNTRIED;
+
     psS32 Nparams = pmModelParameterCount(type);
     if (Nparams == 0) {
@@ -850,6 +853,7 @@
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
         for (psS32 col = 0; col < source->pixels->numCols ; col++) {
-            if ((source->mask != NULL) && (source->mask->data.U8[row][col]))
+            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {
                 continue;
+            }
 
             psF32 xDiff = col + source->pixels->col0 - xPeak;
@@ -858,6 +862,7 @@
             // XXX EAM : calculate xDiff, yDiff up front;
             //           radius is just a function of (xDiff, yDiff)
-            if (!VALID_RADIUS(xDiff, yDiff, R2))
+            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
                 continue;
+            }
 
             psF32 pDiff = source->pixels->data.F32[row][col] - sky;
@@ -865,6 +870,7 @@
             // XXX EAM : check for valid S/N in pixel
             // XXX EAM : should this limit be user-defined?
-            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1)
+            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {
                 continue;
+            }
 
             Sum += pDiff;
@@ -970,6 +976,6 @@
 
 /******************************************************************************
-pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 
-sigma-x, sigma-y plane. return 0,0 clump in case of error. 
+pmSourcePSFClump(source, metadata): Find the likely PSF clump in the
+sigma-x, sigma-y plane. return 0,0 clump in case of error.
 *****************************************************************************/
 
@@ -994,8 +1000,16 @@
         psImage *splane = NULL;
         int binX, binY;
+        bool status;
+
+        psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
+        if (!status)
+            SX_MAX = 10.0;
+        psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
+        if (!status)
+            SY_MAX = 10.0;
 
         // construct a sigma-plane image
         // psImageAlloc does zero the data
-        splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
+        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
         for (int i = 0; i < splane->numRows; i++)
         {
@@ -1131,5 +1145,5 @@
 XXX: How can this function ever return FALSE?
  
-XXX EAM : add the saturated mask value to metadata 
+XXX EAM : add the saturated mask value to metadata
 *****************************************************************************/
 
@@ -1146,4 +1160,5 @@
     int Ncr      = 0;
     int Nsatstar = 0;
+    psRegion allArray = psRegionSet (0, 0, 0, 0);
 
     psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
@@ -1177,10 +1192,5 @@
 
         // XXX EAM : can we use the value of SATURATE if mask is NULL?
-        //
-        // XXX: Must verify this region (the region argument was added to psImageCountPixelMask()
-        // after EAM wrote this code.
-        //
-        psRegion tmpRegion = psRegionSet(0, tmpSrc->mask->numCols, 0, tmpSrc->mask->numRows);
-        int Nsatpix = psImageCountPixelMask(tmpSrc->mask, tmpRegion, PSPHOT_MASK_SATURATED);
+        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED);
 
         // saturated star (size consistent with PSF or larger)
@@ -1254,5 +1264,5 @@
 
 /** pmSourceDefinePixels()
- * 
+ *
  * Define psImage subarrays for the source located at coordinates x,y on the
  * image set defined by readout. The pixels defined by this operation consist of
@@ -1266,7 +1276,7 @@
  * example). This function should be used to define a region of interest around a
  * source, including both source and sky pixels.
- * 
+ *
  * XXX: must code this.
- * 
+ *
  */
 bool pmSourceDefinePixels(
@@ -1360,7 +1370,7 @@
 /******************************************************************************
 pmSourceModelGuess(source, model): This function allocates a new
-pmModel structure based on the given modelType specified in the argument list.  
-The corresponding pmModelGuess function is returned, and used to 
-supply the values of the params array in the pmModel structure.  
+pmModel structure based on the given modelType specified in the argument list.
+The corresponding pmModelGuess function is returned, and used to
+supply the values of the params array in the pmModel structure.
  
 XXX: Many parameters are based on the src->moments structure, which is in
@@ -1598,4 +1608,5 @@
         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        model->status = PM_MODEL_BADARGS;
         return(false);
     }
@@ -1649,19 +1660,10 @@
     }
 
-    // XXX EAM: we need to do something (give an error?) if rc is false
-    // XXX EAM: psMinimizeLMChi2 does not check convergence
-
-    // XXX models can go insane: reject these
-    onPic &= (params->data.F32[2] >= source->pixels->col0);
-    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
-    onPic &= (params->data.F32[3] >= source->pixels->row0);
-    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
-
-    // XXX EAM: save the resulting chisq, nDOF, nIter
+    // save the resulting chisq, nDOF, nIter
     model->chisq = myMin->value;
     model->nIter = myMin->iter;
     model->nDOF  = y->n - nParams;
 
-    // XXX EAM get the Gauss-Newton distance for fixed model parameters
+    // get the Gauss-Newton distance for fixed model parameters
     if (paramMask != NULL) {
         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
@@ -1673,4 +1675,20 @@
         }
         psFree (delta);
+    }
+
+    // set the model success or failure status
+    if (!fitStatus) {
+        model->status = PM_MODEL_NONCONVERGE;
+    } else {
+        model->status = PM_MODEL_SUCCESS;
+    }
+
+    // models can go insane: reject these
+    onPic &= (params->data.F32[2] >= source->pixels->col0);
+    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
+    onPic &= (params->data.F32[3] >= source->pixels->row0);
+    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
+    if (!onPic) {
+        model->status = PM_MODEL_OFFIMAGE;
     }
 
@@ -1732,4 +1750,5 @@
         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        model->status = PM_MODEL_BADARGS;
         return(false);
     }
@@ -1810,10 +1829,4 @@
     // XXX EAM: psMinimizeLMChi2 does not check convergence
 
-    // XXX models can go insane: reject these
-    onPic &= (params->data.F32[2] >= source->pixels->col0);
-    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
-    onPic &= (params->data.F32[3] >= source->pixels->row0);
-    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
-
     // XXX EAM: save the resulting chisq, nDOF, nIter
     model->chisq = myMin->value;
@@ -1832,8 +1845,26 @@
     }
 
-    psFree(paramMask);
+    // set the model success or failure status
+    if (!fitStatus) {
+        model->status = PM_MODEL_NONCONVERGE;
+    } else {
+        model->status = PM_MODEL_SUCCESS;
+    }
+
+    // XXX models can go insane: reject these
+    onPic &= (params->data.F32[2] >= source->pixels->col0);
+    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
+    onPic &= (params->data.F32[3] >= source->pixels->row0);
+    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
+    if (!onPic) {
+        model->status = PM_MODEL_OFFIMAGE;
+    }
+
     psFree(x);
     psFree(y);
+    psFree(yErr);
     psFree(myMin);
+    psFree(covar);
+    psFree(paramMask);
 
     rc = (onPic && fitStatus);
@@ -1846,5 +1877,7 @@
                              pmModel *model,
                              bool center,
-                             psS32 flag)
+                             bool sky,
+                             bool add
+                                )
 {
     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
@@ -1859,4 +1892,6 @@
     psS32 imageCol;
     psS32 imageRow;
+    psF32 skyValue = params->data.F32[0];
+    psF32 pixelValue;
 
     for (psS32 i = 0; i < image->numRows; i++) {
@@ -1864,5 +1899,5 @@
             if ((mask != NULL) && mask->data.U8[i][j])
                 continue;
-            psF32 pixelValue;
+
             // XXX: Should you be adding the pixels for the entire subImage,
             // or a radius of pixels around it?
@@ -1882,15 +1917,21 @@
             x->data.F32[0] = (float) imageCol;
             x->data.F32[1] = (float) imageRow;
-            pixelValue = modelFunc (NULL, params, x);
-            // fprintf (stderr, "%f %f  %d %d  %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue);
-
-            if (flag == 1) {
-                pixelValue = -pixelValue;
-            }
-
-            // XXX: Must figure out how to calculate the image coordinates and
-            // how to use the boolean "center" flag.
-
-            image->data.F32[i][j]+= pixelValue;
+
+            // set the appropriate pixel value for this coordinate
+            if (sky) {
+                pixelValue = modelFunc (NULL, params, x);
+            } else {
+                pixelValue = modelFunc (NULL, params, x) - skyValue;
+            }
+
+
+            // add or subtract the value
+            if (add
+               ) {
+                image->data.F32[i][j] += pixelValue;
+            }
+            else {
+                image->data.F32[i][j] -= pixelValue;
+            }
         }
     }
@@ -1907,8 +1948,9 @@
                       psImage *mask,
                       pmModel *model,
-                      bool center)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0);
+                      bool center,
+                      bool sky)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);
     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
@@ -1920,8 +1962,9 @@
                       psImage *mask,
                       pmModel *model,
-                      bool center)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1);
+                      bool center,
+                      bool sky)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);
     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
Index: /trunk/psModules/src/objects/pmObjects.h
===================================================================
--- /trunk/psModules/src/objects/pmObjects.h	(revision 5843)
+++ /trunk/psModules/src/objects/pmObjects.h	(revision 5844)
@@ -4,5 +4,5 @@
  * images is one of the critical tasks of the IPP or any astronomical software
  * system. This file will define structures and functions related to the task
- * of source detection and measurement. The elements defined in this section 
+ * of source detection and measurement. The elements defined in this section
  * are generally low-level components which can be connected together to
  * construct a complete object measurement suite.
@@ -10,6 +10,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -45,5 +45,5 @@
 
 /** pmPeakType
- * 
+ *
  *  A peak pixel may have several features which may be determined when the
  *  peak is found or measured. These are specified by the pmPeakType enum.
@@ -52,5 +52,5 @@
  *  edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific
  *  number of neighbors at the same value, within some tolarence:
- * 
+ *
  */
 typedef enum {
@@ -63,9 +63,9 @@
 
 /** pmPeak data structure
- *  
+ *
  *  A source has the capacity for several types of measurements. The
  *  simplest measurement of a source is the location and flux of the peak pixel
  *  associated with the source:
- *  
+ *
  */
 typedef struct
@@ -80,9 +80,9 @@
 
 /** pmMoments data structure
- *  
+ *
  * One of the simplest measurements which can be made quickly for an object
  * are the object moments. We specify a structure to carry the moment information
  * for a specific source:
- *  
+ *
  */
 typedef struct
@@ -103,9 +103,9 @@
 
 /** pmPSFClump data structure
- * 
+ *
  * A collection of object moment measurements can be used to determine
  * approximate object classes. The key to this analysis is the location and
  * statistics (in the second-moment plane,
- *  
+ *
  */
 typedef struct
@@ -118,13 +118,17 @@
 pmPSFClump;
 
+// type of model carried by the pmModel structure
 typedef int pmModelType;
-#define PS_MODEL_GAUSS 0
-#define PS_MODEL_PGAUSS 1
-#define PS_MODEL_QGAUSS 2
-#define PS_MODEL_SGAUSS 3
-
+
+typedef enum {
+    PM_MODEL_UNTRIED,               ///< model fit not yet attempted
+    PM_MODEL_SUCCESS,               ///< model fit succeeded
+    PM_MODEL_NONCONVERGE,           ///< model fit did not converge
+    PM_MODEL_OFFIMAGE,              ///< model fit drove out of range
+    PM_MODEL_BADARGS                ///< model fit called with invalid args
+} pmModelStatus;
 
 /** pmModel data structure
- * 
+ *
  * Every source may have two types of models: a PSF model and a FLT (floating)
  * model. The PSF model represents the best fit of the image PSF to the specific
@@ -132,5 +136,5 @@
  * object by the PSF, not by the fit. The FLT model represents the best fit of
  * the given model to the object, with all parameters floating in the fit.
- *  
+ *
  */
 typedef struct
@@ -142,4 +146,5 @@
     int nDOF;    ///< number of degrees of freedom
     int nIter;    ///< number of iterations to reach min
+    int status;         ///< fit status
     float radius;   ///< fit radius actually used
 }
@@ -147,12 +152,12 @@
 
 /** pmSourceType enumeration
- * 
+ *
  * A given source may be identified as most-likely to be one of several source
  * types. The pmSource entry pmSourceType defines the current best-guess for this
  * source.
- * 
+ *
  * XXX: The values given below are currently illustrative and will require
  * some modification as the source classification code is developed. (TBD)
- * 
+ *
  */
 typedef enum {
@@ -178,9 +183,9 @@
 
 /** pmSource data structure
- *  
+ *
  *  This source has the capacity for several types of measurements. The
  *  simplest measurement of a source is the location and flux of the peak pixel
  *  associated with the source:
- *  
+ *
  */
 typedef struct
@@ -194,4 +199,6 @@
     pmModel *modelFLT;   ///< FLT (floating) Model fit (parameters and type).
     pmSourceType type;   ///< Best identification of object.
+    float apMag;
+    float fitMag;
 }
 pmSource;
@@ -211,5 +218,5 @@
 
 /** pmMomentsAlloc()
- * 
+ *
  */
 pmMoments *pmMomentsAlloc();
@@ -217,5 +224,5 @@
 
 /** pmModelAlloc()
- * 
+ *
  */
 pmModel *pmModelAlloc(pmModelType type);
@@ -223,5 +230,5 @@
 
 /** pmSourceAlloc()
- * 
+ *
  */
 pmSource  *pmSourceAlloc();
@@ -229,5 +236,5 @@
 
 /** pmFindVectorPeaks()
- * 
+ *
  * Find all local peaks in the given vector above the given threshold. A peak
  * is defined as any element with a value greater than its two neighbors and with
@@ -243,5 +250,5 @@
  * a vector containing the coordinates (element number) of the detected peaks
  * (type psU32).
- * 
+ *
  */
 psVector *pmFindVectorPeaks(
@@ -252,5 +259,5 @@
 
 /** pmFindImagePeaks()
- * 
+ *
  * Find all local peaks in the given image above the given threshold. This
  * function should find all row peaks using pmFindVectorPeaks, then test each row
@@ -263,5 +270,5 @@
  * peaks at the (+x,+y) corner. When selecting the peaks, their type must also be
  * set. The result of this function is an array of pmPeak entries.
- * 
+ *
  */
 psArray *pmFindImagePeaks(
@@ -272,8 +279,8 @@
 
 /** pmCullPeaks()
- * 
+ *
  * Eliminate peaks from the psList that have a peak value above the given
  * maximum, or fall outside the valid region.
- * 
+ *
  */
 psList *pmCullPeaks(
@@ -285,5 +292,5 @@
 
 /** pmPeaksSubset()
- * 
+ *
  * Create a new peaks array, removing certain types of peaks from the input
  * array of peaks based on the given criteria. Peaks should be eliminated if they
@@ -291,5 +298,5 @@
  * the valid region.  The result of the function is a new array with a reduced
  * number of peaks.
- * 
+ *
  */
 psArray *pmPeaksSubset(
@@ -301,5 +308,5 @@
 
 /** pmSourceDefinePixels()
- * 
+ *
  * Define psImage subarrays for the source located at coordinates x,y on the
  * image set defined by readout. The pixels defined by this operation consist of
@@ -313,7 +320,7 @@
  * example). This function should be used to define a region of interest around a
  * source, including both source and sky pixels.
- * 
+ *
  * XXX: must code this.
- * 
+ *
  */
 // XXX: Uncommenting the pmReadout causes compile errors.
@@ -328,5 +335,5 @@
 
 /** pmSourceLocalSky()
- * 
+ *
  * Measure the local sky in the vicinity of the given source. The Radius
  * defines the square aperture in which the moments will be measured. This
@@ -338,5 +345,5 @@
  * This function allocates the pmMoments structure. The resulting sky is used to
  * set the value of the pmMoments.sky element of the provided pmSource structure.
- * 
+ *
  */
 bool pmSourceLocalSky(
@@ -348,5 +355,5 @@
 
 /** pmSourceMoments()
- * 
+ *
  * Measure source moments for the given source, using the value of
  * source.moments.sky provided as the local background value and the peak
@@ -355,5 +362,5 @@
  * are measured within the given circular radius of the source.peak coordinates.
  * The return value indicates the success (TRUE) of the operation.
- * 
+ *
  */
 bool pmSourceMoments(
@@ -364,5 +371,5 @@
 
 /** pmSourcePSFClump()
- * 
+ *
  * We use the source moments to make an initial, approximate source
  * classification, and as part of the information needed to build a PSF model for
@@ -373,5 +380,5 @@
  * similar. The function returns a pmPSFClump structure, representing the
  * centroid and size of the clump in the sigma_x, sigma_y second-moment plane.
- * 
+ *
  * The goal is to identify and characterize the stellar clump within the
  * sigma_x, sigma_y second-moment plane.  To do this, an image is constructed to
@@ -384,10 +391,10 @@
  *  * used to calculate the median and standard deviation of the sigma_x, sigma_y
  * values. These resulting values are returned via the pmPSFClump structure.
- * 
+ *
  * The return value indicates the success (TRUE) of the operation.
- * 
+ *
  * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD).
  * XXX: Save the clump parameters on the Metadata (TBD)
- * 
+ *
  */
 pmPSFClump pmSourcePSFClump(
@@ -398,5 +405,5 @@
 
 /** pmSourceRoughClass()
- * 
+ *
  * Based on the specified data values, make a guess at the source
  * classification. The sources are provides as a psArray of pmSource entries.
@@ -405,5 +412,5 @@
  * can be extracted from the metadata using the given keywords. Except as noted,
  * the data type for these parameters are psF32.
- * 
+ *
  */
 bool pmSourceRoughClass(
@@ -415,5 +422,5 @@
 
 /** pmSourceModelGuess()
- * 
+ *
  * Convert available data to an initial guess for the given model. This
  * function allocates a pmModel entry for the pmSource structure based on the
@@ -421,5 +428,5 @@
  * are specified for each model below. The guess values are placed in the model
  * parameters. The function returns TRUE on success or FALSE on failure.
- * 
+ *
  */
 pmModel *pmSourceModelGuess(
@@ -430,7 +437,7 @@
 
 /** pmContourType
- * 
+ *
  * Only one type is defined at present.
- * 
+ *
  */
 typedef enum {
@@ -442,5 +449,5 @@
 
 /** pmSourceContour()
- * 
+ *
  * Find points in a contour for the given source at the given level. If type
  * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak,
@@ -452,5 +459,5 @@
  * This function may be used as part of the model guess inputs.  Other contour
  * types may be specified in the future for more refined contours (TBD)
- * 
+ *
  */
 psArray *pmSourceContour(
@@ -463,5 +470,5 @@
 
 /** pmSourceFitModel()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -469,5 +476,5 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel(
@@ -479,12 +486,12 @@
 
 /** pmModelFitStatus()
- * 
+ *
  * This function wraps the call to the model-specific function returned by
  * pmModelFitStatusFunc_GetFunction.  The model-specific function examines the
  * model parameters, parameter errors, Chisq, S/N, and other parameters available
  * from model to decide if the particular fit was successful or not.
- * 
+ *
  * XXX: Must code this.
- * 
+ *
  */
 bool pmModelFitStatus(
@@ -494,5 +501,5 @@
 
 /** pmSourceAddModel()
- * 
+ *
  * Add the given source model flux to/from the provided image. The boolean
  * option center selects if the source is re-centered to the image center or if
@@ -501,5 +508,5 @@
  * is at most the pixel range specified by the source.pixels image. The success
  * status is returned.
- * 
+ *
  */
 bool pmSourceAddModel(
@@ -507,10 +514,11 @@
     psImage *mask,   ///< The image pixel mask (valid == 0)
     pmModel *model,   ///< The input pmModel
-    bool center    ///< A boolean flag that determines whether pixels are centered
+    bool center,    ///< A boolean flag that determines whether pixels are centered
+    bool sky        ///< A boolean flag that determines if the sky is subtracted
 );
 
 
 /** pmSourceSubModel()
- * 
+ *
  * Subtract the given source model flux to/from the provided image. The
  * boolean option center selects if the source is re-centered to the image center
@@ -519,5 +527,5 @@
  * image is at most the pixel range specified by the source.pixels image. The
  * success status is returned.
- * 
+ *
  */
 bool pmSourceSubModel(
@@ -525,10 +533,11 @@
     psImage *mask,   ///< The image pixel mask (valid == 0)
     pmModel *model,   ///< The input pmModel
-    bool center    ///< A boolean flag that determines whether pixels are centered
+    bool center,    ///< A boolean flag that determines whether pixels are centered
+    bool sky        ///< A boolean flag that determines if the sky is subtracted
 );
 
 
 /**
- * 
+ *
  * The function returns both the magnitude of the fit, defined as -2.5log(flux),
  * where the flux is integrated under the model, theoretically from a radius of 0
@@ -537,7 +546,7 @@
  * not excluded by the aperture mask. The model flux is calculated by calling the
  * model-specific function provided by pmModelFlux_GetFunction.
- * 
+ *
  * XXX: must code this.
- * 
+ *
  */
 bool pmSourcePhotometry(
@@ -551,8 +560,8 @@
 
 /**
- * 
+ *
  * This function converts the source classification into the closest available
  * approximation to the Dophot classification scheme:
- * 
+ *
  * PM_SOURCE_DEFECT: 8
  * PM_SOURCE_SATURATED: 8
@@ -569,5 +578,5 @@
  * PM_SOURCE_POOR_FIT_GAL: 2
  * PM_SOURCE_OTHER: ?
- * 
+ *
  */
 int pmSourceDophotType(
@@ -577,11 +586,11 @@
 
 /** pmSourceSextractType()
- * 
+ *
  * This function converts the source classification into the closest available
  * approximation to the Sextractor classification scheme. the correspondence is
  * not yet defined (TBD) .
- * 
+ *
  * XXX: Must code this.
- * 
+ *
  */
 int pmSourceSextractType(
@@ -590,5 +599,5 @@
 
 /** pmSourceFitModel_v5()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -596,5 +605,5 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel_v5(
@@ -606,5 +615,5 @@
 
 /** pmSourceFitModel_v7()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -612,5 +621,5 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel_v7(
@@ -622,7 +631,7 @@
 
 /** pmSourcePhotometry()
- * 
+ *
  * XXX: Need descriptions
- * 
+ *
  */
 bool pmSourcePhotometry(
Index: /trunk/psModules/src/objects/pmPSF.c
===================================================================
--- /trunk/psModules/src/objects/pmPSF.c	(revision 5843)
+++ /trunk/psModules/src/objects/pmPSF.c	(revision 5844)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
Index: /trunk/psModules/src/objects/pmPSFtry.c
===================================================================
--- /trunk/psModules/src/objects/pmPSFtry.c	(revision 5843)
+++ /trunk/psModules/src/objects/pmPSFtry.c	(revision 5844)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -176,5 +176,5 @@
 
     // XXX this function wants aperture radius for pmSourcePhotometry
-    pmPSFtryMetric (psfTry, RADIUS);
+    pmPSFtryMetric_Alt (psfTry, RADIUS);
     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
               modelName,
@@ -288,8 +288,10 @@
 
     // linear fit to rfBin, daBin
-    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
+    psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
+    psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
 
     // XXX EAM : this is the intended API (cycle 7? cycle 8?)
-    poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
+    // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
 
     // XXX EAM : replace this when the above version is implemented
@@ -314,5 +316,78 @@
     psFree (poly);
     psFree (stats);
+    psFree (fitstat);
 
     return true;
 }
+
+bool pmPSFtryMetric_Alt (pmPSFtry *try
+                         , float RADIUS)
+{
+
+    // the measured (aperture - fit) magnitudes (dA == try->metric)
+    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
+    //     dA = dAo + dsky/flux
+    //   where flux is the flux of the star
+    // we fit this trend to find the infinite flux aperture correction (dAo),
+    //   the nominal sky bias (dsky), and the error on dAo
+    // the values of dA are contaminated by stars with close neighbors in the aperture
+    //   we use an outlier rejection to avoid this bias
+
+    // rflux = ten(0.4*fitMag);
+    psVector *rflux = psVectorAlloc (try
+                                     ->sources->n, PS_TYPE_F64);
+    for (int i = 0; i < try
+                ->sources->n; i++) {
+            if (try
+                    ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
+            rflux->data.F64[i] = pow(10.0, 0.4*try
+                                     ->fitMag->data.F64[i]);
+        }
+
+    // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    // XXX EAM
+    stats->min = 1.0;
+    stats->max = 3.0;
+    stats->clipIter = 3;
+
+    // linear clipped fit to rfBin, daBin
+    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+    poly = psVectorClipFitPolynomial1D (poly, stats, try
+                                            ->mask, PSFTRY_MASK_ALL, try
+                                                ->metric, NULL, rflux);
+    fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+
+    try
+        ->psf->ApResid = poly->coeff[0];
+    try
+        ->psf->dApResid = stats->sampleStdev;
+    try
+        ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+    fprintf (stderr, "*******************************************************************************\n");
+
+    FILE *f;
+    f = fopen ("apresid.dat", "w");
+    if (f == NULL)
+        psAbort ("pmPSFtry", "can't open output file");
+
+    for (int i = 0; i < try
+                ->sources->n; i++) {
+            fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try
+                         ->fitMag->data.F64[i], rflux->data.F64[i], try
+                             ->metric->data.F64[i], try
+                                 ->mask->data.U8[i]);
+        }
+    fclose (f);
+
+    psFree (rflux);
+    psFree (poly);
+    psFree (stats);
+
+    // psFree (daFit);
+    // psFree (daResid);
+
+    return true;
+}
Index: /trunk/psModules/src/objects/pmPSFtry.h
===================================================================
--- /trunk/psModules/src/objects/pmPSFtry.h	(revision 5843)
+++ /trunk/psModules/src/objects/pmPSFtry.h	(revision 5844)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 20:32:44 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -18,5 +18,5 @@
 
 /**
- * 
+ *
  * This structure contains a pointer to the collection of sources which will
  * be used to test the PSF model form. It lists the pmModelType type of model
@@ -38,9 +38,9 @@
  * ultimate metric to intercompare multiple types of PSF models is the value of
  * the aperture correction scatter.
- * 
+ *
  * XXX: There are many more members in the SDRS then in the prototype code.
  * I stuck with the prototype code.
- * 
- * 
+ *
+ *
  */
 typedef struct
@@ -58,8 +58,8 @@
 
 /** pmPSFtryMaskValues
- * 
+ *
  * The following datatype defines the masks used by the pmPSFtry analysis to
  * identify sources which should or should not be included in the analysis.
- * 
+ *
  */
 enum {
@@ -74,7 +74,7 @@
 
 /** pmPSFtryAlloc()
- * 
+ *
  * Allocate a pmPSFtry data structure.
- * 
+ *
  */
 pmPSFtry *pmPSFtryAlloc(
@@ -85,9 +85,9 @@
 
 /** pmPSFtryModel()
- * 
+ *
  * This function takes the input collection of sources and performs a complete
  * analysis to determine a PSF model of the given type (specified by model name).
  * The result is a pmPSFtry with the results of the analysis.
- * 
+ *
  */
 pmPSFtry *pmPSFtryModel(
@@ -99,8 +99,8 @@
 
 /** pmPSFtryMetric()
- * 
+ *
  * This function is used to measure the PSF model metric for the set of
  * results contained in the pmPSFtry structure.
- * 
+ *
  */
 bool pmPSFtryMetric(
@@ -109,3 +109,15 @@
 );
 
+/** pmPSFtryMetric_Alt()
+ *
+ * This function is used to measure the PSF model metric for the set of
+ * results contained in the pmPSFtry structure (alternative implementation).
+ *
+ */
+bool pmPSFtryMetric_Alt(
+    pmPSFtry *try
+    ,                      ///< Add comment.
+    float RADIUS                        ///< Add comment.
+);
+
 # endif
Index: /trunk/psModules/test/imsubtract/tst_pmSubtractBias.c
===================================================================
--- /trunk/psModules/test/imsubtract/tst_pmSubtractBias.c	(revision 5843)
+++ /trunk/psModules/test/imsubtract/tst_pmSubtractBias.c	(revision 5844)
@@ -19,5 +19,5 @@
  *  test03: Calculate a row overscan vector and subtract it from each
  *  row in the input image.
- * test05: 
+ * test05:
  *
  *  @author GLG, MHPCC
@@ -25,6 +25,6 @@
  *  XXX: Memory leaks are not being detected.
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-11-23 23:54:30 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:36 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -326,5 +326,5 @@
         rc = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_NONE, stat,
                             0, PM_FIT_NONE, myBias);
-     
+ 
         for (i=0;i<numRows;i++) {
             for (j=0;j<numCols;j++) {
@@ -335,5 +335,5 @@
                     testStatus = 1;
                 }
-     
+ 
                 // Restore myReadout for next test.
                 myReadout->image->data.F32[i][j] = (float) (i + j);
@@ -627,5 +627,5 @@
     // Set overscan axis in the metadata.
     //
-    psBool rc;
+    psBool rc = false;
     if (overscanaxis == PM_OVERSCAN_ROWS) {
         rc = psMetadataAddS32(myReadout->parent->concepts, PS_LIST_HEAD, "CELL.READDIR", 0, NULL, 1);
Index: /trunk/psModules/test/objects/verified/tst_pmObjects01.stdout
===================================================================
--- /trunk/psModules/test/objects/verified/tst_pmObjects01.stdout	(revision 5843)
+++ /trunk/psModules/test/objects/verified/tst_pmObjects01.stdout	(revision 5844)
@@ -1,8 +1,4 @@
 Testing pmPeakAlloc()...
 Testing pmMomentsAlloc()...
-Testing pmModelAlloc(PS_MODEL_GAUSS)...
-Testing pmModelAlloc(PS_MODEL_GAUSS)...
-Testing pmModelAlloc(PS_MODEL_GAUSS)...
-Testing pmModelAlloc(PS_MODEL_GAUSS)...
 ----------------------------------------------------------------------------------
 Calling pmFindVectorPeaks with NULL psVector.  Should generate error and return NULL.