Index: /trunk/psModules/src/objects/Makefile.am
===================================================================
--- /trunk/psModules/src/objects/Makefile.am	(revision 5254)
+++ /trunk/psModules/src/objects/Makefile.am	(revision 5255)
@@ -3,7 +3,20 @@
 libpsmoduleobjects_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS)
 libpsmoduleobjects_la_LDFLAGS  = -release $(PACKAGE_VERSION)
-libpsmoduleobjects_la_SOURCES  = pmObjects.c
+libpsmoduleobjects_la_SOURCES  = \
+    pmObjects.c \
+    pmPSF.c \
+    pmPSFtry.c \
+    pmModelGroup.c \
+    psModulesUtils.c \
+    psLibUtils.c \
+    psEllipse.c
 
 psmoduleincludedir = $(includedir)
 psmoduleinclude_HEADERS = \
-  pmObjects.h
+    pmObjects.h \
+    pmPSF.h \
+    pmPSFtry.h \
+    pmModelGroup.h \
+    psModulesUtils.h \
+    psLibUtils.h \
+    psEllipse.h
Index: /trunk/psModules/src/objects/pmModelGroup.c
===================================================================
--- /trunk/psModules/src/objects/pmModelGroup.c	(revision 5255)
+++ /trunk/psModules/src/objects/pmModelGroup.c	(revision 5255)
@@ -0,0 +1,120 @@
+# include "pmModelGroup.h"
+
+double hypot(double x, double y);
+double sqrt (double x);
+
+#include "psEllipse.h"
+#include "models/pmModel_GAUSS.c"
+#include "models/pmModel_PGAUSS.c"
+#include "models/pmModel_QGAUSS.c"
+#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},
+                               };
+
+pmModelFunc pmModelFunc_GetFunction (pmModelType type)
+{
+    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (NULL);
+    }
+    return (models[type].modelFunc);
+}
+
+pmModelFlux pmModelFlux_GetFunction (pmModelType type)
+{
+    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (NULL);
+    }
+    return (models[type].modelFlux);
+}
+
+pmModelRadius pmModelRadius_GetFunction (pmModelType type)
+{
+    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (NULL);
+    }
+    return (models[type].modelRadius);
+}
+
+pmModelLimits pmModelLimits_GetFunction (pmModelType type)
+{
+    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (NULL);
+    }
+    return (models[type].modelLimits);
+}
+
+pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type)
+{
+    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (NULL);
+    }
+    return (models[type].modelGuessFunc);
+}
+
+pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type)
+{
+    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (NULL);
+    }
+    return (models[type].modelFitStatusFunc);
+}
+
+pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type)
+{
+    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (NULL);
+    }
+    return (models[type].modelFromPSFFunc);
+}
+
+psS32 pmModelParameterCount (pmModelType type)
+{
+    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (0);
+    }
+    return (models[type].nParams);
+}
+
+psS32 pmModelSetType (char *name)
+{
+
+    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
+    for (int i = 0; i < Nmodels; i++) {
+        if (!strcmp(models[i].name, name)) {
+            return (i);
+        }
+    }
+    return (-1);
+}
+
+char *pmModelGetType (pmModelType type)
+{
+
+    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (NULL);
+    }
+    return (models[type].name);
+}
Index: /trunk/psModules/src/objects/pmModelGroup.h
===================================================================
--- /trunk/psModules/src/objects/pmModelGroup.h	(revision 5255)
+++ /trunk/psModules/src/objects/pmModelGroup.h	(revision 5255)
@@ -0,0 +1,211 @@
+/** @file  pmModelGroup.h
+ *
+ *  The object model function types are defined to allow for the flexible addition
+ *  of new object models. Every object model, with parameters represented by
+ *  pmModel, has an associated set of functions which provide necessary support
+ *  operations. A set of abstract functions allow the programmer to select the
+ *  approriate function or property for a specific named object model.
+ * 
+ *  @author EAM, IfA
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-10-10 19:53:40 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ * 
+ */
+#include "pmObjects.h"
+#include "pmPSF.h"
+/**
+ * 
+ * This function returns the number of parameters used by the listed function.
+ * 
+ */
+int pmModelParameterCount(
+    pmModelType type                    ///< Add comment.
+);
+
+
+/**
+ * 
+ * This function returns the user-space model names for the specified model type.
+ * 
+ */
+char *pmModelGetType(
+    pmModelType type                    ///< Add comment.
+);
+
+
+/**
+ * 
+ * This function returns the internal model type code for the user-space model names.
+ * 
+ */
+pmModelType pmModelSetType(
+    char *name                          ///< Add comment.
+);
+
+
+#ifndef PM_MODEL_GROUP_H
+#define PM_MODEL_GROUP_H
+
+/**
+ * 
+ *  This function is the model chi-square minimization function for this model.
+ * 
+ */
+typedef psMinimizeLMChi2Func pmModelFunc;
+
+
+/**
+ * 
+ * This function returns the integrated flux for the given model parameters.
+ */
+typedef psF64 (*pmModelFlux)(const psVector *params);
+
+
+/**
+ * 
+ *  This function returns the radius at which the given model and parameters
+ *  achieves the given flux.
+ * 
+ */
+typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
+
+/**
+ * 
+ *  This function sets the model parameter limits vectors for the given model
+ * 
+ */
+typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max);
+
+/**
+ * 
+ *  This function provides the model guess parameters based on the details of
+ *   the given source.
+ * 
+ */
+typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
+
+
+/**
+ * 
+ *  This function constructs the PSF model for the given source based on the
+ *  supplied psf and the FLT model for the object.
+ * 
+ */
+typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
+
+/**
+ * 
+ *  This function returns the success / failure status of the given model fit
+ * 
+ */
+typedef bool (*pmModelFitStatusFunc)(pmModel *model);
+
+/**
+ * 
+ *  Each of the function types above has a corresponding function which returns
+ *  the function given the model type:
+ * 
+ */
+
+
+/**
+ * 
+ * pmModelFunc is the function used to determine the value of the model at a
+ * specific coordinate, and is the one used by psMinimizeLMChi2.
+ * 
+ */
+pmModelFunc          pmModelFunc_GetFunction (pmModelType type);
+
+
+/**
+ * 
+ * pmModelFlux returns the total integrated flux for the given input parameters.
+ * 
+ */
+pmModelFlux          pmModelFlux_GetFunction (pmModelType type);
+
+
+/**
+ * 
+ * pmModelRadius returns the scaling radius at which the flux of the model
+ * matches the specified flux. This presumes that the model is a function of an
+ * elliptical contour.
+ * 
+ */
+pmModelRadius        pmModelRadius_GetFunction (pmModelType type);
+
+
+/**
+ * 
+ * pmModelLimits sets the parameter limit vectors for the function.
+ * 
+ */
+pmModelLimits        pmModelLimits_GetFunction (pmModelType type);
+
+
+/**
+ * 
+ * pmModelGuessFunc generates an initial guess for the model based on the
+ * provided source statistics (moments and pixel values as needed).
+ * 
+ */
+pmModelGuessFunc     pmModelGuessFunc_GetFunction (pmModelType type);
+
+
+/**
+ * 
+ * pmModelFromPSFFunc takes as input a representation of the psf and a value
+ * for the model and fills in the PSF parameters of the model. The input
+ * primarily relies upon the centroid coordinates of the input model, though the
+ * normalization may potentially be used.
+ * 
+ */
+pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction (pmModelType type);
+
+
+/**
+ * 
+ * pmModelFitStatusFunc returns a true or false values based on the success or
+ * failure of a model fit.  The success is determined by quantities such as the
+ * chisq or the signal-to-noise.
+ * 
+ */
+pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type);
+
+
+
+
+/**
+ * 
+ * Every model instance belongs to a class of models, defined by the value of
+ * the pmModelType type entry. Various functions need access to information about
+ * each of the models. Some of this information varies from model to model, and
+ * may depend on the current parameter values or other data quantities. In order
+ * to keep the code from requiring the information about each model to be coded
+ * into the low-level fitting routines, we define a collection of functions which
+ * allow us to abstract this type of model-dependent information. These generic
+ * functions take the model type and return the corresponding function pointer
+ * for the specified model. Each model is defined by creating this collection of
+ * specific functions, and placing them in a single file for each model. We
+ * define the following structure to carry the collection of information about
+ * the models.
+ * 
+ */
+typedef struct
+{
+    char *name;
+    int nParams;
+    pmModelFunc          modelFunc;
+    pmModelFlux          modelFlux;
+    pmModelRadius        modelRadius;
+    pmModelLimits        modelLimits;
+    pmModelGuessFunc     modelGuessFunc;
+    pmModelFromPSFFunc   modelFromPSFFunc;
+    pmModelFitStatusFunc modelFitStatusFunc;
+}
+pmModelGroup;
+
+# endif
Index: /trunk/psModules/src/objects/pmObjects.c
===================================================================
--- /trunk/psModules/src/objects/pmObjects.c	(revision 5254)
+++ /trunk/psModules/src/objects/pmObjects.c	(revision 5255)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-28 20:43:52 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-10-10 19:53:40 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -13,10 +13,10 @@
  */
 
-#include<stdio.h>
-#include<math.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
 #include "pslib.h"
-#include "psConstants.h"
 #include "pmObjects.h"
-
+#include "pmModelGroup.h"
 /******************************************************************************
 pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members.
@@ -27,4 +27,5 @@
                     pmPeakType class)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak));
     tmp->x = x;
@@ -33,4 +34,5 @@
     tmp->class = class;
 
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmp);
 }
@@ -42,9 +44,10 @@
 pmMoments *pmMomentsAlloc()
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments));
     tmp->x = 0.0;
     tmp->y = 0.0;
     tmp->Sx = 0.0;
-    tmp->Sx = 0.0;
+    tmp->Sy = 0.0;
     tmp->Sxy = 0.0;
     tmp->Sum = 0.0;
@@ -53,4 +56,5 @@
     tmp->nPixels = 0;
 
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmp);
 }
@@ -58,6 +62,8 @@
 static void modelFree(pmModel *tmp)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     psFree(tmp->params);
     psFree(tmp->dparams);
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
 }
 
@@ -65,41 +71,22 @@
 pmModelAlloc(): Allocate the pmModel structure, along with its parameters,
 and initialize the type member.  Initialize the params to 0.0.
-XXX EAM: changing params and dparams to psVector
+XXX EAM: simplifying code with pmModelParameterCount
 *****************************************************************************/
 pmModel *pmModelAlloc(pmModelType type)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel));
 
     tmp->type = type;
     tmp->chisq = 0.0;
-    switch (type) {
-    case PS_MODEL_GAUSS:
-        tmp->params  = psVectorAlloc(7, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(7, PS_TYPE_F32);
-        break;
-    case PS_MODEL_PGAUSS:
-        tmp->params  = psVectorAlloc(7, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(7, PS_TYPE_F32);
-        break;
-    case PS_MODEL_TWIST_GAUSS:
-        tmp->params  = psVectorAlloc(11, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(11, PS_TYPE_F32);
-        break;
-    case PS_MODEL_WAUSS:
-        tmp->params  = psVectorAlloc(9, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(9, PS_TYPE_F32);
-        break;
-    case PS_MODEL_SERSIC:
-        tmp->params  = psVectorAlloc(8, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(8, PS_TYPE_F32);
-        break;
-    case PS_MODEL_SERSIC_CORE:
-        tmp->params  = psVectorAlloc(12, PS_TYPE_F32);
-        tmp->dparams = psVectorAlloc(12, PS_TYPE_F32);
-        break;
-    default:
+    tmp->nIter = 0;
+    psS32 Nparams = pmModelParameterCount(type);
+    if (Nparams == 0) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
         return(NULL);
     }
+
+    tmp->params  = psVectorAlloc(Nparams, PS_TYPE_F32);
+    tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32);
 
     for (psS32 i = 0; i < tmp->params->n; i++) {
@@ -109,19 +96,22 @@
 
     psMemSetDeallocator(tmp, (psFreeFunc) modelFree);
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmp);
 }
 
 /******************************************************************************
-XXX: We don't free pixels and mask since that caused a memory error.
-We might need to increase the reference counter and decrease it here.
+XXX EAM : we can now free these pixels - memory ref is incremented now
 *****************************************************************************/
 static void sourceFree(pmSource *tmp)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     psFree(tmp->peak);
-    //    psFree(tmp->pixels);
-    //    psFree(tmp->mask);
+    psFree(tmp->pixels);
+    psFree(tmp->weight);
+    psFree(tmp->mask);
     psFree(tmp->moments);
     psFree(tmp->modelPSF);
     psFree(tmp->modelFLT);
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
 }
 
@@ -132,7 +122,9 @@
 pmSource *pmSourceAlloc()
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource));
     tmp->peak = NULL;
     tmp->pixels = NULL;
+    tmp->weight = NULL;
     tmp->mask = NULL;
     tmp->moments = NULL;
@@ -142,4 +134,5 @@
     psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
 
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmp);
 }
@@ -159,4 +152,5 @@
                             psF32 threshold)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_VECTOR_NON_NULL(vector, NULL);
     PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL);
@@ -177,4 +171,5 @@
             tmpVector = psVectorAlloc(0, PS_TYPE_U32);
         }
+        psTrace(__func__, 3, "---- %s() end ----\n", __func__);
         return(tmpVector);
     }
@@ -241,4 +236,5 @@
     }
 
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmpVector);
 }
@@ -248,9 +244,11 @@
 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???
 *****************************************************************************/
 static psVector *getRowVectorFromImage(psImage *image,
                                        psU32 row)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
@@ -260,4 +258,5 @@
         tmpVector->data.F32[col] = image->data.F32[row][col];
     }
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return(tmpVector);
 }
@@ -276,4 +275,5 @@
                               pmPeakType type)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);
 
@@ -283,5 +283,8 @@
     }
     psArrayAdd(list, 100, tmpPeak);
-
+    psFree (tmpPeak);
+    // XXX EAM : is this free appropriate?  (does psArrayAdd increment memory counter?)
+
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return(list);
 }
@@ -308,8 +311,10 @@
                           psF32 threshold)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
     if ((image->numRows == 1) || (image->numCols == 1)) {
         psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns.");
+        psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
         return(NULL);
     }
@@ -364,4 +369,6 @@
         }
     }
+    psFree (tmpRow);
+    psFree (row1);
 
     //
@@ -369,4 +376,5 @@
     //
     if (image->numRows == 1) {
+        psTrace(__func__, 3, "---- %s() end ----\n", __func__);
         return(list);
     }
@@ -442,8 +450,10 @@
                 }
             } else {
-                psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");
+                psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
             }
 
         }
+        psFree (tmpRow);
+        psFree (row1);
     }
 
@@ -484,7 +494,10 @@
             }
         } else {
-            psError(PS_ERR_UNKNOWN, true, "peak specified valid colum range.");
-        }
-    }
+            psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
+        }
+    }
+    psFree (tmpRow);
+    psFree (row1);
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(list);
 }
@@ -495,12 +508,13 @@
                              psS32 y)
 {
-
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     if ((x >= valid.x0) &&
             (x <= valid.x1) &&
             (y >= valid.y0) &&
             (y <= valid.y1)) {
+        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
         return(true);
     }
-
+    psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
     return(false);
 }
@@ -516,4 +530,8 @@
  
 XXX: changed API to create a NEW output psArray (should change name as well)
+ 
+XXX: Do we free the psList elements of those culled peaks?
+ 
+XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset?
 *****************************************************************************/
 psList *pmCullPeaks(psList *peaks,
@@ -521,6 +539,6 @@
                     const psRegion valid)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(peaks, NULL);
-    //    PS_ASSERT_PTR_NON_NULL(valid, NULL);
 
     psListElem *tmpListElem = (psListElem *) peaks->head;
@@ -539,4 +557,5 @@
     }
 
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(peaks);
 }
@@ -544,6 +563,11 @@
 // XXX EAM: I changed this to return a new, subset array
 //          rather than alter the existing one
-psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion valid)
-{
+// XXX: Fix the *valid pointer.
+psArray *pmPeaksSubset(
+    psArray *peaks,
+    psF32 maxValue,
+    const psRegion valid)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(peaks, NULL);
 
@@ -561,4 +585,5 @@
         psArrayAdd (output, 200, tmpPeak);
     }
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(output);
 }
@@ -601,158 +626,51 @@
      members.
 *****************************************************************************/
-pmSource *pmSourceLocalSky(const psImage *image,
-                           const pmPeak *peak,
-                           psStatsOptions statsOptions,
-                           psF32 innerRadius,
-                           psF32 outerRadius)
-{
-    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
-    PS_ASSERT_PTR_NON_NULL(peak, NULL);
-    PS_FLOAT_COMPARE(0.0, innerRadius, NULL);
-    PS_FLOAT_COMPARE(innerRadius, outerRadius, NULL);
-    psS32 innerRadiusS32 = (psS32) innerRadius;
-    psS32 outerRadiusS32 = (psS32) outerRadius;
-
-    //
-    // We define variables for code readability.
-    //
-    // XXX: Since the peak->xy coords are in image, not subImage coords,
-    // these variables should be renamed for clarity (imageCenterRow, etc).
-    //
-    // peak->x,y is guaranteed to be on image
-    psS32 SubImageCenterRow = peak->y;
-    psS32 SubImageCenterCol = peak->x;
-
-    // XXX EAM : I added this code to stay on the image. So did George
-    psS32 SubImageStartRow  = PS_MAX(0, SubImageCenterRow - outerRadiusS32);
-    psS32 SubImageEndRow    = PS_MIN(image->numRows - 1, SubImageCenterRow + outerRadiusS32);
-    psS32 SubImageStartCol  = PS_MAX(0, SubImageCenterCol - outerRadiusS32);
-    psS32 SubImageEndCol    = PS_MIN(image->numCols - 1, SubImageCenterCol + outerRadiusS32);
-    // AnulusWidth == number of pixels width in the annulus.  We add one since
-    // the pixels at the inner AND outher radius are included.
-    psS32 AnulusWidth = 1 + (outerRadiusS32 - innerRadiusS32);
-    // Example: assume an outer/inner radius of 20/10.  Then the subimage
-    // should have width/length of 40.  An 18-by-18 interior region will
-    // be masked.
-    //    printf("pmSourceLocalSky(): innerRadiusS32 is %d\n", innerRadiusS32);
-    //    printf("pmSourceLocalSky(): outerRadiusS32 is %d\n", outerRadiusS32);
-    //    printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth);
-
-    // XXX EAM : these tests should not be needed: we can never hit this error because of above
-    # if (1)
-
-        if (SubImageStartRow < 0) {
-            psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
-                    SubImageStartRow);
-            return(NULL);
-        }
-    if (SubImageEndRow >= image->numRows) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
-                SubImageEndRow);
-        return(NULL);
-    }
-    if (SubImageStartCol < 0) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",
-                SubImageStartCol);
-        return(NULL);
-    }
-    if (SubImageEndCol >= image->numCols) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
-                SubImageEndCol);
-        return(NULL);
-    }
-    # endif
-
-    //
-    // Grab a subimage of the original image of size (2 * outerRadius).
-    //
-    // XXX: Must fix for new psImageSubset
-    //    psImage *subImage = psImageSubset((psImage *) image,
-    //                                      SubImageStartCol,
-    //                                      SubImageStartRow,
-    //                                      SubImageEndCol,
-    //                                      SubImageEndRow);
-    //    printf("pmSourceLocalSky: subimage width/length is (%d, %d)\n", subImage->numCols, subImage->numRows);
-    psRegion tmpRegion = psRegionSet(SubImageStartCol,
-                                     SubImageEndCol,
-                                     SubImageStartRow,
-                                     SubImageEndRow);
-    psImage *subImage = psImageSubset((psImage *) image, tmpRegion);
-
-    psImage *subImageMask = psImageAlloc(subImage->numCols,
-                                         subImage->numRows,
-                                         PS_TYPE_U8);
-
-    //
-    // Loop through the subimage mask, initialize mask to 0.
-    //
-    for (psS32 row = 0 ; row < subImageMask->numRows; row++) {
-        for (psS32 col = 0 ; col < subImageMask->numCols; col++) {
-            subImageMask->data.U8[row][col] = 0;
-        }
-    }
-
-    //
-    // Loop through the subimage, mask off pixels in the inner square.
-    // XXX this uses a static mask value of 1
-    //
-    for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) {
-        for (psS32 col = AnulusWidth; col <= (subImageMask->numCols - AnulusWidth) - 1; col++) {
-            subImageMask->data.U8[row][col] = 1;
-        }
-    }
-
-
-    //    for (psS32 row = 0 ; row < subImage->numRows; row++) {
-    //        for (psS32 col = 0 ; col < subImage->numCols; col++) {
-    //            printf("(%d) ", subImageMask->data.U8[row][col]);
-    //        }
-    //        printf("\n");
-    //    }
-
-    //
-    // Allocate the myStats structure, then call psImageStats(), which will
-    // calculate the specified statistic.
-    //
+
+
+
+
+
+
+
+
+
+bool pmSourceLocalSky(
+    pmSource *source,
+    psStatsOptions statsOptions,
+    psF32 Radius)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->pixels, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_INT_POSITIVE(Radius, false);
+    PS_ASSERT_INT_NONNEGATIVE(Radius, false);
+
+    psImage *image = source->pixels;
+    psImage *mask  = source->mask;
+    pmPeak *peak  = source->peak;
+    psRegion srcRegion;
+
+    srcRegion = psRegionForSquare(peak->x, peak->y, Radius);
+    srcRegion = psRegionForImage(mask, srcRegion);
+
+    psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED);
     psStats *myStats = psStatsAlloc(statsOptions);
-    myStats = psImageStats(myStats, subImage, subImageMask, 1);
-
-    //
-    // Create the output mySource, and set appropriate members.
-    //
-    pmSource *mySource = pmSourceAlloc();
-    mySource->peak = (pmPeak *) peak;
-    mySource->moments = pmMomentsAlloc();
+    myStats = psImageStats(myStats, image, mask, 0xff);
+    psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED);
+
     psF64 tmpF64;
     p_psGetStatValue(myStats, &tmpF64);
-    mySource->moments->Sky = (psF32) tmpF64;
-    mySource->pixels = subImage;
-    mySource->mask = subImageMask;
-
-    //
-    // Free things.  XXX: This should be static memory.
-    //
     psFree(myStats);
 
-    return(mySource);
-}
-
-/******************************************************************************
-bool checkRadius(*peak, radius, x, y): private function which simply
-determines if the (x, y) point is within the radius of the specified peak.
- 
-XXX: macro this for performance.
-*****************************************************************************/
-static bool checkRadius(pmPeak *peak,
-                        psF32 radius,
-                        psS32 x,
-                        psS32 y)
-{
-    if (PS_SQR(radius) >= (psF32) (PS_SQR(x - peak->x) + PS_SQR(y - peak->y))) {
-        return(true);
-    }
-
-    return(false);
+    if (isnan(tmpF64)) {
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return(false);
+    }
+    source->moments = pmMomentsAlloc();
+    source->moments->Sky = (psF32) tmpF64;
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
+    return (true);
 }
 
@@ -770,4 +688,5 @@
                          psF32 y)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     /// XXX EAM should compare with hypot (x,y) for speed
     if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) {
@@ -775,4 +694,5 @@
     }
 
+    psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
     return(false);
 }
@@ -787,16 +707,23 @@
     pmSource->peak
     pmSource->pixels
+    pmSource->weight
+    pmSource->mask
  
 XXX: The peak calculations are done in image coords, not subImage coords.
  
-XXX: mask values?
-*****************************************************************************/
-pmSource *pmSourceMoments(pmSource *source,
-                          psF32 radius)
-{
-    PS_ASSERT_PTR_NON_NULL(source, NULL);
-    PS_ASSERT_PTR_NON_NULL(source->peak, NULL);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, NULL);
-    PS_FLOAT_COMPARE(0.0, radius, NULL);
+XXX EAM : this version clips input pixels on S/N
+XXX EAM : this version returns false for several reasons
+*****************************************************************************/
+# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
+
+bool pmSourceMoments(pmSource *source,
+                     psF32 radius)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
     //
@@ -816,7 +743,7 @@
     // XY  = SUM (x - xc)*(y - yc)*(z - sky)
     //
-    psF32 Sum = 0.0;
     psF32 peakPixel = -PS_MAX_F32;
     psS32 numPixels = 0;
+    psF32 Sum = 0.0;
     psF32 X1 = 0.0;
     psF32 Y1 = 0.0;
@@ -824,7 +751,11 @@
     psF32 Y2 = 0.0;
     psF32 XY = 0.0;
-    psF32 x = 0;
-    psF32 y = 0;
-    //
+    psF32 x  = 0;
+    psF32 y  = 0;
+    psF32 R2 = PS_SQR(radius);
+
+    psF32 xPeak = source->peak->x;
+    psF32 yPeak = source->peak->y;
+
     // XXX why do I get different results for these two methods of finding Sx?
     // XXX Sx, Sy would be better measured if we clip pixels close to sky
@@ -834,34 +765,48 @@
     // proceed with the moments calculation.  need to do two loops for a
     // numerically stable result.  first loop: get the sums.
-    //
+    // XXX EAM : mask == 0 is valid
+
     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] != 0)) {
-                psS32 imgColCoord = col + source->pixels->col0;
-                psS32 imgRowCoord = row + source->pixels->row0;
-                if (checkRadius(source->peak,
-                                radius,
-                                imgColCoord,
-                                imgRowCoord)) {
-                    psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
-                    psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
-                    psF32 pDiff = source->pixels->data.F32[row][col] - sky;
-
-                    Sum+= pDiff;
-                    X1+= xDiff * pDiff;
-                    Y1+= yDiff * pDiff;
-                    XY+= xDiff * yDiff * pDiff;
-
-                    X2+= PS_SQR(xDiff) * pDiff;
-                    Y2+= PS_SQR(yDiff) * pDiff;
-
-                    if (source->pixels->data.F32[row][col] > peakPixel) {
-                        peakPixel = source->pixels->data.F32[row][col];
-                    }
-                    numPixels++;
-                }
-            }
-        }
-    }
+            if ((source->mask != NULL) && (source->mask->data.U8[row][col]))
+                continue;
+
+            psF32 xDiff = col + source->pixels->col0 - xPeak;
+            psF32 yDiff = row + source->pixels->row0 - yPeak;
+
+            // XXX EAM : calculate xDiff, yDiff up front;
+            //           radius is just a function of (xDiff, yDiff)
+            if (!VALID_RADIUS(xDiff, yDiff, R2))
+                continue;
+
+            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
+
+            // 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)
+                continue;
+
+            Sum += pDiff;
+            X1  += xDiff * pDiff;
+            Y1  += yDiff * pDiff;
+            XY  += xDiff * yDiff * pDiff;
+
+            X2  += PS_SQR(xDiff) * pDiff;
+            Y2  += PS_SQR(yDiff) * pDiff;
+
+            peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel);
+            numPixels++;
+        }
+    }
+    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
+    if ((numPixels < 3) || (Sum <= 0)) {
+        psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return (false);
+    }
+
+    psTrace (".psModules.pmSourceMoments", 5,
+             "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
+             sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
 
     //
@@ -872,8 +817,17 @@
     x = X1/Sum;
     y = Y1/Sum;
-    source->moments->x = x + ((psF32) source->peak->x);
-    source->moments->y = y + ((psF32) source->peak->y);
-
-    source->moments->Sxy = XY/Sum;
+    if ((fabs(x) > radius) || (fabs(y) > radius)) {
+        psTrace (".psModules.pmSourceMoments", 5,
+                 "large centroid swing; invalid peak %d, %d\n",
+                 source->peak->x, source->peak->y);
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return (false);
+    }
+
+    source->moments->x = x + xPeak;
+    source->moments->y = y + yPeak;
+
+    // XXX EAM : Sxy needs to have x*y subtracted
+    source->moments->Sxy = XY/Sum - x*y;
     source->moments->Sum = Sum;
     source->moments->Peak = peakPixel;
@@ -883,40 +837,12 @@
     source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
     source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
-    return(source);
-
-    // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it
-    # if (0)
-        //
-        // second loop: get the difference sums
-        //
-        X2 = Y2 = 0;
-    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] != 0)) {
-                psS32 imgColCoord = col + source->pixels->col0;
-                psS32 imgRowCoord = row + source->pixels->row0;
-                if (checkRadius(source->peak,
-                                radius,
-                                imgColCoord,
-                                imgRowCoord)) {
-                    psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
-                    psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
-                    psF32 pDiff = source->pixels->data.F32[row][col] - sky;
-
-                    Sum+= pDiff;
-                    X2+= PS_SQR(xDiff - x) * pDiff;
-                    Y2+= PS_SQR(yDiff - y) * pDiff;
-                }
-            }
-        }
-    }
-
-    //
-    // second moment X = sqrt (X2/Sum)
-    //
-    source->moments->Sx = (X2/Sum);
-    source->moments->Sy = (Y2/Sum);
-    return(source);
-    # endif
+
+    psTrace (".psModules.pmSourceMoments", 4,
+             "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
+             sky, Sum, source->moments->x, source->moments->y,
+             source->moments->Sx, source->moments->Sy, source->moments->Sxy);
+
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
+    return(true);
 }
 
@@ -924,4 +850,5 @@
 int pmComparePeakAscend (const void **a, const void **b)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmPeak *A = *(pmPeak **)a;
     pmPeak *B = *(pmPeak **)b;
@@ -930,8 +857,12 @@
 
     diff = A->counts - B->counts;
-    if (diff < FLT_EPSILON)
+    if (diff < FLT_EPSILON) {
+        psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
         return (-1);
-    if (diff > FLT_EPSILON)
+    } else if (diff > FLT_EPSILON) {
+        psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
         return (+1);
+    }
+    psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);
     return (0);
 }
@@ -939,4 +870,5 @@
 int pmComparePeakDescend (const void **a, const void **b)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     pmPeak *A = *(pmPeak **)a;
     pmPeak *B = *(pmPeak **)b;
@@ -945,37 +877,34 @@
 
     diff = A->counts - B->counts;
-    if (diff < FLT_EPSILON)
+    if (diff < FLT_EPSILON) {
+        psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
         return (+1);
-    if (diff > FLT_EPSILON)
+    } else if (diff > FLT_EPSILON) {
+        psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
         return (-1);
+    }
+    psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);
     return (0);
 }
 
 /******************************************************************************
-pmSourceRoughClass(source, metadata): make a guess at the source
-classification.
- 
-XXX: This is not useable code, as of the release date.  There remains a fair
-bit of coding to be completed.
- 
-XXX: The sigX and sigY stuff in the SDRS is unclear.
- 
-XXX: How can this function ever return FALSE?
-*****************************************************************************/
-
-# define NPIX 10
-# define SCALE 0.1
-
-// XXX I am ignore memory freeing issues (EAM)
-bool pmSourceRoughClass(psArray *sources, psMetadata *metadata)
-{
-    PS_ASSERT_PTR_NON_NULL(sources, false);
-    PS_ASSERT_PTR_NON_NULL(metadata, false);
-    psBool rc = true;
+pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 
+sigma-x, sigma-y plane. return 0,0 clump in case of error. 
+*****************************************************************************/
+
+// XXX EAM include a S/N cutoff in selecting the sources?
+pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+
+    # define NPIX 10
+    # define SCALE 0.1
+
     psArray *peaks  = NULL;
-    psF32 clumpX = 0.0;
-    psF32 clumpDX = 0.0;
-    psF32 clumpY = 0.0;
-    psF32 clumpDY = 0.0;
+    pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0};
+    pmPSFClump psfClump = emptyClump;
+
+    PS_ASSERT_PTR_NON_NULL(sources, emptyClump);
+    PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);
 
     // find the sigmaX, sigmaY clump
@@ -986,5 +915,10 @@
 
         // construct a sigma-plane image
+        // psImageAlloc does zero the data
         splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
+        for (int i = 0; i < splane->numRows; i++)
+        {
+            memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
+        }
 
         // place the sources in the sigma-plane image (ignore 0,0 values?)
@@ -992,6 +926,10 @@
         {
             pmSource *tmpSrc = (pmSource *) sources->data[i];
-            PS_ASSERT_PTR_NON_NULL(tmpSrc, false); // just skip this one?
-            PS_ASSERT_PTR_NON_NULL(tmpSrc->moments, false); // just skip this one?
+            if (tmpSrc == NULL) {
+                continue;
+            }
+            if (tmpSrc->moments == NULL) {
+                continue;
+            }
 
             // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
@@ -1022,5 +960,11 @@
         psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);
 
-    }
+        psFree (splane);
+        psFree (stats);
+
+    }
+    // XXX EAM : possible errors:
+    //           1) no peak in splane
+    //           2) no significant peak in splane
 
     // measure statistics on Sx, Sy if Sx, Sy within range of clump
@@ -1036,5 +980,5 @@
         psArraySort (peaks, pmComparePeakDescend);
         clump = peaks->data[0];
-        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d\n", clump->x, clump->y);
+        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);
 
         // define section window for clump
@@ -1077,24 +1021,58 @@
 
         stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);
-        clumpX  = stats->clippedMean;
-        clumpDX = stats->clippedStdev;
+        psfClump.X  = stats->clippedMean;
+        psfClump.dX = stats->clippedStdev;
 
         stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);
-        clumpY  = stats->clippedMean;
-        clumpDY = stats->clippedStdev;
-
-        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", clumpX, clumpY);
-        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", clumpDX, clumpDY);
+        psfClump.Y  = stats->clippedMean;
+        psfClump.dY = stats->clippedStdev;
+
+        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
+        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
         // these values should be pushed on the metadata somewhere
-    }
-
-    int Nsat   = 0;
-    int Ngal   = 0;
-    int Nfaint = 0;
-    int Nstar  = 0;
-    int Npsf   = 0;
-    int Ncr    = 0;
+
+        psFree (stats);
+        psFree (peaks);
+        psFree (tmpSx);
+        psFree (tmpSy);
+    }
+
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
+    return (psfClump);
+}
+
+/******************************************************************************
+pmSourceRoughClass(source, metadata): make a guess at the source
+classification.
+ 
+XXX: push the clump info into the metadata?
+ 
+XXX: How can this function ever return FALSE?
+ 
+XXX EAM : add the saturated mask value to metadata 
+*****************************************************************************/
+
+bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+
+    psBool rc = true;
+
+    int Nsat     = 0;
+    int Ngal     = 0;
+    int Nstar    = 0;
+    int Npsf     = 0;
+    int Ncr      = 0;
+    int Nsatstar = 0;
+
     psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
     starsn->n = 0;
+
+    // check return status value (do these exist?)
+    bool status;
+    psF32 RDNOISE    = psMetadataLookupF32 (&status, metadata, "RDNOISE");
+    psF32 GAIN       = psMetadataLookupF32 (&status, metadata, "GAIN");
+    psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
+    // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
 
     // XXX allow clump size to be scaled relative to sigmas?
@@ -1109,30 +1087,48 @@
         psF32 sigY = tmpSrc->moments->Sy;
 
-        // check return status value (do these exist?)
-        bool status;
-        psF32 RDNOISE  = psMetadataLookupF32 (&status, metadata, "RDNOISE");
-        psF32 GAIN     = psMetadataLookupF32 (&status, metadata, "GAIN");
-        psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
-
-        psF32 PSF_SN_LIM   = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
-        psF32 FAINT_SN_LIM = psMetadataLookupF32 (&status, metadata, "FAINT_SN_LIM");
-
-        // saturated object (star or single pixel not distinguished)
-        if (tmpSrc->moments->Peak > SATURATE) {
-            tmpSrc->type |= PS_SOURCE_SATURATED;
+        // calculate and save signal-to-noise estimates
+        psF32 S  = tmpSrc->moments->Sum;
+        psF32 A  = 4 * M_PI * sigX * sigY;
+        psF32 B  = tmpSrc->moments->Sky;
+        psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));
+        psF32 SN = (S * sqrt(GAIN) / RT);
+        tmpSrc->moments->SN = SN;
+
+        // 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);
+
+        // saturated star (size consistent with PSF or larger)
+        // Nsigma should be user-configured parameter
+        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
+        if ((Nsatpix > 1) && big) {
+            tmpSrc->type |= PM_SOURCE_SATSTAR;
+            Nsatstar ++;
+            continue;
+        }
+
+        // saturated object (not a star, eg bleed trails, hot pixels)
+        if (Nsatpix > 1) {
+            tmpSrc->type |= PM_SOURCE_SATURATED;
             Nsat ++;
             continue;
         }
 
-        // too small to be stellar
-        if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) {
-            tmpSrc->type |= PS_SOURCE_DEFECT;
+        // likely defect (too small to be stellar) (push out to 3 sigma)
+        // low S/N objects which are small are probably stellar
+        // only set candidate defects if
+        if ((sigX < 0.05) || (sigY < 0.05)) {
+            tmpSrc->type |= PM_SOURCE_DEFECT;
             Ncr ++;
             continue;
         }
 
-        // possible galaxy
-        if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) {
-            tmpSrc->type |= PS_SOURCE_GALAXY;
+        // likely unsaturated galaxy (too large to be stellar)
+        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
+            tmpSrc->type |= PM_SOURCE_GALAXY;
             Ngal ++;
             continue;
@@ -1140,24 +1136,12 @@
 
         // the rest are probable stellar objects
-        psF32 S  = tmpSrc->moments->Sum;
-        psF32 A  = M_PI * sigX * sigY;
-        psF32 B  = tmpSrc->moments->Sky;
-        psF32 RT = sqrtf(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrtf(GAIN)));
-        psF32 SN = (S * sqrtf(GAIN) / RT);
-
         starsn->data.F32[starsn->n] = SN;
         starsn->n ++;
         Nstar ++;
 
-        // faint star
-        if (SN < FAINT_SN_LIM) {
-            tmpSrc->type |= PS_SOURCE_FAINTSTAR;
-            Nfaint ++;
-            continue;
-        }
-
-        // PSF star
-        if (SN > PSF_SN_LIM) {
-            tmpSrc->type |= PS_SOURCE_PSFSTAR;
+        // PSF star (within 1.5 sigma of clump center
+        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
+        if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
+            tmpSrc->type |= PM_SOURCE_PSFSTAR;
             Npsf ++;
             continue;
@@ -1165,5 +1149,5 @@
 
         // random type of star
-        tmpSrc->type |= PS_SOURCE_OTHER;
+        tmpSrc->type |= PM_SOURCE_OTHER;
     }
 
@@ -1173,26 +1157,53 @@
         stats = psVectorStats (stats, starsn, NULL, NULL, 0);
         psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
-    }
-
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:  %3d\n", Nstar);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:   %3d\n", Npsf);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nfaint: %3d\n", Nfaint);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:   %3d\n", Ngal);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:   %3d\n", Nsat);
-    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:    %3d\n", Ncr);
-
+        psFree (stats);
+        psFree (starsn);
+    }
+
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:    %3d\n", Nstar);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:     %3d\n", Npsf);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:     %3d\n", Ngal);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:     %3d\n", Nsat);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:      %3d\n", Ncr);
+
+    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
 }
 
+/** 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
+ * a square window (of full width 2Radius+1) centered on the pixel which contains
+ * the given coordinate, in the frame of the readout. The window is defined to
+ * have limits which are valid within the boundary of the readout image, thus if
+ * the radius would fall outside the image pixels, the subimage is truncated to
+ * only consist of valid pixels. If readout->mask or readout->weight are not
+ * NULL, matching subimages are defined for those images as well. This function
+ * fails if no valid pixels can be defined (x or y less than Radius, for
+ * 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(
+    pmSource *mySource,                 ///< Add comment.
+    pmReadout *readout,                 ///< Add comment.
+    psF32 x,                            ///< Add comment.
+    psF32 y,                            ///< Add comment.
+    psF32 Radius)                       ///< Add comment.
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented.  Returning FALSE.\n");
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
+    return(false);
+}
+
 /******************************************************************************
 pmSourceSetPixelsCircle(source, image, radius)
  
-XXX: Why boolean output?
- 
-XXX: Why are we checking source->moments for NULL?  Should the circle be
-     centered on the centroid or the peak?
- 
-XXX: The circle will have a diameter of (1+radius).  This is different from
-     the pmSourceSetLocal() function.
+XXX: This was replaced by DefinePixels in SDRS.  Remove it.
 *****************************************************************************/
 bool pmSourceSetPixelsCircle(pmSource *source,
@@ -1200,9 +1211,10 @@
                              psF32 radius)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_IMAGE_NON_NULL(image, false);
     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(source->moments, false);
-    // PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_FLOAT_COMPARE(0.0, radius, false);
 
@@ -1217,33 +1229,9 @@
     psS32 SubImageCenterCol = source->peak->x;
     // XXX EAM : for the circle to stay on the image
+    // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
     psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - radiusS32);
-    psS32 SubImageEndRow    = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32);
+    psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1);
     psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
-    psS32 SubImageEndCol    = PS_MIN (image->numCols - 1, SubImageCenterCol + radiusS32);
-
-    // XXX EAM : this should not be needed: we can never hit this error
-    # if (1)
-
-        if (SubImageStartRow < 0) {
-            psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
-                    SubImageStartRow);
-            return(false);
-        }
-    if (SubImageEndRow >= image->numRows) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
-                SubImageEndRow);
-        return(false);
-    }
-    if (SubImageStartCol < 0) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",
-                SubImageStartCol);
-        return(false);
-    }
-    if (SubImageEndCol >= image->numCols) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
-                SubImageEndCol);
-        return(false);
-    }
-    # endif
+    psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);
 
     // XXX: Must recycle image.
@@ -1255,10 +1243,8 @@
         psFree(source->pixels);
     }
-    // XXX: Must fix this.  psImageSubset() has different parameters in latest CVS.
-    //    source->pixels = psImageSubset((psImage *) image,
-    //                                   SubImageStartCol,
-    //                                   SubImageStartRow,
-    //                                   SubImageEndCol,
-    //                                   SubImageEndRow);
+    source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
+                                   SubImageStartRow,
+                                   SubImageEndCol,
+                                   SubImageEndRow));
 
     // XXX: Must recycle image.
@@ -1268,5 +1254,5 @@
     source->mask = psImageAlloc(source->pixels->numCols,
                                 source->pixels->numRows,
-                                PS_TYPE_F32);
+                                PS_TYPE_U8); // XXX EAM : type was F32
 
     //
@@ -1287,135 +1273,31 @@
         }
     }
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     return(true);
 }
 
 /******************************************************************************
-pmSourceModelGuess(source, image, model): This function allocates a new
-pmModel structure and stores it in the pmSource data structure specified in
-the argument list.  The model type is specified in the argument list.  The
-params array in that pmModel structure are allocated, and then set to the
-appropriate values.  This function returns true if everything was successful.
- 
-XXX: Many of the initial parameters are set to 0.0 since I don't know what
-the appropiate initial guesses are.
- 
-XXX: The image argument is redundant.
+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.  
  
 XXX: Many parameters are based on the src->moments structure, which is in
 image, not subImage coords.  Therefore, the calls to the model evaluation
 functions will be in image, not subImage coords.  Remember this.
- 
-XXX: The source->models member used to be allocated here.  Now I assume
-->modelPSF should be allocated
-*****************************************************************************/
-bool pmSourceModelGuess(pmSource *source,
-                        const psImage *image,
-                        pmModelType model)
-{
-    PS_ASSERT_PTR_NON_NULL(source, false);
+*****************************************************************************/
+pmModel *pmSourceModelGuess(pmSource *source,
+                            pmModelType modelType)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(source->moments, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_IMAGE_NON_NULL(image, false);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
-    if (source->modelPSF != NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->modelPSF was non-NULL; calling psFree(source->modelPSF).\n");
-        psFree(source->modelPSF);
-    }
-    if (!((model == PS_MODEL_GAUSS) ||
-            (model == PS_MODEL_PGAUSS) ||
-            (model == PS_MODEL_WAUSS) ||
-            (model == PS_MODEL_TWIST_GAUSS) ||
-            (model == PS_MODEL_SERSIC) ||
-            (model == PS_MODEL_SERSIC_CORE))) {
-        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
-        return(false);
-    }
-
-    source->modelPSF = pmModelAlloc(model);
-
-    psVector *params = source->modelPSF->params;
-
-    switch (model) {
-    case PS_MODEL_GAUSS:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
-        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
-        params->data.F32[6] = source->moments->Sxy;
-        return(true);
-
-    case PS_MODEL_PGAUSS:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
-        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
-        params->data.F32[6] = source->moments->Sxy;
-        return(true);
-
-    case PS_MODEL_WAUSS:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
-        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
-        params->data.F32[6] = source->moments->Sxy;
-        // XXX: What are these?
-        // source->modelPSF->params[7] = B2;
-        // source->modelPSF->params[8] = B3;
-        return(true);
-
-        // XXX EAM : I might drop this model (or rather, replace it)
-    case PS_MODEL_TWIST_GAUSS:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        // XXX: What are these?
-        // params->data.F32[4] = SxInner;
-        // params->data.F32[5] = SyInner;
-        // params->data.F32[6] = SxyInner;
-        // params->data.F32[7] = SxOuter;
-        // params->data.F32[8] = SyOuter;
-        // params->data.F32[9] = SxyOuter;
-        // params->data.F32[10] = N;
-        return(true);
-
-    case PS_MODEL_SERSIC:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
-        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
-        params->data.F32[6] = source->moments->Sxy;
-        // XXX: What are these?
-        //params->data.F32[7] = Nexp;
-        return(true);
-
-    case PS_MODEL_SERSIC_CORE:
-        params->data.F32[0] = source->moments->Sky;
-        params->data.F32[1] = source->peak->counts - source->moments->Sky;
-        params->data.F32[2] = source->moments->x;
-        params->data.F32[3] = source->moments->y;
-        // XXX: What are these?
-        // params->data.F32[4] SxInner;
-        // params->data.F32[5] SyInner;
-        // params->data.F32[6] SxyInner;
-        // params->data.F32[7] Zd;
-        // params->data.F32[8] SxOuter;
-        // params->data.F32[9] SyOuter;
-        // params->data.F32[10] = SxyOuter;
-        // params->data.F32[11] = Nexp;
-        return(true);
-
-    default:
-        psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
-        return(false);
-    }
+
+    pmModel *model = pmModelAlloc(modelType);
+
+    pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType);
+    modelGuessFunc(model, source);
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
+    return(model);
 }
 
@@ -1441,50 +1323,27 @@
 testing.  Try to reproduce that and debug.
 *****************************************************************************/
-static psF32 evalModel(pmSource *src,
-                       psU32 row,
-                       psU32 col)
-{
-    PS_ASSERT_PTR_NON_NULL(src, false);
-    PS_ASSERT_PTR_NON_NULL(src->modelPSF, false);
-    PS_ASSERT_PTR_NON_NULL(src->modelPSF->params, false);
-
-    // XXX: The following step will not be necessary if the modelPSF->params
-    // member is a psVector.  Suggest to IfA.
-
-    // XXX EAM: done: modelPSF->params is now a vector
-    psVector *params = src->modelPSF->params;
-
-    //
+
+// XXX EAM : I have made this a public function
+// XXX EAM : this now uses a pmModel as the input
+// XXX EAM : it was using src->type to find the model, not model->type
+psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(image, false);
+    PS_ASSERT_PTR_NON_NULL(model, false);
+    PS_ASSERT_PTR_NON_NULL(model->params, false);
+
     // Allocate the x coordinate structure and convert row/col to image space.
     //
     psVector *x = psVectorAlloc(2, PS_TYPE_F32);
-    x->data.F32[0] = (psF32) (col + src->pixels->col0);
-    x->data.F32[1] = (psF32) (row + src->pixels->row0);
+    x->data.F32[0] = (psF32) (col + image->col0);
+    x->data.F32[1] = (psF32) (row + image->row0);
     psF32 tmpF;
-
-    switch (src->modelPSF->type) {
-    case PS_MODEL_GAUSS:
-        tmpF = pmMinLM_Gauss2D(NULL, params, x);
-        break;
-    case PS_MODEL_PGAUSS:
-        tmpF = pmMinLM_PsuedoGauss2D(NULL, params, x);
-        break;
-    case PS_MODEL_TWIST_GAUSS:
-        tmpF = pmMinLM_TwistGauss2D(NULL, params, x);
-        break;
-    case PS_MODEL_WAUSS:
-        tmpF = pmMinLM_Wauss2D(NULL, params, x);
-        break;
-    case PS_MODEL_SERSIC:
-        tmpF = pmMinLM_Sersic(NULL, params, x);
-        break;
-    case PS_MODEL_SERSIC_CORE:
-        tmpF = pmMinLM_SersicCore(NULL, params, x);
-        break;
-    default:
-        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
-        return(NAN);
-    }
+    pmModelFunc modelFunc;
+
+    modelFunc = pmModelFunc_GetFunction (model->type);
+    tmpF = modelFunc (NULL, model->params, x);
     psFree(x);
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmpF);
 }
@@ -1508,4 +1367,5 @@
                        psU32 dir)
 {
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     //
     // Convert coords to subImage space.
@@ -1517,13 +1377,18 @@
     if (!((0 <= subCol) && (subCol < source->pixels->numCols))) {
         psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range");
+        psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
         return(NAN);
     }
     if (!((0 <= subRow) && (subRow < source->pixels->numRows))) {
+        psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
         psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range");
         return(NAN);
     }
 
-    psF32 oldValue = evalModel(source, subRow, subCol);
+    // XXX EAM : i changed this to match pmModelEval above, but see
+    // XXX EAM   the note below in pmSourceContour
+    psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
     if (oldValue == level) {
+        psTrace(__func__, 4, "---- %s() end ----\n", __func__);
         return(((psF32) (subCol + source->pixels->col0)));
     }
@@ -1545,6 +1410,7 @@
 
     while (subCol != lastColumn) {
-        psF32 newValue = evalModel(source, subRow, subCol);
+        psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow);
         if (oldValue == level) {
+            psTrace(__func__, 4, "---- %s() end ----\n", __func__);
             return((psF32) (subCol + source->pixels->col0));
         }
@@ -1552,4 +1418,5 @@
         if ((newValue <= level) && (level <= oldValue)) {
             // This is simple linear interpolation.
+            psTrace(__func__, 4, "---- %s() end ----\n", __func__);
             return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) );
         }
@@ -1557,4 +1424,5 @@
         if ((oldValue <= level) && (level <= newValue)) {
             // This is simple linear interpolation.
+            psTrace(__func__, 4, "---- %s() end ----\n", __func__);
             return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) );
         }
@@ -1563,4 +1431,5 @@
     }
 
+    psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
     return(NAN);
 }
@@ -1576,4 +1445,6 @@
 XXX: What is mode?
 XXX: The top, bottom of the contour is not correctly determined.
+XXX EAM : this functions is using the model for the contour, but it should
+          be using only the image counts
 *****************************************************************************/
 psArray *pmSourceContour(pmSource *source,
@@ -1582,4 +1453,5 @@
                          pmContourType mode)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(image, false);
@@ -1587,5 +1459,6 @@
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->modelPSF, false);
+    PS_ASSERT_PTR_NON_NULL(source->modelFLT, false);
+    // XXX EAM : what is the purpose of modelPSF/modelFLT?
 
     //
@@ -1610,4 +1483,5 @@
             psFree(xVec);
             psFree(yVec);
+            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
             return(NULL);
             //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
@@ -1622,8 +1496,9 @@
             psFree(xVec);
             psFree(yVec);
+            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
             return(NULL);
             //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
         }
-        //printf("The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);
+        psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);
         xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;
 
@@ -1646,4 +1521,5 @@
             psFree(xVec);
             psFree(yVec);
+            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
             return(NULL);
             //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
@@ -1657,4 +1533,5 @@
             psFree(xVec);
             psFree(yVec);
+            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
             return(NULL);
             //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
@@ -1672,37 +1549,43 @@
     tmpArray->data[0] = (psPtr *) yVec;
     tmpArray->data[1] = (psPtr *) xVec;
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     return(tmpArray);
 }
 
-#if 0
-static psVector *minLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
-static psVector *minLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
-static psVector *minLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
-static psVector *minLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x);
-static psVector *minLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x);
-static psVector *minLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x);
-#endif
-
 // XXX EAM : these are better starting values, but should be available from metadata?
-#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 20
+#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
 /******************************************************************************
-pmSourceFitModel(source, image): must create the appropiate arguments to the
+pmSourceFitModel(source, model): must create the appropiate arguments to the
 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
  
 XXX: should there be a mask value?
-XXX: Probably should remove the "image" argument.
-*****************************************************************************/
-bool pmSourceFitModel(pmSource *source,
-                      const psImage *image)
-{
+XXX EAM : fit the specified model (not necessarily the one in source)
+*****************************************************************************/
+bool pmSourceFitModel_v5(pmSource *source,
+                         pmModel *model,
+                         const bool PSF)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(source->moments, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->modelPSF, false);
-    PS_ASSERT_IMAGE_NON_NULL(image, false);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
-    psBool rc;
+    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->weight, false);
+    psBool fitStatus = true;
+    psBool onPic     = true;
+    psBool rc        = true;
+
+    // XXX EAM : is it necessary for the mask & weight to exist?  the
+    //           tests below could be conditions (!NULL)
+
+    psVector *params = model->params;
+    psVector *dparams = model->dparams;
+    psVector *paramMask = NULL;
+
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+
+    int nParams = PSF ? params->n - 4 : params->n;
 
     // find the number of valid pixels
@@ -1716,4 +1599,9 @@
             }
         }
+    }
+    if (count <  nParams + 1) {
+        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return(false);
     }
 
@@ -1733,7 +1621,7 @@
                 x->data[tmpCnt] = (psPtr *) coord;
                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
-
-                // XXX EAM : this is approximate: need to apply the gain and rdnoise
-                yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j]));
+                yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]);
+                // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
+                //           the minimization function calculates sq()
                 tmpCnt++;
             }
@@ -1744,59 +1632,241 @@
                             PM_SOURCE_FIT_MODEL_TOLERANCE);
 
-    psVector *params = source->modelPSF->params;
-
-    switch (source->modelPSF->type) {
-    case PS_MODEL_GAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D);
-        break;
-    case PS_MODEL_PGAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_PsuedoGauss2D);
-        break;
-    case PS_MODEL_TWIST_GAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D);
-        break;
-    case PS_MODEL_WAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D);
-        break;
-    case PS_MODEL_SERSIC:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Sersic);
-        break;
-    case PS_MODEL_SERSIC_CORE:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore);
-        break;
-    default:
-        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
-        rc = false;
-    }
+    // PSF model only fits first 4 parameters, FLT model fits all
+    if (PSF) {
+        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
+        for (int i = 0; i < 4; i++) {
+            paramMask->data.U8[i] = 0;
+        }
+        for (int i = 4; i < paramMask->n; i++) {
+            paramMask->data.U8[i] = 1;
+        }
+    }
+
+    // XXX EAM : covar must be F64?
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+
+    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
+    for (int i = 0; i < dparams->n; i++) {
+        if ((paramMask != NULL) && paramMask->data.U8[i])
+            continue;
+        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
+    }
+
     // 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
-    source->modelPSF->chisq = myMin->value;
-    source->modelPSF->nDOF  = y->n - params->n;
-    source->modelPSF->nIter = myMin->iter;
-
+    model->chisq = myMin->value;
+    model->nIter = myMin->iter;
+    model->nDOF  = y->n - nParams;
+
+    // XXX EAM get the Gauss-Newton distance for fixed model parameters
+    if (paramMask != NULL) {
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
+        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
+        for (int i = 0; i < dparams->n; i++) {
+            if (!paramMask->data.U8[i])
+                continue;
+            dparams->data.F32[i] = delta->data.F64[i];
+        }
+        psFree (delta);
+    }
+
+    psFree(x);
+    psFree(y);
+    psFree(yErr);
+    psFree(myMin);
+    psFree(covar);
+    psFree(paramMask);
+
+    rc = (onPic && fitStatus);
+    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
+    return(rc);
+}
+
+// XXX EAM : new version with parameter range limits and weight enhancement
+bool pmSourceFitModel (pmSource *source,
+                       pmModel *model,
+                       const bool PSF)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->moments, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->weight, false);
+
+    // XXX EAM : is it necessary for the mask & weight to exist?  the
+    //           tests below could be conditions (!NULL)
+
+    psBool fitStatus = true;
+    psBool onPic     = true;
+    psBool rc        = true;
+    psF32  Ro, ymodel;
+
+    psVector *params = model->params;
+    psVector *dparams = model->dparams;
+    psVector *paramMask = NULL;
+
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+
+    // XXX EAM : I need to use the sky value to constrain the weight model
+    int nParams = PSF ? params->n - 4 : params->n;
+    psF32 So = params->data.F32[0];
+
+    // find the number of valid pixels
+    // XXX EAM : this loop and the loop below could just be one pass
+    //           using the psArrayAdd and psVectorExtend functions
+    psS32 count = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            if (source->mask->data.U8[i][j] == 0) {
+                count++;
+            }
+        }
+    }
+    if (count <  nParams + 1) {
+        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return(false);
+    }
+
+    // construct the coordinate and value entries
+    psArray *x = psArrayAlloc(count);
+    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
+    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
+    psS32 tmpCnt = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            if (source->mask->data.U8[i][j] == 0) {
+                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+                // XXX: Convert i/j to image space:
+                // XXX EAM: coord order is (x,y) == (col,row)
+                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
+                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+                x->data[tmpCnt] = (psPtr *) coord;
+                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
+
+                // compare observed flux to model flux to adjust weight
+                ymodel = modelFunc (NULL, model->params, coord);
+
+                // this test enhances the weight based on deviation from the model flux
+                Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
+
+                // psMinimizeLMChi2_EAM takes wt = 1/dY^2
+                if (source->weight->data.F32[i][j] == 0) {
+                    yErr->data.F32[tmpCnt] = 0.0;
+                } else {
+                    yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
+                }
+                tmpCnt++;
+            }
+        }
+    }
+
+    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
+                            PM_SOURCE_FIT_MODEL_TOLERANCE);
+
+    // PSF model only fits first 4 parameters, FLT model fits all
+    if (PSF) {
+        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
+        for (int i = 0; i < 4; i++) {
+            paramMask->data.U8[i] = 0;
+        }
+        for (int i = 4; i < paramMask->n; i++) {
+            paramMask->data.U8[i] = 1;
+        }
+    }
+
+    // XXX EAM : I've added three types of parameter range checks
+    // XXX EAM : this requires my new psMinimization functions
+    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
+    psVector *beta_lim = NULL;
+    psVector *params_min = NULL;
+    psVector *params_max = NULL;
+
+    // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
+    //           in the SDRS, I define a new psMinimization which will take these in
+    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
+    modelLimits (&beta_lim, &params_min, &params_max);
+    for (int i = 0; i < params->n; i++) {
+        covar->data.F64[0][i] = beta_lim->data.F32[i];
+        covar->data.F64[1][i] = params_min->data.F32[i];
+        covar->data.F64[2][i] = params_max->data.F32[i];
+    }
+
+    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
+    for (int i = 0; i < dparams->n; i++) {
+        if ((paramMask != NULL) && paramMask->data.U8[i])
+            continue;
+        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
+    }
+
+    // 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
+    model->chisq = myMin->value;
+    model->nIter = myMin->iter;
+    model->nDOF  = y->n - nParams;
+
+    // XXX EAM get the Gauss-Newton distance for fixed model parameters
+    if (paramMask != NULL) {
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
+        psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);
+        for (int i = 0; i < dparams->n; i++) {
+            if (!paramMask->data.U8[i])
+                continue;
+            dparams->data.F32[i] = delta->data.F64[i];
+        }
+    }
+
+    psFree(paramMask);
     psFree(x);
     psFree(y);
     psFree(myMin);
+
+    rc = (onPic && fitStatus);
+    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
 }
 
-static bool sourceAddOrSubModel(psImage *image,
-                                pmSource *src,
-                                bool center,
-                                psS32 flag)
-{
-    PS_ASSERT_PTR_NON_NULL(src, false);
-    PS_ASSERT_PTR_NON_NULL(src->moments, false);
-    PS_ASSERT_PTR_NON_NULL(src->peak, false);
-    PS_ASSERT_PTR_NON_NULL(src->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(src->modelPSF, false);
+bool p_pmSourceAddOrSubModel(psImage *image,
+                             psImage *mask,
+                             pmModel *model,
+                             bool center,
+                             psS32 flag)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+
+    PS_ASSERT_PTR_NON_NULL(model, false);
     PS_ASSERT_IMAGE_NON_NULL(image, false);
     PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
 
     psVector *x = psVectorAlloc(2, PS_TYPE_F32);
-    psVector *params = src->modelPSF->params;
-
-    for (psS32 i = 0; i < src->pixels->numRows; i++) {
-        for (psS32 j = 0; j < src->pixels->numCols; j++) {
+    psVector *params = model->params;
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+    psS32 imageCol;
+    psS32 imageRow;
+
+    for (psS32 i = 0; i < image->numRows; i++) {
+        for (psS32 j = 0; j < image->numCols; j++) {
+            if ((mask != NULL) && mask->data.U8[i][j])
+                continue;
             psF32 pixelValue;
             // XXX: Should you be adding the pixels for the entire subImage,
@@ -1805,33 +1875,19 @@
             // Convert i/j to imace coord space:
             // XXX: Make sure you have col/row order correct.
-            psS32 imageRow = i + src->pixels->row0;
-            psS32 imageCol = j + src->pixels->col0;
+            // XXX EAM : 'center' option changes this
+            // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2]
+            if (center) {
+                imageCol = j - 0.5*image->numCols + model->params->data.F32[2];
+                imageRow = i - 0.5*image->numRows + model->params->data.F32[3];
+            } else {
+                imageCol = j + image->col0;
+                imageRow = i + image->row0;
+            }
 
             x->data.F32[0] = (float) imageCol;
             x->data.F32[1] = (float) imageRow;
-            switch (src->modelPSF->type) {
-            case PS_MODEL_GAUSS:
-                pixelValue = pmMinLM_Gauss2D(NULL, params, x);
-                break;
-            case PS_MODEL_PGAUSS:
-                pixelValue = pmMinLM_PsuedoGauss2D(NULL, params, x);
-                break;
-            case PS_MODEL_TWIST_GAUSS:
-                pixelValue = pmMinLM_TwistGauss2D(NULL, params, x);
-                break;
-            case PS_MODEL_WAUSS:
-                pixelValue = pmMinLM_Wauss2D(NULL, params, x);
-                break;
-            case PS_MODEL_SERSIC:
-                pixelValue = pmMinLM_Sersic(NULL, params, x);
-                break;
-            case PS_MODEL_SERSIC_CORE:
-                pixelValue = pmMinLM_SersicCore(NULL, params, x);
-                break;
-            default:
-                psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
-                psFree(x);
-                return(false);
-            }
+            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;
@@ -1841,8 +1897,9 @@
             // how to use the boolean "center" flag.
 
-            image->data.F32[imageRow][imageCol]+= pixelValue;
+            image->data.F32[i][j]+= pixelValue;
         }
     }
     psFree(x);
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     return(true);
 }
@@ -1853,8 +1910,12 @@
  *****************************************************************************/
 bool pmSourceAddModel(psImage *image,
-                      pmSource *src,
+                      psImage *mask,
+                      pmModel *model,
                       bool center)
 {
-    return(sourceAddOrSubModel(image, src, center, 0));
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0);
+    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
+    return(rc);
 }
 
@@ -1862,259 +1923,13 @@
  *****************************************************************************/
 bool pmSourceSubModel(psImage *image,
-                      pmSource *src,
+                      psImage *mask,
+                      pmModel *model,
                       bool center)
 {
-    return(sourceAddOrSubModel(image, src, center, 1));
-}
-
-
-// XXX: Put this is psConstants.h
-#define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \
-if (VEC1->n != N) { \
-    psError(PS_ERR_BAD_PARAMETER_SIZE, true, \
-            "psVector %s has size %d, should be %d.", \
-            #VEC1, VEC1->n, N); \
-    return(RVAL); \
-}
-
-
-/**
-   all of these object representation functions have the same form : func(*deriv, *params, *x)
- 
-   the argument "x" contains a single "x,y" coordinate pair.  The function computes the object
-   model, based on the parameters in "params" at the x,y point specified by *x, and returns the value.
-   The derivatives are also caculated and returned in the "deriv" argument.  parameter error checking is
-   skipped because speed is most important.
-**/
-
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = sqrt(2.0) / SigmaX;
-    params->data.F32[5] = sqrt(2.0) / SigmaY;
-    params->data.F32[6] = Sxy;
-*****************************************************************************/
-float pmMinLM_Gauss2D(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psF32 X  = x->data.F32[0] - params->data.F32[2];
-    psF32 Y  = x->data.F32[1] - params->data.F32[3];
-    psF32 px = params->data.F32[4]*X;
-    psF32 py = params->data.F32[5]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 r  = exp(-z);
-    psF32 q  = params->data.F32[1]*r;
-    psF32 f  = q + params->data.F32[0];
-
-    if (deriv != NULL) {
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y);
-        deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-    }
-    return(f);
-}
-
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = sqrt(2) / SigmaX;
-    params->data.F32[5] = sqrt(2) / SigmaY;
-    params->data.F32[6] = Sxy;
-*****************************************************************************/
-float pmMinLM_PsuedoGauss2D(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psF32 X  = x->data.F32[0] - params->data.F32[2];
-    psF32 Y  = x->data.F32[1] - params->data.F32[3];
-    psF32 px = params->data.F32[4]*X;
-    psF32 py = params->data.F32[5]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 t  = 1 + z + 0.5*z*z;
-    psF32 r  = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
-    psF32 f  = params->data.F32[1]*r + params->data.F32[0];
-
-    if (deriv != NULL) {
-        // note difference from a pure gaussian: q = params->data.F32[1]*r
-        psF32 q = params->data.F32[1]*r*r*t;
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
-        deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-    }
-    return(f);
-}
-
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = Sx;
-    params->data.F32[5] = Sy;
-    params->data.F32[6] = Sxy;
-    params->data.F32[7] = B2;
-    params->data.F32[8] = B3;
-*****************************************************************************/
-float pmMinLM_Wauss2D(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psF32 X = x->data.F32[0] - params->data.F32[2];
-    psF32 Y = x->data.F32[1] - params->data.F32[2];
-    psF32 px = params->data.F32[4]*X;
-    psF32 py = params->data.F32[5]*Y;
-    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0);
-    psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */
-    psF32 f = params->data.F32[1]*r + params->data.F32[0];
-
-    if (deriv != NULL) {
-        // note difference from gaussian: q = params->data.F32[1]*r
-        psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0));
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
-        deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-        deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t;
-        deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0;
-        // The values of 100 dampen the swing of params->data.F32[7,8] */
-    }
-    return(f);
-}
-
-// XXX: What should these be?
-#define FFACTOR 1.0
-#define FSCALE 1.0
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = SxInner;
-    params->data.F32[5] = SyInner;
-    params->data.F32[6] = SxyInner;
-    params->data.F32[7] = SxOuter;
-    params->data.F32[8] = SyOuter;
-    params->data.F32[9] = SxyOuter;
-    params->data.F32[10] = N;
-*****************************************************************************/
-float pmMinLM_TwistGauss2D(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psF32 X = x->data.F32[0] - params->data.F32[2];
-    psF32 Y = x->data.F32[1] - params->data.F32[3];
-    psF32 px1 = params->data.F32[4]*X;
-    psF32 py1 = params->data.F32[5]*Y;
-    psF32 px2 = params->data.F32[7]*X;
-    psF32 py2 = params->data.F32[8]*Y;
-    psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y;
-    psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y;
-    psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10]));
-
-    psF32 f = params->data.F32[5]*r + params->data.F32[6];
-
-    if (deriv != NULL) {
-        psF32 q1 = params->data.F32[5]*PS_SQR(r);
-        psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0));
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y);
-        deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X);
-
-        // These fudge factors impede the growth of params->data.F32[4] beyond
-        // params->data.F32[7].
-        psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]);
-        psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
-        deriv->data.F32[4] = -2.0*q1*px1*X*f2;
-
-        // These fudge factors impede the growth of params->data.F32[5] beyond
-        // params->data.F32[8].
-        f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]);
-        f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
-        deriv->data.F32[5] = -2.0*q1*py1*Y*f2;
-        deriv->data.F32[6] = -q1*X*Y;
-        deriv->data.F32[7] = -2.0*q2*px2*X;
-        deriv->data.F32[8] = -2.0*q2*py2*Y;
-        deriv->data.F32[9] = -q2*X*Y;
-        deriv->data.F32[10] = -q1*log(z2);
-    }
-
-    return(f);
-}
-
-/******************************************************************************
-    float Sersic()
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = Sx;
-    params->data.F32[5] = Sy;
-    params->data.F32[6] = Sxy;
-    params->data.F32[7] = Nexp;
-*****************************************************************************/
-float pmMinLM_Sersic(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
-    return(0.0);
-}
-
-/******************************************************************************
-    float SersicBulge()
-    params->data.F32[0] So;
-    params->data.F32[1] Zo;
-    params->data.F32[2] Xo;
-    params->data.F32[3] Yo;
-    params->data.F32[4] SxInner;
-    params->data.F32[5] SyInner;
-    params->data.F32[6] SxyInner;
-    params->data.F32[7] Zd;
-    params->data.F32[8] SxOuter;
-    params->data.F32[9] SyOuter;
-    params->data.F32[10] = SxyOuter;
-    params->data.F32[11] = Nexp;
-*****************************************************************************/
-float pmMinLM_SersicCore(
-    psVector *deriv,
-    const psVector *params,
-    const psVector *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
-    return(0.0);
-}
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1);
+    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 5254)
+++ /trunk/psModules/src/objects/pmObjects.h	(revision 5255)
@@ -1,10 +1,15 @@
 /** @file  pmObjects.h
  *
- *  This file will ...
+ * The process of finding, measuring, and classifying astronomical sources on
+ * 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 
+ * are generally low-level components which can be connected together to
+ * construct a complete object measurement suite.
  *
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-28 20:43:52 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-10-10 19:53:40 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -19,7 +24,23 @@
 #endif
 
-#include<stdio.h>
-#include<math.h>
+#include <stdio.h>
+#include <math.h>
 #include "pslib.h"
+#include "pmAstrometry.h"
+/**
+ * In the object analysis process, we will use specific mask values to mark the
+ * image pixels. The following structure defines the relevant mask values.
+ *
+ * XXX: This is probably a bad solution: we will want to set mask values
+ * outside of the PSPHOT code.  Perhaps we can set up a registered set of mask
+ * values with specific meanings that other functions can add to or define?
+ */
+enum {
+    PSPHOT_MASK_CLEAR     = 0x00,
+    PSPHOT_MASK_INVALID   = 0x01,
+    PSPHOT_MASK_SATURATED = 0x02,
+    PSPHOT_MASK_MARKED    = 0x08,
+} psphotMaskValues;
+
 
 /** pmPeakType
@@ -40,5 +61,10 @@
 } pmPeakType;
 
+
 /** 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:
  *  
  */
@@ -52,63 +78,101 @@
 pmPeak;
 
+
 /** 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
 {
-    float x;                            ///< X-coord of centroid.
-    float y;                            ///< Y-coord of centroid.
-    float Sx;                           ///< x-second moment.
-    float Sy;                           ///< y-second moment.
-    float Sxy;                          ///< xy cross moment.
-    float Sum;                          ///< Pixel sum above sky (background).
-    float Peak;                         ///< Peak counts above sky.
-    float Sky;                          ///< Sky level (background).
-    int nPixels;                        ///< Number of pixels used.
+    float x;    ///< X-coord of centroid.
+    float y;    ///< Y-coord of centroid.
+    float Sx;    ///< x-second moment.
+    float Sy;    ///< y-second moment.
+    float Sxy;    ///< xy cross moment.
+    float Sum;    ///< Pixel sum above sky (background).
+    float Peak;    ///< Peak counts above sky.
+    float Sky;    ///< Sky level (background).
+    float SN;    ///< approx signal-to-noise
+    int nPixels;   ///< Number of pixels used.
 }
 pmMoments;
 
-/** pmModelType enumeration
- *  
- */
-typedef enum {
-    PS_MODEL_GAUSS,                     ///< Regular 2-D Gaussian
-    PS_MODEL_PGAUSS,                    ///< Psuedo 2-D Gaussian
-    PS_MODEL_TWIST_GAUSS,               ///< 2-D Twisted Gaussian
-    PS_MODEL_WAUSS,                     ///< 2-D Waussian
-    PS_MODEL_SERSIC,                    ///< Sersic
-    PS_MODEL_SERSIC_CORE,               ///< Sersic Core
-    PS_MODEL_UNDEFINED                  ///< Undefined
-} pmModelType;
-
-/** pmModel data structure
- *  
- */
-// XXX: The SDRS has the "type" member of type psS32.
+
+/** 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
 {
-    pmModelType type;                   ///< Model to be used.
-    psVector *params;                   ///< Paramater values.
-    psVector *dparams;                  ///< Parameter errors.
-    float chisq;                        ///< Fit chi-squared.
-    int nDOF;                           ///< number of degrees of freedom
-    int nIter;                          ///< number of iterations to reach min
+    float X;
+    float dX;
+    float Y;
+    float dY;
+}
+pmPSFClump;
+
+typedef int pmModelType;
+#define PS_MODEL_GAUSS 0
+#define PS_MODEL_PGAUSS 1
+#define PS_MODEL_QGAUSS 2
+#define PS_MODEL_SGAUSS 3
+
+
+/** 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
+ * object. In this case, the PSF-dependent parameters are specified for the
+ * 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
+{
+    pmModelType type;   ///< Model to be used.
+    psVector *params;   ///< Paramater values.
+    psVector *dparams;   ///< Parameter errors.
+    float chisq;   ///< Fit chi-squared.
+    int nDOF;    ///< number of degrees of freedom
+    int nIter;    ///< number of iterations to reach min
+    float radius;   ///< fit radius actually used
 }
 pmModel;
 
 /** 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 {
-    PS_SOURCE_PSFSTAR,
-    PS_SOURCE_GALAXY,
-    PS_SOURCE_DEFECT,
-    PS_SOURCE_SATURATED,
-    PS_SOURCE_SATSTAR,
-    PS_SOURCE_FAINTSTAR,
-    PS_SOURCE_BRIGHTSTAR,
-    PS_SOURCE_OTHER
+    PM_SOURCE_DEFECT,                   ///< a cosmic-ray
+    PM_SOURCE_SATURATED,                ///< random saturated pixels
+
+    PM_SOURCE_SATSTAR,                  ///< a saturated star
+    PM_SOURCE_PSFSTAR,                  ///< a PSF star
+    PM_SOURCE_GOODSTAR,                 ///< a good-quality star
+
+    PM_SOURCE_POOR_FIT_PSF,             ///< poor quality PSF fit
+    PM_SOURCE_FAIL_FIT_PSF,             ///< failed to get a good PSF fit
+    PM_SOURCE_FAINTSTAR,                ///< below S/N cutoff
+
+    PM_SOURCE_GALAXY,                   ///< an extended object (galaxy)
+    PM_SOURCE_FAINT_GALAXY,             ///< a galaxy below S/N cutoff
+    PM_SOURCE_DROP_GALAXY,              ///< ?
+    PM_SOURCE_FAIL_FIT_GAL,             ///< failed on the galaxy fit
+    PM_SOURCE_POOR_FIT_GAL,             ///< poor quality galaxy fit
+
+    PM_SOURCE_OTHER,                    ///< unidentified
 } pmSourceType;
 
@@ -122,207 +186,344 @@
 typedef struct
 {
-    pmPeak *peak;                       ///< Description of peak pixel.
-    psImage *pixels;                    ///< Rectangular region including object pixels.
-    psImage *mask;                      ///< Mask which marks pixels associated with objects.
-    pmMoments *moments;                 ///< Basic moments measure for the object.
-    pmModel *modelPSF;                  ///< PSF model parameters and type
-    pmModel *modelFLT;                  ///< FLT model parameters and type
-    pmSourceType type;                  ///< Best identification of object.
+    pmPeak *peak;   ///< Description of peak pixel.
+    psImage *pixels;   ///< Rectangular region including object pixels.
+    psImage *weight;   ///< Image variance.
+    psImage *mask;   ///< Mask which marks pixels associated with objects.
+    pmMoments *moments;   ///< Basic moments measure for the object.
+    pmModel *modelPSF;   ///< PSF Model fit (parameters and type)
+    pmModel *modelFLT;   ///< FLT (floating) Model fit (parameters and type).
+    pmSourceType type;   ///< Best identification of object.
 }
 pmSource;
 
-/** pmPeak data structure
- *  
- *  
- *  
- */
-typedef struct
-{
-    psS32 type;                         ///< PSF Model in use
-    psArray *params;                    ///< Model parameters (psPolynomial2D)
-    psF32 chisq;                        ///< PSF goodness statistic
-    psS32 nPSFstars;                    ///< number of stars used to measure PSF
-}
-pmPSF;
-
-
-
+
+/** pmPeakAlloc()
+ *
+ *  @return pmPeak*    newly allocated pmPeak with all internal pointers set to NULL
+ */
 pmPeak *pmPeakAlloc(
-    int x,                              ///< Row-coordinate in image space
-    int y,                              ///< Col-coordinate in image space
-    float counts,                       ///< The value of the peak pixel
-    pmPeakType class                    ///< The type of peak pixel
-);
-
+    int x,    ///< Row-coordinate in image space
+    int y,    ///< Col-coordinate in image space
+    float counts,   ///< The value of the peak pixel
+    pmPeakType class   ///< The type of peak pixel
+);
+
+
+/** pmMomentsAlloc()
+ * 
+ */
 pmMoments *pmMomentsAlloc();
+
+
+/** pmModelAlloc()
+ * 
+ */
 pmModel *pmModelAlloc(pmModelType type);
-pmSource *pmSourceAlloc();
-
-/******************************************************************************
-pmFindVectorPeaks(vector, threshold): Find all local peaks in the given vector
-above the given threshold.  Returns a vector of type PS_TYPE_U32 containing
-the location (x value) of all peaks.
- *****************************************************************************/
+
+
+/** pmSourceAlloc()
+ * 
+ */
+pmSource  *pmSourceAlloc();
+
+
+/** 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
+ * a value above the threshold. Two types of special cases must be addressed.
+ * Equal value elements: If an element has the same value as the following
+ * element, it is not considered a peak. If an element has the same value as the
+ * preceding element (but not the following), then it is considered a peak. Note
+ * that this rule (arbitrarily) identifies flat regions by their trailing edge.
+ * Edge cases: At start of the vector, the element must be higher than its
+ * neighbor. At the end of the vector, the element must be higher or equal to its
+ * neighbor. These two rules again places the peak associated with a flat region
+ * which touches the image edge at the image edge. The result of this function is
+ * a vector containing the coordinates (element number) of the detected peaks
+ * (type psU32).
+ * 
+ */
 psVector *pmFindVectorPeaks(
-    const psVector *vector,             ///< The input vector (float)
-    float threshold                     ///< Threshold above which to find a peak
-);
-
-/******************************************************************************
-pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
-above the given threshold.  Returns a psList containing the location (x/y
-value) of all peaks.
- *****************************************************************************/
+    const psVector *vector,  ///< The input vector (float)
+    float threshold   ///< Threshold above which to find a peak
+);
+
+
+/** 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
+ * peak and exclude peaks which are not local peaks. A peak is a local peak if it
+ * has a higher value than all 8 neighbors. If the peak has the same value as its
+ * +y neighbor or +x neighbor, it is NOT a local peak. If any other neighbors
+ * have an equal value, the peak is considered a valid peak. Note two points:
+ * first, the +x neighbor condition is already enforced by pmFindVectorPeaks.
+ * Second, these rules have the effect of making flat-topped regions have single
+ * 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(
-    const psImage *image,               ///< The input image where peaks will be found (float)
-    float threshold                     ///< Threshold above which to find a peak
-);
-
-/******************************************************************************
-psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psList that have
-a peak value above the given maximum, or fall outside the valid region.
- *****************************************************************************/
+    const psImage *image,  ///< The input image where peaks will be found (float)
+    float threshold   ///< Threshold above which to find a peak
+);
+
+
+/** pmCullPeaks()
+ * 
+ * Eliminate peaks from the psList that have a peak value above the given
+ * maximum, or fall outside the valid region.
+ * 
+ */
 psList *pmCullPeaks(
-    psList *peaks,                      ///< The psList of peaks to be culled
-    float maxValue,                     ///< Cull peaks above this value
-    const psRegion valid               ///< Cull peaks otside this psRegion
-);
-
-/******************************************************************************
-pmSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius):
- 
- *****************************************************************************/
-pmSource *pmSourceLocalSky(
-    const psImage *image,               ///< The input image (float)
-    const pmPeak *peak,                 ///< The peak for which the psSource struct is created.
-    psStatsOptions statsOptions,        ///< The statistic used in calculating the background sky
-    float innerRadius,                  ///< The inner radius of the suqare annulus for calculating sky
-    float outerRadius                   ///< The outer radius of the suqare annulus for calculating sky
-);
-
-/******************************************************************************
- *****************************************************************************/
-pmSource *pmSourceMoments(
-    pmSource *source,                   ///< The input pmSource for which moments will be computed
-    float radius                        ///< Use a circle of pixels around the peak
-);
-
-/******************************************************************************
-pmSourceRoughClass(pmArray *source, psMetaDeta *metadata): make a guess at the
-source classification.
- *****************************************************************************/
+    psList *peaks,   ///< The psList of peaks to be culled
+    float maxValue,   ///< Cull peaks above this value
+    const psRegion valid                ///< Cull peaks otside this psRegion
+);
+
+
+/** 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
+ * have a peak value above the given maximum value limit or if the fall outside
+ * the valid region.  The result of the function is a new array with a reduced
+ * number of peaks.
+ * 
+ */
+psArray *pmPeaksSubset(
+    psArray *peaks,                     ///< Add comment.
+    float maxvalue,                     ///< Add comment.
+    const psRegion valid                ///< Add comment.
+);
+
+
+/** 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
+ * a square window (of full width 2Radius+1) centered on the pixel which contains
+ * the given coordinate, in the frame of the readout. The window is defined to
+ * have limits which are valid within the boundary of the readout image, thus if
+ * the radius would fall outside the image pixels, the subimage is truncated to
+ * only consist of valid pixels. If readout->mask or readout->weight are not
+ * NULL, matching subimages are defined for those images as well. This function
+ * fails if no valid pixels can be defined (x or y less than Radius, for
+ * 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.
+bool pmSourceDefinePixels(
+    pmSource *mySource,                 ///< Add comment.
+    pmReadout *readout,                 ///< Add comment.
+    psF32 x,                            ///< Add comment.
+    psF32 y,                            ///< Add comment.
+    psF32 Radius                        ///< Add comment.
+);
+
+
+/** 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
+ * function assumes the source pixels have been defined, and that the value of
+ * Radius here is smaller than the value of Radius used to define the pixels. The
+ * annular region not contained within the radius defined here is used to measure
+ * the local background in the vicinity of the source. The local background
+ * measurement uses the specified statistic passed in via the statsOptions entry.
+ * 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(
+    pmSource *source,   ///< The input image (float)
+    psStatsOptions statsOptions, ///< The statistic used in calculating the background sky
+    float Radius   ///< The inner radius of the square annulus to exclude
+);
+
+
+/** pmSourceMoments()
+ * 
+ * Measure source moments for the given source, using the value of
+ * source.moments.sky provided as the local background value and the peak
+ * coordinates as the initial source location. The resulting moment values are
+ * applied to the source.moments entry, and the source is returned. The moments
+ * are measured within the given circular radius of the source.peak coordinates.
+ * The return value indicates the success (TRUE) of the operation.
+ * 
+ */
+bool pmSourceMoments(
+    pmSource *source,   ///< The input pmSource for which moments will be computed
+    float radius   ///< Use a circle of pixels around the peak
+);
+
+
+/** 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
+ * the image. As long as the PSF shape does not vary excessively across the
+ * image, the sources which are represented by a PSF (the start) will have very
+ * similar second moments. The function pmSourcePSFClump searches a collection of
+ * sources with measured moments for a group with moments which are all very
+ * 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
+ * represent this plane.  The units of sigma_x and sigma_y are in image pixels. A
+ * pixel in this analysis image represents 0.1 pixels in the input image. The
+ * dimensions of the image need only be 10 pixels. The peak pixel in this image
+ * (above a threshold of half of the image maximum) is found. The coordinates of
+ * this peak pixel represent the 2D mode of the sigma_x, sigma_y distribution.
+ * The sources with sigma_x, sigma_y within 0.2 pixels of this value are then
+ *  * 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(
+    psArray *source,   ///< The input pmSource
+    psMetadata *metadata  ///< Contains classification parameters
+);
+
+
+/** pmSourceRoughClass()
+ * 
+ * Based on the specified data values, make a guess at the source
+ * classification. The sources are provides as a psArray of pmSource entries.
+ * Definable parameters needed to make the classification are provided to the
+ * routine with the psMetadata structure. The rules (in SDRS) refer to values which
+ * can be extracted from the metadata using the given keywords. Except as noted,
+ * the data type for these parameters are psF32.
+ * 
+ */
 bool pmSourceRoughClass(
-    psArray *source,                    ///< The input pmSource
-    psMetadata *metadata                ///< Contains classification parameters
-);
-/******************************************************************************
-pmSourceSetPixelCircle(source, image, radius)
- *****************************************************************************/
-bool pmSourceSetPixelsCircle(
-    pmSource *source,                   ///< The input pmSource
-    const psImage *image,               ///< The input image (float)
-    float radius                        ///< The radius of the circle
-);
-
-/******************************************************************************
- *****************************************************************************/
-bool pmSourceModelGuess(
-    pmSource *source,                   ///< The input pmSource
-    const psImage *image,               ///< The input image (float)
-    pmModelType model                   ///< The type of model to be created.
-);
-
-/******************************************************************************
- *****************************************************************************/
+    psArray *source,   ///< The input pmSource
+    psMetadata *metadata,  ///< Contains classification parameters
+    pmPSFClump clump   ///< Statistics about the PSF clump
+);
+
+
+/** 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
+ * provided model selection. The method of defining the model parameter guesses
+ * 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(
+    pmSource *source,   ///< The input pmSource
+    pmModelType model   ///< The type of model to be created.
+);
+
+
+/** pmContourType
+ * 
+ * Only one type is defined at present.
+ * 
+ */
 typedef enum {
     PS_CONTOUR_CRUDE,
+    PS_CONTOUR_UNKNOWN01,
+    PS_CONTOUR_UNKNOWN02
 } pmContourType;
 
+
+/** 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,
+ * running along each pixel row until the level is crossed, then interpolating to
+ * the level coordinate for that row. This is done for each row, with the
+ * starting point determined by the midpoint of the previous row, until the
+ * starting point has a value below the contour level. The returned contour
+ * consists of two vectors giving the x and y coordinates of the contour levels.
+ * 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(
-    pmSource *source,                   ///< The input pmSource
-    const psImage *image,               ///< The input image (float) (this arg should be removed)
-    float level,                        ///< The level of the contour
-    pmContourType mode                  ///< Currently this must be PS_CONTOUR_CRUDE
-);
-
-/******************************************************************************
- *****************************************************************************/
+    pmSource *source,   ///< The input pmSource
+    const psImage *image,  ///< The input image (float) (this arg should be removed)
+    float level,   ///< The level of the contour
+    pmContourType mode   ///< Currently this must be PS_CONTOUR_CRUDE
+);
+
+
+/** 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
+ * interest are specified by the source.pixelsand source.maskentries. This
+ * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
+ * on success or FALSE on failure.
+ * 
+ */
 bool pmSourceFitModel(
-    pmSource *source,                   ///< The input pmSource
-    const psImage *image                ///< The input image (float)
-);
-
-/******************************************************************************
- *****************************************************************************/
+    pmSource *source,   ///< The input pmSource
+    pmModel *model,   ///< model to be fitted
+    const bool PSF   ///< Treat model as PSF or FLT?
+);
+
+
+/** 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(
+    pmModel *model                      ///< Add comment.
+);
+
+
+/** 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
+ * it is placed at its centroid location. The boolean option sky selects if the
+ * background sky is applied (TRUE) or not. The pixel range in the target image
+ * is at most the pixel range specified by the source.pixels image. The success
+ * status is returned.
+ * 
+ */
 bool pmSourceAddModel(
-    psImage *image,                     ///< The opuut image (float)
-    pmSource *source,                   ///< The input pmSource
-    bool center                         ///< A boolean flag that determines whether pixels are centered
-);
-
-/******************************************************************************
- *****************************************************************************/
+    psImage *image,   ///< The output image (float)
+    psImage *mask,   ///< The image pixel mask (valid == 0)
+    pmModel *model,   ///< The input pmModel
+    bool center    ///< A boolean flag that determines whether pixels are centered
+);
+
+
+/** 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
+ * or if it is placed at its centroid location. The boolean option sky selects if
+ * the background sky is applied (TRUE) or not. The pixel range in the target
+ * image is at most the pixel range specified by the source.pixels image. The
+ * success status is returned.
+ * 
+ */
 bool pmSourceSubModel(
-    psImage *image,                     ///< The output image (float)
-    pmSource *source,                   ///< The input pmSource
-    bool center                         ///< A boolean flag that determines whether pixels are centered
-);
-
-/******************************************************************************
-XXX: Why only *x argument?
-XXX EAM: psMinimizeLMChi2Func returns psF64, not float
- *****************************************************************************/
-float pmMinLM_Gauss2D(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_PsuedoGauss2D(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_Wauss2D(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_TwistGauss2D(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_Sersic(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_SersicCore(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
-);
-
-/******************************************************************************
- *****************************************************************************/
-float pmMinLM_PsuedoSersic(
-    psVector *deriv,                    ///< A possibly-NULL structure for the output derivatives
-    const psVector *params,             ///< A psVector which holds the parameters of this function
-    const psVector *x                   ///< A psVector which holds the row/col coordinate
+    psImage *image,   ///< The output image (float)
+    psImage *mask,   ///< The image pixel mask (valid == 0)
+    pmModel *model,   ///< The input pmModel
+    bool center    ///< A boolean flag that determines whether pixels are centered
 );
 
@@ -330,68 +531,92 @@
 /**
  * 
- *  The object model functions are defined to allow for the flexible addition
- *  of new object models. Every object model, with parameters represented by
- *  pmModel, has an associated set of functions which provide necessary support
- *  operations. A set of abstract functions allow the programmer to select the
- *  approriate function or property for a specific named object model.
- * 
- */
+ * 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
+ * to infinity. In practice, we integrate the model beyond 50sigma.  The aperture magnitude is
+ * defined as -2.5log(flux) , where the flux is summed for all pixels which are
+ * 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(
+    float *fitMag,                      ///< integrated fit magnitude
+    float *obsMag,   ///< aperture flux magnitude
+    pmModel *model,                     ///< model used for photometry
+    psImage *image,                     ///< image pixels to be used
+    psImage *mask                       ///< mask of pixels to ignore
+);
+
 
 /**
  * 
- *  This function is the model chi-square minimization function for this model.
- * 
- */
-typedef psMinimizeLMChi2Func pmModelFunc;
-
-
-/**
- * 
- * This function returns the integrated flux for the given model parameters.
- */
-typedef psF64 (*pmModelFlux)(const psVector *params);
-
-
-/**
- * 
- *  This function provides the model guess parameters based on the details of
- *   the given source.
- * 
- */
-typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
-
-
-/**
- * 
- *  This function constructs the PSF model for the given source based on the
- *  supplied psf and the FLT model for the object.
- * 
- */
-typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
-
-
-/**
- * 
- *  This function returns the radius at which the given model and parameters
- *  achieves the given flux.
- * 
- */
-typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
-
-
-/**
- * 
- *  Each of the function types above has a corresponding function which returns
- *  the function given the model type:
- * 
- */
-pmModelFunc pmModelFunc_GetFunction (pmModelType type);
-pmModelFlux pmModelFlux_GetFunction (pmModelType type);
-pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type);
-pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type);
-pmModelRadius pmModelRadius_GetFunction (pmModelType type);
-psS32 pmModelParameterCount(pmModelType type);
-psS32 pmModelSetType(char *name);
-char *pmModelGetType(pmModelType type);
+ * This function converts the source classification into the closest available
+ * approximation to the Dophot classification scheme:
+ * 
+ * PM_SOURCE_DEFECT: 8
+ * PM_SOURCE_SATURATED: 8
+ * PM_SOURCE_SATSTAR: 10
+ * PM_SOURCE_PSFSTAR: 1
+ * PM_SOURCE_GOODSTAR: 1
+ * PM_SOURCE_POOR_FIT_PSF: 7
+ * PM_SOURCE_FAIL_FIT_PSF: 4
+ * PM_SOURCE_FAINTSTAR: 4
+ * PM_SOURCE_GALAXY: 2
+ * PM_SOURCE_FAINT_GALAXY: 2
+ * PM_SOURCE_DROP_GALAXY: 2
+ * PM_SOURCE_FAIL_FIT_GAL: 2
+ * PM_SOURCE_POOR_FIT_GAL: 2
+ * PM_SOURCE_OTHER: ?
+ * 
+ */
+int pmSourceDophotType(
+    pmSource *source                    ///< Add comment.
+);
+
+
+/** 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(
+    pmSource *source                    ///< Add comment.
+);
+
+/** 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
+ * interest are specified by the source.pixelsand source.maskentries. This
+ * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
+ * on success or FALSE on failure.
+ * 
+ */
+bool pmSourceFitModel_v5(
+    pmSource *source,   ///< The input pmSource
+    pmModel *model,   ///< model to be fitted
+    const bool PSF   ///< Treat model as PSF or FLT?
+);
+
+
+/** 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
+ * interest are specified by the source.pixelsand source.maskentries. This
+ * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
+ * on success or FALSE on failure.
+ * 
+ */
+bool pmSourceFitModel_v7(
+    pmSource *source,   ///< The input pmSource
+    pmModel *model,   ///< model to be fitted
+    const bool PSF   ///< Treat model as PSF or FLT?
+);
 
 #endif
Index: /trunk/psModules/src/objects/pmPSF.c
===================================================================
--- /trunk/psModules/src/objects/pmPSF.c	(revision 5255)
+++ /trunk/psModules/src/objects/pmPSF.c	(revision 5255)
@@ -0,0 +1,191 @@
+/** @file  pmPSF.c
+ *
+ * This file contains typedefs for the Point-Spread Function and prototypes
+ * for functions that calculate the PSF.
+ *
+ *  @author EAM, IfA
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-10-10 19:53:40 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+
+/*****************************************************************************/
+/* INCLUDE FILES                                                             */
+/*****************************************************************************/
+
+#include <pslib.h>
+#include "psLibUtils.h"
+#include "pmObjects.h"
+#include "pmPSF.h"
+#include "pmModelGroup.h"
+
+/*****************************************************************************/
+/* DEFINE STATEMENTS                                                         */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* TYPE DEFINITIONS                                                          */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* GLOBAL VARIABLES                                                          */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* FILE STATIC VARIABLES                                                     */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* FUNCTION IMPLEMENTATION - LOCAL                                           */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
+/*****************************************************************************/
+
+
+
+/*****************************************************************************
+pmPSFFree(psf): function to free a pmPSF structure
+ *****************************************************************************/
+static void pmPSFFree (pmPSF *psf)
+{
+
+    if (psf == NULL)
+        return;
+
+    psFree (psf->params);
+    return;
+}
+
+
+
+/*****************************************************************************
+pmPSFAlloc (type): allocate a pmPSF.
+    NOTE: a PSF always has 4 parameters fewer than the equivalent model.
+      This is because those 4 parameters are
+ X-center
+ Y-center
+ ???: Sky background value?
+ ???: Zo?
+ *****************************************************************************/
+pmPSF *pmPSFAlloc (pmModelType type)
+{
+    int Nparams;
+
+    pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF));
+
+    psf->type     = type;
+    psf->chisq    = 0.0;
+    psf->ApResid  = 0.0;
+    psf->dApResid = 0.0;
+    psf->skyBias  = 0.0;
+
+    Nparams = pmModelParameterCount (type);
+    if (!Nparams) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return(NULL);
+    }
+
+    psf->params = psArrayAlloc(Nparams - 4);
+    for (int i = 0; i < psf->params->n; i++) {
+        // XXX EAM : make this a user-defined value?
+        // XXX EAM : future version (0.7.0?) psf->params->data[i] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1);
+        psf->params->data[i] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1);
+    }
+
+    psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
+    return(psf);
+}
+
+
+
+/*****************************************************************************
+pmPSFFromModels (*psf, *models, *mask): build a PSF from a collection of
+free-fitted models.  The PSF ignores the first 4 (independent) model
+parameters and constructs a polynomial fit to the remaining as a function of
+image coordinate.
+    input: an array of pmModels, pre-allocated psf
+Note: some of the array entries may be NULL (failed fits); ignore them.
+ *****************************************************************************/
+bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask)
+{
+
+    // construct the fit vectors from the collection of objects
+    // use the mask to ignore missing fits
+    psVector *x  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *y  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *z  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
+
+    for (int i = 0; i < models->n; i++) {
+        pmModel *model = models->data[i];
+        if (model == NULL)
+            continue;
+
+        // XXX EAM : this is fragile: x and y must be stored consistently at 2,3
+        // note that the data is provided in the F64 array
+        x->data.F64[i] = model->params->data.F32[2];
+        y->data.F64[i] = model->params->data.F32[3];
+    }
+
+    // we are doing a robust fit.  after each pass, we drop points which are
+    // more deviant than three sigma.
+    // mask is currently updated for each pass. this is inconsistent?
+
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+
+    for (int i = 0; i < psf->params->n; i++) {
+        for (int j = 0; j < models->n; j++) {
+            pmModel *model = models->data[j];
+            if (model == NULL)
+                continue;
+            z->data.F64[j] = model->params->data.F32[i + 4];
+            dz->data.F64[j] = 1;
+            // XXX EAM : need to use actual errors?
+            // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
+        }
+
+        // XXX EAM : this is the expected API (cycle 7? cycle 8?)
+        psf->params->data[i] = psVectorClipFitPolynomial2D(psf->params->data[i], stats, mask, 0xff, z, dz, x, y);
+
+        // XXX EAM : drop this when above is implemented...
+        // psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
+
+        // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
+        // psPolynomial2DDump (psf->params->data[i]);
+    }
+
+    psFree (stats);
+    psFree (x);
+    psFree (y);
+    psFree (z);
+    psFree (dz);
+    return (true);
+}
+
+
+
+/*****************************************************************************
+pmModelFromPSF (*modelFLT, *psf):  use the model position parameters to
+construct a realization of the PSF model at the object coordinates
+ *****************************************************************************/
+pmModel *pmModelFromPSF (pmModel *modelFLT, pmPSF *psf)
+{
+
+    // need to define the relationship between the modelFLT and modelPSF ?
+
+    // find function used to set the model parameters
+    pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type);
+
+    // allocate a new pmModel to hold the PSF version
+    pmModel *modelPSF = pmModelAlloc (psf->type);
+
+    // set model parameters for this source based on PSF information
+    modelFromPSFFunc (modelPSF, modelFLT, psf);
+
+    return (modelPSF);
+}
Index: /trunk/psModules/src/objects/pmPSF.h
===================================================================
--- /trunk/psModules/src/objects/pmPSF.h	(revision 5255)
+++ /trunk/psModules/src/objects/pmPSF.h	(revision 5255)
@@ -0,0 +1,87 @@
+/** @file  pmPSF.h
+ *
+ * This file contains typedefs for the Point-Spread Function and prototypes
+ * for functions that calculate the PSF.
+ *
+ *  @author EAM, IfA
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-10-10 19:53:40 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+
+# ifndef PM_PSF_H
+# define PM_PSF_H
+
+
+/** pmPSF data structure
+ * 
+ * It is useful to generate a model to define the point-spread-function which
+ * describes the flux distribution for unresolved sources in an image. In
+ * general, the PSF varies with position in the image. We allow any of the source
+ * models defined for the pmModel to represent the PSF. For a given source model,
+ * the 2D spatial variation of all of the source parameters, except the first
+ * four PSF-independent parameters, are represented as polynomial, stored in a
+ * psArray. The other elements of the structure define the quality of the PSF
+ * determination.
+ * 
+ */
+typedef struct
+{
+    pmModelType type;   ///< PSF Model in use
+    psArray *params;   ///< Model parameters (psPolynomial2D)
+    float chisq;   ///< PSF goodness statistic
+    float ApResid;                      ///< ???
+    float dApResid;                     ///< ???
+    float skyBias;                      ///< ???
+    int nPSFstars;   ///< number of stars used to measure PSF
+}
+pmPSF;
+
+
+/**
+ * 
+ * Allocator for the pmPSF structure.
+ * 
+ */
+pmPSF *pmPSFAlloc(
+    pmModelType type                    ///< Add comment
+);
+
+
+/**
+ * 
+ * This function takes a collection of pmModel fitted models from across a
+ * single image and builds a pmPSF representation of the PSF. The input array of
+ * model fits may consist of entries to be ignored (noted by a non-zero mask
+ * entry). The analysis of the models fits a 2D polynomial for each parameter to
+ * the collection of model parameters as a function of position (and
+ * normalization?). In this process, some of the input models may be marked as
+ * outliers and excluded from the fit. These elements will be marked with a
+ * specific mask value (1 == PSFTRY_MASK_OUTLIER).
+ * 
+ */
+bool pmPSFFromModels(
+    pmPSF *psf,                         ///< Add comment
+    psArray *models,                    ///< Add comment
+    psVector *mask                      ///< Add comment
+);
+
+
+/**
+ * 
+ * This function constructs a pmModel instance based on the pmPSF description
+ * of the PSF. The input is a pmModel with at least the values of the centroid
+ * coordinates (possibly normalization if this is needed) defined. The values of
+ * the PSF-dependent parameters are specified for the specific realization based
+ * on the coordinates of the object.
+ * 
+ */
+pmModel *pmModelFromPSF(
+    pmModel *model,                     ///< Add comment
+    pmPSF *psf                          ///< Add comment
+);
+
+# endif
Index: /trunk/psModules/src/objects/pmPSFtry.c
===================================================================
--- /trunk/psModules/src/objects/pmPSFtry.c	(revision 5255)
+++ /trunk/psModules/src/objects/pmPSFtry.c	(revision 5255)
@@ -0,0 +1,338 @@
+# include <pslib.h>
+# include "psLibUtils.h"
+# include "pmObjects.h"
+# include "psModulesUtils.h"
+# include "pmPSF.h"
+# include "pmPSFtry.h"
+# include "pmModelGroup.h"
+
+// ********  pmPSFtry functions  **************************************************
+// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
+// * fitted version of the model, the pmPSF fit to the fitted model parameters,
+// * and the PSF fits to the source. It also includes the statistics from the
+// * fits, both the individual sources, and the collection
+
+// free a pmPSFtry structure
+static void pmPSFtryFree (pmPSFtry *test)
+{
+
+    if (test == NULL)
+        return;
+
+    psFree (test->psf);
+    psFree (test->sources);
+    psFree (test->modelFLT);
+    psFree (test->modelPSF);
+    psFree (test->metric);
+    psFree (test->fitMag);
+    psFree (test->mask);
+    return;
+}
+
+// allocate a pmPSFtry based on the desired sources and the model (identified by name)
+pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName)
+{
+
+    pmModelType type;
+
+    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
+
+    // XXX probably need to increment ref counter
+    type           = pmModelSetType (modelName);
+    test->psf      = pmPSFAlloc (type);
+    test->sources  = psMemCopy(sources);
+    test->modelFLT = psArrayAlloc (sources->n);
+    test->modelPSF = psArrayAlloc (sources->n);
+    test->metric   = psVectorAlloc (sources->n, PS_TYPE_F64);
+    test->fitMag   = psVectorAlloc (sources->n, PS_TYPE_F64);
+    test->mask     = psVectorAlloc (sources->n, PS_TYPE_U8);
+
+    for (int i = 0; i < test->modelFLT->n; i++) {
+        test->mask->data.U8[i]  = 0;
+        test->modelFLT->data[i] = NULL;
+        test->modelPSF->data[i] = NULL;
+        test->metric->data.F64[i] = 0;
+        test->fitMag->data.F64[i] = 0;
+    }
+
+    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
+    return (test);
+}
+
+// build a pmPSFtry for the given model:
+// - fit each source with the free-floating model
+// - construct the pmPSF from the collection of models
+// - fit each source with the PSF-parameter models
+// - measure the pmPSF quality metric (dApResid)
+
+// sources used in for pmPSFtry may be masked by the analysis
+// mask values indicate the reason the source was rejected:
+
+pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS)
+{
+    bool status;
+    float obsMag;
+    float fitMag;
+    float x;
+    float y;
+    int Nflt = 0;
+    int Npsf = 0;
+
+    pmPSFtry *try
+    = pmPSFtryAlloc (sources, modelName);
+
+    // stage 1:  fit an independent model (freeModel) to all sources
+    psTimerStart ("fit");
+    for (int i = 0; i < try
+                ->sources->n; i++) {
+
+            pmSource *source = try
+                                   ->sources->data[i];
+            pmModel  *model  = pmSourceModelGuess (source, try
+                                                       ->psf->type);
+            x = source->peak->x;
+            y = source->peak->y;
+
+            // set temporary object mask and fit object
+            // fit model as FLT, not PSF
+            psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
+            status = pmSourceFitModel (source, model, false);
+            psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
+
+            // exclude the poor fits
+            if (!status) {
+                try
+                    ->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
+                psFree (model);
+                continue;
+            }
+            try
+                ->modelFLT->data[i] = model;
+            Nflt ++;
+        }
+    psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
+
+    // make this optional?
+    // DumpModelFits (try->modelFLT, "modelsFLT.dat");
+
+    // stage 2: construct a psf (pmPSF) from this collection of model fits
+    pmPSFFromModels (try
+                     ->psf, try
+                         ->modelFLT, try
+                             ->mask);
+
+    // stage 3: refit with fixed shape parameters
+    psTimerStart ("fit");
+    for (int i = 0; i < try
+                ->sources->n; i++) {
+            // masked for: bad model fit, outlier in parameters
+            if (try
+                    ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
+
+            pmSource *source = try
+                                   ->sources->data[i];
+            pmModel  *modelFLT = try
+                                     ->modelFLT->data[i];
+
+            // set shape for this model based on PSF
+            pmModel *modelPSF = pmModelFromPSF (modelFLT, try
+                                                    ->psf);
+            x = source->peak->x;
+            y = source->peak->y;
+
+            psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
+            status = pmSourceFitModel (source, modelPSF, true);
+
+            // skip poor fits
+            if (!status) {
+                try
+                    ->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL;
+                psFree (modelPSF);
+                goto next_source;
+            }
+
+            // otherwise, save the resulting model
+            try
+                ->modelPSF->data[i] = modelPSF;
+
+            // XXX : use a different aperture radius from the fit radius?
+            // XXX : use a different estimator for the local sky?
+            // XXX : pass 'source' as input?
+            if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
+                try
+                    ->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
+                goto next_source;
+            }
+
+            try
+                ->metric->data.F64[i] = obsMag - fitMag;
+            try
+                ->fitMag->data.F64[i] = fitMag;
+            Npsf ++;
+
+next_source:
+            psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
+
+        }
+    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
+
+    // make this optional
+    // DumpModelFits (try->modelPSF, "modelsPSF.dat");
+
+    // XXX this function wants aperture radius for pmSourcePhotometry
+    pmPSFtryMetric (try
+                    , RADIUS);
+    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
+              modelName, try
+                  ->psf->ApResid, try
+                      ->psf->dApResid, try
+                          ->psf->skyBias);
+
+    return (try
+           );
+}
+
+
+bool pmPSFtryMetric (pmPSFtry *try
+                     , float RADIUS)
+{
+
+    float dBin;
+    int   nKeep, nSkip;
+
+    // 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]);
+        }
+
+    // find min and max of (1/flux):
+    psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+    psVectorStats (stats, rflux, NULL, try
+                       ->mask, PSFTRY_MASK_ALL);
+
+    // build binned versions of rflux, metric
+    dBin = (stats->max - stats->min) / 10.0;
+    psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64);
+    psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
+    psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
+    psFree (stats);
+
+    psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
+
+    // group data in daBin bins, measure lower 50% mean
+    for (int i = 0; i < daBin->n; i++) {
+
+        psVector *tmp = psVectorAlloc (try
+                                       ->sources->n, PS_TYPE_F64);
+        tmp->n = 0;
+
+        // accumulate data within bin range
+        for (int j = 0; j < try
+                    ->sources->n; j++) {
+                // masked for: bad model fit, outlier in parameters
+                if (try
+                        ->mask->data.U8[j] & PSFTRY_MASK_ALL) continue;
+
+                // skip points with extreme dA values
+                if (fabs(try
+                         ->metric->data.F64[j]) > 0.5) continue;
+
+                // skip points outside of this bin
+                if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
+                    continue;
+                if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
+                    continue;
+
+                tmp->data.F64[tmp->n] = try
+                                            ->metric->data.F64[j];
+                tmp->n ++;
+            }
+
+        // is this a valid point?
+        maskB->data.U8[i] = 0;
+        if (tmp->n < 2) {
+            maskB->data.U8[i] = 1;
+            psFree (tmp);
+            continue;
+        }
+
+        // dA values are contaminated with low outliers
+        // measure statistics only on upper 50% of points
+        // this would be easier if we could sort in reverse:
+        //
+        // psVectorSort (tmp, tmp);
+        // tmp->n = 0.5*tmp->n;
+        // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+        // psVectorStats (stats, tmp, NULL, NULL, 0);
+        // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
+
+        psVectorSort (tmp, tmp);
+        nKeep = 0.5*tmp->n;
+        nSkip = tmp->n - nKeep;
+
+        psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
+        for (int j = 0; j < tmp2->n; j++) {
+            tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
+        }
+
+        stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+        psVectorStats (stats, tmp2, NULL, NULL, 0);
+        psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
+
+        daBin->data.F64[i] = stats->sampleMedian;
+
+        psFree (stats);
+        psFree (tmp);
+        psFree (tmp2);
+    }
+
+    // linear fit to rfBin, daBin
+    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
+
+    // XXX EAM : this is the intended API (cycle 7? cycle 8?)
+    poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
+
+    // XXX EAM : replace this when the above version is implemented
+    // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL);
+
+    psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
+    psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
+
+    stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
+    stats = psVectorStats (stats, daResid, NULL, maskB, 1);
+
+    try
+        ->psf->ApResid = poly->coeff[0];
+    try
+        ->psf->dApResid = stats->clippedStdev;
+    try
+        ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+    psFree (rflux);
+    psFree (rfBin);
+    psFree (daBin);
+    psFree (maskB);
+    psFree (daBinFit);
+    psFree (daResid);
+    psFree (poly);
+    psFree (stats);
+
+    return true;
+}
Index: /trunk/psModules/src/objects/pmPSFtry.h
===================================================================
--- /trunk/psModules/src/objects/pmPSFtry.h	(revision 5255)
+++ /trunk/psModules/src/objects/pmPSFtry.h	(revision 5255)
@@ -0,0 +1,112 @@
+/** @file  pmPSFtry.h
+ *
+ * This file contains code that allows the user to try to fit several
+ * PSF models to an image.
+ *
+ *  @author EAM, IfA
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-10-10 19:53:40 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+
+# ifndef PM_PSF_TRY_H
+# define PM_PSF_TRY_H
+
+
+/**
+ * 
+ * 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
+ * being tests, and contains an element to store the resulting psf
+ * representation. In addition, this structure carries the complete collection of
+ * FLT (floating parameter) and PSF (fixed parameter) model fits to each of the
+ * sources modelFLT and modelPSF. It also contains a mask which is set by the
+ * model fitting and psf fitting steps. For each model, the value of the quality
+ * metric is stored in the vector metric and the fitted instrumental magnitude is
+ * stored in fitMag. The quality metric for the PSF model is the aperture
+ * magnitude minus the fitted magnitude for each source. This collection of
+ * aperture residuals is examined in the analysis process, and a linear trend of
+ * the residual with the inverse object flux (ie, 100:4his 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 modmag) is fitted. The result of this
+ * fit is a measured sky bias (systematic error in the sky measured by the fits),
+ * an effective infinite-magnitude aperture correction (ApResid), and the scatter
+ * of the aperture correction for the ensemble of PSF stars (dApResid). The
+ * 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
+{
+    pmPSF      *psf;                    ///< Add comment.
+    psArray    *sources;                ///< pointers to the original sources
+    psArray    *modelFLT;               ///< model fits, floating parameters
+    psArray    *modelPSF;               ///< model fits, PSF parameters
+    psVector   *mask;                   ///< Add comment.
+    psVector   *metric;                 ///< Add comment.
+    psVector   *fitMag;                 ///< Add comment.
+}
+pmPSFtry;
+
+
+/** 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 {
+    PSFTRY_MASK_CLEAR    = 0x00,        ///< Add comment.
+    PSFTRY_MASK_OUTLIER  = 0x01,        ///< 1: outlier in psf polynomial fit (provided by psPolynomials)
+    PSFTRY_MASK_FLT_FAIL = 0x02,        ///< 2: flt model failed to converge
+    PSFTRY_MASK_PSF_FAIL = 0x04,        ///< 3: psf model failed to converge
+    PSFTRY_MASK_BAD_PHOT = 0x08,        ///< 4: invalid source photometry
+    PSFTRY_MASK_ALL      = 0x0f,        ///< Add comment.
+} pmPSFtryMaskValues;
+
+
+/** pmPSFtryAlloc()
+ * 
+ * Allocate a pmPSFtry data structure.
+ * 
+ */
+pmPSFtry *pmPSFtryAlloc(
+    psArray *stars,                     ///< Add comment.
+    char *modelName                     ///< Add comment.
+);
+
+
+/** 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(
+    psArray *sources,                   ///< Add comment.
+    char *modelName,                    ///< Add comment.
+    float radius                        ///< Add comment.
+);
+
+
+/** pmPSFtryMetric()
+ * 
+ * This function is used to measure the PSF model metric for the set of
+ * results contained in the pmPSFtry structure.
+ * 
+ */
+bool pmPSFtryMetric(
+    pmPSFtry *try
+    ,                      ///< Add comment.
+    float RADIUS                        ///< Add comment.
+);
+
+# endif
Index: /trunk/psModules/src/objects/psEllipse.c
===================================================================
--- /trunk/psModules/src/objects/psEllipse.c	(revision 5255)
+++ /trunk/psModules/src/objects/psEllipse.c	(revision 5255)
@@ -0,0 +1,62 @@
+# include "pslib.h"
+# include "psEllipse.h"
+
+EllipseAxes EllipseMomentsToAxes (EllipseMoments moments)
+{
+
+    EllipseAxes axes;
+
+    double f = sqrt (0.25*PS_SQR(moments.x2 - moments.y2) + PS_SQR(moments.xy));
+    if (f > (moments.x2 + moments.y2) / 2.0) {
+        f = 0.98*(moments.x2 + moments.y2) / 2.0;
+    }
+
+    axes.major = sqrt (0.5*(moments.x2 + moments.y2) + f);
+    axes.minor = sqrt (0.5*(moments.x2 + moments.y2) - f);
+    axes.theta = atan2 (2*moments.xy, moments.x2 - moments.y2) / 2;
+    // theta in radians
+
+    return (axes);
+}
+
+EllipseShape EllipseAxesToShape (EllipseAxes axes)
+{
+
+    EllipseShape shape;
+
+    double r1 = 1.0 / PS_SQR(axes.major) + 1.0 / PS_SQR(axes.minor);
+    double r2 = 1.0 / PS_SQR(axes.major) - 1.0 / PS_SQR(axes.minor);
+
+    double sxr = r1 + r2*cos(2*axes.theta);
+    double syr = r1 - r2*cos(2*axes.theta);
+
+    shape.sx = 1.0 / sqrt(sxr);
+    shape.sy = 1.0 / sqrt(syr);
+    shape.sxy = r2*sin(2*axes.theta);
+
+    return (shape);
+}
+
+EllipseAxes EllipseShapeToAxes (EllipseShape shape)
+{
+
+    EllipseAxes axes;
+
+    double f1 = 1.0 / PS_SQR(shape.sx) + 1.0 / PS_SQR(shape.sy);
+    double f2 = 1.0 / PS_SQR(shape.sx) - 1.0 / PS_SQR(shape.sy);
+
+    // force the axis ratio to be less than 10
+    double r1 = 0.5*0.95*sqrt (PS_SQR(f1) - PS_SQR(f2));
+
+    shape.sxy = PS_MIN(PS_MAX(shape.sxy, -r1), r1);
+
+    axes.theta = atan2 (-2.0*shape.sxy, f2) / 2.0;
+
+    double Ar = 0.25*f1 + 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
+    double Br = 0.25*f1 - 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
+
+    axes.minor = 1.0 / sqrt (Ar);
+    axes.major = 1.0 / sqrt (Br);
+
+    return (axes);
+}
Index: /trunk/psModules/src/objects/psEllipse.h
===================================================================
--- /trunk/psModules/src/objects/psEllipse.h	(revision 5255)
+++ /trunk/psModules/src/objects/psEllipse.h	(revision 5255)
@@ -0,0 +1,30 @@
+// strucures to define elliptical shape parameters
+typedef struct
+{
+    double major;
+    double minor;
+    double theta;
+}
+EllipseAxes;
+
+typedef struct
+{
+    double x2;
+    double y2;
+    double xy;
+}
+EllipseMoments;
+
+typedef struct
+{
+    double sx;
+    double sy;
+    double sxy;
+}
+EllipseShape;
+
+// conversions between elliptical shape representations
+EllipseAxes EllipseMomentsToAxes (EllipseMoments moments);
+EllipseShape EllipseAxesToShape (EllipseAxes axes);
+EllipseAxes EllipseShapeToAxes (EllipseShape shape);
+
Index: /trunk/psModules/src/objects/psLibUtils.c
===================================================================
--- /trunk/psModules/src/objects/psLibUtils.c	(revision 5255)
+++ /trunk/psModules/src/objects/psLibUtils.c	(revision 5255)
@@ -0,0 +1,509 @@
+# include <strings.h>  // for strncasecmp
+# include <pslib.h>
+# include "psLibUtils.h"
+
+// XXX EAM : this is NOT included in the C99 headers ??
+FILE *fdopen(int fildes, const char *mode);
+
+// XXX EAM : these utility functions should be added back into PSLib
+
+static psHash *timers = NULL;
+
+/* XXX: remove
+bool psTimerClear (char *name) {
+ 
+  bool status;
+ 
+  if (name == NULL) return false;
+ 
+  status = psHashRemove (timers, name);
+  return (status);
+}
+*/
+
+void psTimerFree ()
+{
+
+    psFree (timers);
+    p_psTimeFinalize();
+    return;
+}
+
+// start/restart a named timer
+bool psTimerStart (char *name)
+{
+
+    psTime *start;
+
+    if (timers == NULL)
+        timers = psHashAlloc (16);
+
+    start = psTimeGetNow (PS_TIME_UTC);
+    psHashAdd (timers, name, start);
+    psFree (start);
+    return (TRUE);
+}
+
+// get current elapsed time on named timer
+psF64 psTimerMark (char *name)
+{
+
+    psTime *start;
+    psTime *mark;
+    psF64   delta;
+
+    if (timers == NULL)
+        return (0);
+
+    start = psHashLookup (timers, name);
+    if (start == NULL)
+        return (0);
+
+    mark = psTimeGetNow (PS_TIME_UTC);
+    delta = psTimeDelta (mark, start);
+    psFree (mark);
+    // psFree (start); -- XXX is psHashLookup not psMemCopying?
+
+    return (delta);
+}
+
+/* XXX: remove
+// find the location of the specified argument
+int psArgumentGet (int argc, char **argv, char *arg) {
+ 
+    int i;
+ 
+    for (i = 0; i < argc; i++) {
+ if (!strcmp(argv[i], arg))
+     return (i);
+    }
+  
+    return ((int)NULL);
+}
+*/
+
+/* XXX: remove
+// remove the specified argument (by location)
+int psArgumentRemove (int N, int *argc, char **argv) {
+ 
+    int i;
+ 
+    if ((N != (int)NULL) && (N != 0)) {
+ (*argc)--;
+ for (i = N; i < *argc; i++) {
+     argv[i] = argv[i+1];
+ }
+    }
+ 
+    return (N);
+}
+*/
+
+// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae)
+// 2 = default, -v = 3, -vv = 4, -vvv = 5
+psS32 psLogArguments (int *argc, char **argv)
+{
+
+    int N, level;
+
+    // default log level is 2
+    level = 2;
+
+    // set in order, so that -vvv overrides -vv overrides -v
+    if ((N = psArgumentGet (*argc, argv, "-v"))) {
+        psArgumentRemove (N, argc, argv);
+        level = 3;
+    }
+    if ((N = psArgumentGet (*argc, argv, "-vv"))) {
+        psArgumentRemove (N, argc, argv);
+        level = 4;
+    }
+    if ((N = psArgumentGet (*argc, argv, "-vvv"))) {
+        psArgumentRemove (N, argc, argv);
+        level = 5;
+    }
+
+    if ((N = psArgumentGet (*argc, argv, "-logfmt"))) {
+        if (*argc < N + 2) {
+            psAbort ("psLogArguments", "USAGE: -logfmt (format)");
+        }
+        psArgumentRemove (N, argc, argv);
+        psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid
+        psArgumentRemove (N, argc, argv);
+    }
+
+    // set the level, return the level
+    psLogSetLevel (level);
+    return (level);
+}
+
+// set trace levels by facility
+psS32 psTraceArguments (int *argc, char **argv)
+{
+
+    int N;
+
+    // argument format is: -trace (facil) (level)
+    while ((N = psArgumentGet (*argc, argv, "-trace"))) {
+        if (*argc < N + 3) {
+            psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)");
+        }
+        psArgumentRemove (N, argc, argv);
+        psTraceSetLevel (argv[N], atoi(argv[N+1]));
+        psArgumentRemove (N, argc, argv);
+        psArgumentRemove (N, argc, argv);
+    }
+    if ((N = psArgumentGet (*argc, argv, "-trace-levels"))) {
+        psTracePrintLevels ();
+        exit (2);
+    }
+    return (TRUE);
+}
+
+# if (0)
+    // alternate implementation of this function from pmObjects.c
+    // now defined in psLib SDRS as psImageRow
+    psVector *psGetRowVectorFromImage(psImage *image,
+                                      psU32 row)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
+
+    psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32);
+    memcpy (tmpVector->data.F32, image->data.F32[row], image->numCols*sizeof(psF32));
+    return(tmpVector);
+}
+# endif
+
+// XXX EAM : this is now in psLib
+void psImageSmooth_EAM (psImage *image, float sigma, float Nsigma)
+{
+
+    int Nx, Ny, Npixel, Nrange;
+    float factor, g, s;
+    psVector *temp;
+
+    // relevant terms
+    Nrange = sigma*Nsigma + 0.5;
+    Npixel = 2*Nrange + 1;
+    factor = -0.5/(sigma*sigma);
+
+    Nx = image->numCols;
+    Ny = image->numRows;
+
+    // generate gaussian
+    psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_F32);
+    for (int i = -Nrange; i < Nrange + 1; i++) {
+        gaussnorm->data.F32[i+Nrange] = exp (factor*i*i);
+    }
+    psF32 *gauss = &gaussnorm->data.F32[Nrange];
+
+    // smooth in X direction
+    temp = psVectorAlloc (Nx, PS_TYPE_F32);
+    for (int j = 0; j < Ny; j++) {
+        psF32 *vi = image->data.F32[j];
+        psF32 *vo = temp->data.F32;
+        for (int i = 0; i < Nx; i++) {
+            g = s = 0;
+            for (int n = -Nrange; n < Nrange + 1; n++) {
+                if (i+n < 0)
+                    continue;
+                if (i+n >= Nx)
+                    continue;
+                s += gauss[n]*vi[i+n];
+                g += gauss[n];
+            }
+            vo[i] = s / g;
+        }
+        memcpy (image->data.F32[j], temp->data.F32, Nx*sizeof(psF32));
+    }
+    psFree (temp);
+
+    // smooth in Y direction
+    temp = psVectorAlloc (image->numRows, PS_TYPE_F32);
+    for (int i = 0; i < Nx; i++) {
+        psF32  *vo = temp->data.F32;
+        psF32 **vi = image->data.F32;
+        for (int j = 0; j < Ny; j++) {
+            g = s = 0;
+            for (int n = -Nrange; n < Nrange + 1; n++) {
+                if (j+n < 0)
+                    continue;
+                if (j+n >= Ny)
+                    continue;
+                s += gauss[n]*vi[j+n][i];
+                g += gauss[n];
+            }
+            vo[j] = s / g;
+        }
+        // replace temp in image
+        for (int j = 0; j < Ny; j++) {
+            vi[j][i] = vo[j];
+        }
+    }
+    psFree (temp);
+    psFree (gaussnorm);
+}
+
+bool psImageInit (psImage *image,...)
+{
+
+    va_list argp;
+    psU8  vU8;
+    psF32 vF32;
+    psF64 vF64;
+
+    if (image == NULL)
+        return (false);
+
+    va_start (argp, image);
+
+    switch (image->type.type) {
+    case PS_TYPE_U8:
+        vU8 = va_arg (argp, psU32);
+
+        for (int iy = 0; iy < image->numRows; iy++) {
+            for (int ix = 0; ix < image->numCols; ix++) {
+                image->data.U8[iy][ix] = vU8;
+            }
+        }
+        break;
+
+    case PS_TYPE_F32:
+        vF32 = va_arg (argp, psF64);
+
+        for (int iy = 0; iy < image->numRows; iy++) {
+            for (int ix = 0; ix < image->numCols; ix++) {
+                image->data.F32[iy][ix] = vF32;
+            }
+        }
+        return (true);
+
+    case PS_TYPE_F64:
+        vF64 = va_arg (argp, psF64);
+
+        for (int iy = 0; iy < image->numRows; iy++) {
+            for (int ix = 0; ix < image->numCols; ix++) {
+                image->data.F64[iy][ix] = vF64;
+            }
+        }
+        return (true);
+
+    default:
+        psAbort ("psphot.psUtils", "datatype %d not defined in psImageInit\n", image->type);
+        return (false);
+    }
+    return (false);
+}
+
+/* XXX: remove
+// count number of pixels with given mask value
+int psImageCountPixelMask (psImage *mask, psU8 value) 
+{
+    int Npixels = 0;
+  
+    for (int i = 0; i < mask->numRows; i++) {
+ for (int j = 0; j < mask->numCols; j++) {
+     if (mask->data.U8[i][j] & value) {
+  Npixels ++;
+     }
+ }
+    }
+    return (Npixels);
+}
+*/
+
+// define a square region centered on the given coordinate
+// XXX EAM : this is now in psLib
+# if (0)
+    psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius)
+{
+    psRegion *region;
+    region = psRegionAlloc (x - radius, x + radius + 1,
+                            y - radius, y + radius + 1);
+    return (region);
+}
+# endif
+
+// set actual region based on image parameters:
+// compensate for negative upper limits
+// XXX this is inconsistent: the coordindates should always be in the parent
+//     frame, which means the negative values should subtract from Nx,Ny of
+//     the parent, not the child.  but, we don't carry the dimensions of the
+//     parent in the psImage container.  for now, use the child Nx,Ny
+//     force range to be on this subimage
+// XXX EAM : this needs to be changed to use psRegion rather than psRegion*
+// XXX EAM : this is now in psLib
+# if (0)
+    psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in)
+{
+
+    // x0,y0, x1,y1 are in *parent* units
+
+    if (out == NULL) {
+        out = psRegionAlloc(in->x0, in->x1, in->y0, in->y1);
+    } else {
+        *out = *in;
+    }
+    // XXX these are probably wrong (see above)
+    if (out->x1 <= 0) {
+        out->x1 = image->col0 + image->numCols + out->x1;
+    }
+    if (out->y1 <= 0) {
+        out->y1 = image->row0 + image->numRows + out->y1;
+    }
+
+    // force the lower-limits to be on the child
+    out->x0 = PS_MAX(image->col0, out->x0);
+    out->y0 = PS_MAX(image->row0, out->y0);
+
+    // force the upper-limits to be on the child
+    out->x1 = PS_MIN(image->col0 + image->numCols, out->x1);
+    out->y1 = PS_MIN(image->row0 + image->numRows, out->y1);
+    return (out);
+}
+# endif
+
+// mask the area contained by the region
+// the region is defined wrt the parent image
+// XXX EAM : this is now in psLib
+# if (0)
+    void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue)
+{
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+        for (int ix = 0; ix < image->numCols; ix++) {
+            if (ix + image->col0 <  region->x0)
+                continue;
+            if (ix + image->col0 >= region->x1)
+                continue;
+            if (iy + image->row0 <  region->y0)
+                continue;
+            if (iy + image->row0 >= region->y1)
+                continue;
+            if (logical_and) {
+                image->data.U8[iy][ix] &= maskValue;
+            } else {
+                image->data.U8[iy][ix] |= maskValue;
+            }
+        }
+    }
+}
+# endif
+
+// mask the area not contained by the region
+// the region is defined wrt the parent image
+// XXX EAM : this is now in psLib
+# if (0)
+    void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue)
+{
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+        for (int ix = 0; ix < image->numCols; ix++) {
+            if (ix + image->col0 <  region->x0)
+                goto maskit;
+            if (ix + image->col0 >= region->x1)
+                goto maskit;
+            if (iy + image->row0 <  region->y0)
+                goto maskit;
+            if (iy + image->row0 >= region->y1)
+                goto maskit;
+            continue;
+maskit:
+            if (logical_and) {
+                image->data.U8[iy][ix] &= maskValue;
+            } else {
+                image->data.U8[iy][ix] |= maskValue;
+            }
+        }
+    }
+}
+# endif
+
+// mask the area contained by the region
+// the region is defined wrt the parent image
+// XXX EAM : this is now in psLib
+# if (0)
+    void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue)
+{
+
+    double dx, dy, r2, R2;
+
+    R2 = PS_SQR(radius);
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+        for (int ix = 0; ix < image->numCols; ix++) {
+            dx = ix + image->col0 - x;
+            dy = iy + image->row0 - y;
+            r2 = PS_SQR(dx) + PS_SQR(dy);
+            if (r2 > R2)
+                continue;
+            if (logical_and) {
+                image->data.U8[iy][ix] &= maskValue;
+            } else {
+                image->data.U8[iy][ix] |= maskValue;
+            }
+        }
+    }
+}
+# endif
+
+// mask the area contained by the region
+// the region is defined wrt the parent image
+// XXX EAM : this is now in psLib
+# if (0)
+    void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue)
+{
+
+    double dx, dy, r2, R2;
+
+    R2 = PS_SQR(radius);
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+        for (int ix = 0; ix < image->numCols; ix++) {
+            dx = ix + image->col0 - x;
+            dy = iy + image->row0 - y;
+            r2 = PS_SQR(dx) + PS_SQR(dy);
+            if (r2 < R2)
+                continue;
+            if (logical_and) {
+                image->data.U8[iy][ix] &= maskValue;
+            } else {
+                image->data.U8[iy][ix] |= maskValue;
+            }
+        }
+    }
+}
+# endif
+
+/* XXX: remove
+psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) {
+ 
+  int nBin = (upper - lower) / delta;
+ 
+  psVector *out = psVectorAlloc (nBin, type);
+  
+  for (int i = 0; i < nBin; i++) {
+    out->data.F64[i] = lower + i * delta;
+  }
+ 
+  return (out);
+}
+*/
+
+// XXX EAM a utility function
+bool p_psVectorPrintRow (int fd, psVector *a, char *name)
+{
+
+    FILE *f;
+    f = fdopen(fd, "a+");
+    fprintf (f, "vector: %s\n", name);
+
+    for (int i = 0; i < a[0].n; i++) {
+        fprintf (f, "%f  ", p_psVectorGetElementF64(a, i));
+    }
+    fprintf (f, "\n");
+    fclose(f);
+    return (true);
+}
+// XXX EAM is the use of fdopen here a good way to do this?
Index: /trunk/psModules/src/objects/psLibUtils.h
===================================================================
--- /trunk/psModules/src/objects/psLibUtils.h	(revision 5255)
+++ /trunk/psModules/src/objects/psLibUtils.h	(revision 5255)
@@ -0,0 +1,59 @@
+
+# ifndef PS_LIB_UTILS
+# define PS_LIB_UTILS
+
+// XXX EAM : psEllipse needs to be be included in psLib
+# include "psEllipse.h"
+
+// structure to carry a dynamic string
+typedef struct
+{
+    int NLINE;
+    int Nline;
+    char *line;
+}
+psLine;
+
+# define psMemCopy(A)(psMemIncrRefCounter((A)))
+
+// XXX EAM : my version using varience instead of stdev
+bool psMinimizeGaussNewtonDelta_EAM (psVector *delta,
+                                     const psVector *params,
+                                     const psVector *paramMask,
+                                     const psArray  *x,
+                                     const psVector *y,
+                                     const psVector *yErr,
+                                     psMinimizeLMChi2Func func);
+
+// minimize : using varience vs sigma and parameter limits
+psBool       p_psMinLM_GuessABP_EAM (psImage  *Alpha, psVector *Beta, psVector *Params, const psImage  *alpha, const psVector *beta, const psVector *params, const psVector *paramMask, const psVector *beta_lim, const psVector *params_min, const psVector *params_max, psF64 lambda);
+psBool       psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func);
+psF64        p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda);
+
+// psLib extra utilities
+bool       psTimerStart (char *name);   // added to SDRS
+//XXX: I removed this: bool       psTimerClear (char *name);   // added to SDRS
+psF64       psTimerMark (char *name);    // added to SDRS
+void       psTimerFree ();              // added to SDRS (as psTimerStop)
+psS32       psLogArguments (int *argc, char **argv);   // added to SDRS (part of psArgumentVerbosity)
+psS32       psTraceArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity)
+//XXX: remove: int      psArgumentGet (int argc, char **argv, char *arg); // added to SDRS
+//XXX: remove: int      psArgumentRemove (int N, int *argc, char **argv); // added to SDRS
+//XXX: remove: psVector    *psVectorCreate (double lower, double upper, double delta, psElemType type); // added to SDRS
+// psVector    *psGetRowVectorFromImage(psImage *image, psU32 row); // added to SDRS (as psImageRow)
+//XXX: remove: int          psImageCountPixelMask (psImage *mask, psU8 value); // added to SDRS
+
+// basic image functions
+bool         psImageInit (psImage *image,...); // added to SDRS (v.15)
+void      psImageSmooth_EAM (psImage *image, float sigma, float Nsigma); // added to SDRS (v.15)
+
+// psLine functions -- keep out for now?
+psLine      *psLineAlloc (int Nline);
+bool      psLineInit (psLine *line);
+bool      psLineAdd (psLine *line, char *format, ...);
+
+// not included in the .h file -- these are fairly lame, keep out?
+bool p_psVectorPrint (int fd, psVector *a, char *name);
+bool p_psVectorPrintRow (int fd, psVector *a, char *name);
+
+# endif
Index: /trunk/psModules/src/objects/psModulesUtils.c
===================================================================
--- /trunk/psModules/src/objects/psModulesUtils.c	(revision 5255)
+++ /trunk/psModules/src/objects/psModulesUtils.c	(revision 5255)
@@ -0,0 +1,95 @@
+# include "psModulesUtils.h"
+
+// XXXX: Must figure out the right names for these defines, then replace!
+#define PS_META_STR 0
+#define PS_META_F32 0
+
+
+
+// extract config informatin from config data or from image header, if specified
+// XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme
+psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name)
+{
+
+    char *source;
+    char *keyword;
+    psF32 value;
+    psMetadataItem *item;
+
+    // look for the entry in the config collection
+    item = psMetadataLookup (config, name);
+    if (item == NULL) {
+        psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name);
+        value = psMetadataLookupF32 (status, header, name);
+        if (*status == false) {
+            psAbort ("pmConfigLookupF32", "no key %s in header", name);
+        }
+        *status = true;
+        return (value);
+    }
+
+    // I'm either expecting a string, with the name "HD:keyword"...
+    if (item->type == PS_META_STR) {
+        source = item->data.V;
+        if (!strncasecmp (source, "HD:", 3)) {
+            keyword = &source[3];
+            value = psMetadataLookupF32 (status, header, keyword);
+            if (*status == false) {
+                psAbort ("pmConfigLookupF32", "no key %s in config", name);
+            }
+            *status = true;
+            // psFree (item);
+            return (value);
+        }
+    }
+
+    //  ... or a value (F32?)
+    if (item->type == PS_META_F32) {
+        value = item->data.F32;
+        // psFree (item);
+        return (value);
+    }
+
+    *status = false;
+    psError(PS_ERR_UNKNOWN, true, "invalid item");
+    return (0);
+}
+
+bool pmModelFitStatus (pmModel *model)
+{
+
+    bool status;
+
+    pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type);
+    status = statusFunc (model);
+
+    return (status);
+}
+
+bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask)
+{
+
+    float obsSum = 0;
+    float fitSum = 0;
+    float sky = model->params->data.F32[0];
+
+    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
+    fitSum = modelFluxFunc (model->params);
+
+    for (int ix = 0; ix < image->numCols; ix++) {
+        for (int iy = 0; iy < image->numRows; iy++) {
+            if (mask->data.U8[iy][ix])
+                continue;
+            obsSum += image->data.F32[iy][ix] - sky;
+        }
+    }
+    if (obsSum <= 0)
+        return false;
+    if (fitSum <= 0)
+        return false;
+
+    *fitMag = -2.5*log10(fitSum);
+    *obsMag = -2.5*log10(obsSum);
+    return (true);
+}
+
Index: /trunk/psModules/src/objects/psModulesUtils.h
===================================================================
--- /trunk/psModules/src/objects/psModulesUtils.h	(revision 5255)
+++ /trunk/psModules/src/objects/psModulesUtils.h	(revision 5255)
@@ -0,0 +1,20 @@
+
+#ifndef PS_MODULE_UTILS
+#define PS_MODULE_UTILS
+
+#include <strings.h>  // for strcasecmp
+#include <pslib.h>
+#include "psLibUtils.h"
+#include "pmObjects.h"
+#include "pmModelGroup.h"
+
+// psModule extra utilities
+// bool      pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF);
+bool       pmModelFitStatus (pmModel *model);
+int      pmSourceDophotType (pmSource *source);
+bool      pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask);
+
+// XXX: unify with paul's image/header/metadata functions
+psF32        pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
+
+#endif
Index: /trunk/psModules/test/objects/tst_pmObjects01.c
===================================================================
--- /trunk/psModules/test/objects/tst_pmObjects01.c	(revision 5254)
+++ /trunk/psModules/test/objects/tst_pmObjects01.c	(revision 5255)
@@ -17,8 +17,14 @@
  *
  * XXX: Memory leaks are not being caught.  If I allocated a psVector in these functions
- * abd never deallocate, no error is generated.
+ * and never deallocate, no error is generated.
+Fully Tested:
+    pmPeakAlloc()
+    pmMomentsAlloc()
+Weakly Tested:
+    pmSourceMoments()
+ 
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-28 20:42:52 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-10-10 19:53:54 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27,4 +33,6 @@
 #include "pslib.h"
 #include "pmObjects.h"
+#include "pmModelGroup.h"
+
 #define NUM_ROWS 10
 #define NUM_COLS 10
@@ -37,6 +45,7 @@
 static int test04(void);
 static int test05(void);
-static int test06(void);
+//static int test06(void);
 static int test07(void);
+/*
 static int test08(void);
 static int test09(void);
@@ -44,4 +53,5 @@
 static int test16(void);
 static int test20(void);
+*/
 testDescription tests[] = {
                               {test00, 000, "pmObjects: structure allocators and deallocators", true, false},
@@ -50,12 +60,14 @@
                               {test03, 001, "pmObjects: pmCullPeaks()", true, false},
                               {test04, 001, "pmObjects: pmSourceLocalSky()", true, false},
-                              {test06, 001, "pmObjects: pmSourceSetPixelsCircle()", true, false},
                               {test05, 001, "pmObjects: pmSourceMoments()", true, false},
+                              //                              {test06, 001, "pmObjects: pmSourceSetPixelsCircle()", true, false},
                               {test07, 001, "pmObjects: pmMin()", true, false},
-                              {test08, 001, "pmObjects: pmSourceModelGuess()", true, false},
-                              {test09, 001, "pmObjects: pmSourceContour()", true, false},
-                              {test15, 001, "pmObjects: pmSourceAddModel()", true, false},
-                              {test16, 001, "pmObjects: pmSourceSubModel()", true, false},
-                              {test20, 001, "pmObjects: pmSourceSubModel()", true, false},
+                              /*
+                                                            {test08, 001, "pmObjects: pmSourceModelGuess()", true, false},
+                                                            {test09, 001, "pmObjects: pmSourceContour()", true, false},
+                                                            {test15, 001, "pmObjects: pmSourceAddModel()", true, false},
+                                                            {test16, 001, "pmObjects: pmSourceSubModel()", true, false},
+                                                            {test20, 001, "pmObjects: pmSourceSubModel()", true, false},
+                              */
                               {NULL}
                           };
@@ -64,4 +76,41 @@
 {
     psLogSetFormat("HLNM");
+    //
+    // We include the function names here in psTraceSetLevel() commands for
+    // debugging convenience.  There is no guarantee that this list of functions
+    // is complete.
+    //
+    psTraceSetLevel(".", 0);
+    psTraceSetLevel("pmPeakAlloc", 0);
+    psTraceSetLevel("pmMomentsAlloc", 0);
+    psTraceSetLevel("modelFree", 0);
+    psTraceSetLevel("pmModelAlloc", 0);
+    psTraceSetLevel("sourceFree", 0);
+    psTraceSetLevel("pmSourceAlloc", 0);
+    psTraceSetLevel("pmFindVectorPeaks", 0);
+    psTraceSetLevel("getRowVectorFromImage", 0);
+    psTraceSetLevel("myListAddPeak", 0);
+    psTraceSetLevel("pmFindImagePeaks", 0);
+    psTraceSetLevel("isItInThisRegion", 0);
+    psTraceSetLevel("pmCullPeaks", 0);
+    psTraceSetLevel("pmPeaksSubset", 0);
+    psTraceSetLevel("pmSourceLocalSky", 0);
+    psTraceSetLevel("checkRadius2", 0);
+    psTraceSetLevel("pmSourceMoments", 0);
+    psTraceSetLevel("pmComparePeakAscend", 0);
+    psTraceSetLevel("pmComparePeakDescend", 0);
+    psTraceSetLevel("pmSourcePSFClump", 0);
+    psTraceSetLevel("pmSourceRoughClass", 0);
+    psTraceSetLevel("pmSourceDefinePixels", 0);
+    psTraceSetLevel("pmSourceModelGuess", 0);
+    psTraceSetLevel("pmModelEval", 0);
+    psTraceSetLevel("findValue", 0);
+    psTraceSetLevel("pmSourceContour", 0);
+    psTraceSetLevel("pmSourceFitModel_v5", 0);
+    psTraceSetLevel("pmSourceFitModel", 0);
+    psTraceSetLevel("p_pmSourceAddOrSubModel", 0);
+    psTraceSetLevel("pmSourceAddModel", 0);
+    psTraceSetLevel("pmSourceSubModel", 0);
+
     return !runTestSuite(stderr, "Test Point Driver", tests, argc, argv);
 }
@@ -116,4 +165,13 @@
                 (tmpMoments->nPixels != 0)) {
             printf("TEST ERROR: pmMomentsAlloc() did not properly initialize the pmMoments structure.\n");
+            printf("    tmpMoments->x is %f\n", tmpMoments->x);
+            printf("    tmpMoments->y is %f\n", tmpMoments->y);
+            printf("    tmpMoments->Sx is %f\n", tmpMoments->Sx);
+            printf("    tmpMoments->Sy is %f\n", tmpMoments->Sy);
+            printf("    tmpMoments->Sxy is %f\n", tmpMoments->Sxy);
+            printf("    tmpMoments->Sum is %f\n", tmpMoments->Sum);
+            printf("    tmpMoments->Peak is %f\n", tmpMoments->Peak);
+            printf("    tmpMoments->Sky is %f\n", tmpMoments->Sky);
+            printf("    tmp    Moments->nPixels is %d\n", tmpMoments->nPixels);
             testStatus = false;
         }
@@ -121,135 +179,37 @@
     psFree(tmpMoments);
 
-    printf("Testing pmModelAlloc(PS_MODEL_GAUSS)...\n");
-    pmModel *tmpModel = pmModelAlloc(PS_MODEL_GAUSS);
-    if (tmpModel == NULL) {
-        printf("TEST ERROR: pmModelAlloc(PS_MODEL_GAUSS) returned a NULL pmModel\n");
-        testStatus = false;
-    } else {
-        if ((tmpModel->params->n != 7) || (tmpModel->dparams->n != 7)) {
-            printf("TEST ERROR: pmModelAlloc(PS_MODEL_GAUSS) allocated an incorrect number of params (%ld, %ld)\n",
-                   tmpModel->params->n, tmpModel->dparams->n);
+
+    //
+    // Loop through each type of model
+    //
+    psS32 i = 0;
+    while (0 != pmModelParameterCount(i)) {
+        printf("Testing pmModelAlloc(%s)...\n", pmModelGetType(0));
+        pmModel *tmpModel = pmModelAlloc(i);
+        if (tmpModel == NULL) {
+            printf("TEST ERROR: pmModelAlloc(%s) returned a NULL pmModel\n", pmModelGetType(0));
             testStatus = false;
         } else {
-            for (psS32 i = 0 ; i < 7 ; i++) {
-                if ((tmpModel->params->data.F32[i] != 0.0) ||
-                        (tmpModel->dparams->data.F32[i] != 0.0)) {
-                    printf("TEST ERROR: pmModelAlloc(PS_MODEL_GAUSS) did not ininitialize the params/dparams array to 0.0.\n");
-                    testStatus = false;
-                }
-            }
-        }
-    }
-    psFree(tmpModel);
-
-    printf("Testing pmModelAlloc(PS_MODEL_PGAUSS)...\n");
-    tmpModel = pmModelAlloc(PS_MODEL_PGAUSS);
-    if (tmpModel == NULL) {
-        printf("TEST ERROR: pmModelAlloc(PS_MODEL_PGAUSS) returned a NULL pmModel\n");
-        testStatus = false;
-    } else {
-        if ((tmpModel->params->n != 7) || (tmpModel->dparams->n != 7)) {
-            printf("TEST ERROR: pmModelAlloc(PS_MODEL_PGAUSS) allocated an incorrect number of params (%ld, %ld)\n",
-                   tmpModel->params->n, tmpModel->dparams->n);
-            testStatus = false;
-        } else {
-            for (psS32 i = 0 ; i < 7 ; i++) {
-                if ((tmpModel->params->data.F32[i] != 0.0) ||
-                        (tmpModel->dparams->data.F32[i] != 0.0)) {
-                    printf("TEST ERROR: pmModelAlloc(PS_MODEL_PGAUSS) did not ininitialize the params/dparams array to 0.0.\n");
-                    testStatus = false;
-                }
-            }
-        }
-    }
-    psFree(tmpModel);
-
-    printf("Testing pmModelAlloc(PS_MODEL_TWIST_GAUSS)...\n");
-    tmpModel = pmModelAlloc(PS_MODEL_TWIST_GAUSS);
-    if (tmpModel == NULL) {
-        printf("TEST ERROR: pmModelAlloc(PS_MODEL_TWIST_GAUSS) returned a NULL pmModel\n");
-        testStatus = false;
-    } else {
-        if ((tmpModel->params->n != 11) || (tmpModel->dparams->n != 11)) {
-            printf("TEST ERROR: pmModelAlloc(PS_MODEL_TWIST_GAUSS) allocated an incorrect number of params (%ld, %ld)\n",
-                   tmpModel->params->n, tmpModel->dparams->n);
-            testStatus = false;
-        } else {
-            for (psS32 i = 0 ; i < 11 ; i++) {
-                if ((tmpModel->params->data.F32[i] != 0.0) ||
-                        (tmpModel->dparams->data.F32[i] != 0.0)) {
-                    printf("TEST ERROR: pmModelAlloc(PS_MODEL_TWIST_GAUSS) did not ininitialize the params/dparams array to 0.0.\n");
-                    testStatus = false;
-                }
-            }
-        }
-    }
-    psFree(tmpModel);
-
-    printf("Testing pmModelAlloc(PS_MODEL_WAUSS)...\n");
-    tmpModel = pmModelAlloc(PS_MODEL_WAUSS);
-    if (tmpModel == NULL) {
-        printf("TEST ERROR: pmModelAlloc(PS_MODEL_WAUSS) returned a NULL pmModel\n");
-        testStatus = false;
-    } else {
-        if ((tmpModel->params->n != 9) || (tmpModel->dparams->n != 9)) {
-            printf("TEST ERROR: pmModelAlloc(PS_MODEL_WAUSS) allocated an incorrect number of params (%ld, %ld)\n",
-                   tmpModel->params->n, tmpModel->dparams->n);
-            testStatus = false;
-        } else {
-            for (psS32 i = 0 ; i < 9 ; i++) {
-                if ((tmpModel->params->data.F32[i] != 0.0) ||
-                        (tmpModel->dparams->data.F32[i] != 0.0)) {
-                    printf("TEST ERROR: pmModelAlloc(PS_MODEL_WAUSS) did not ininitialize the params/dparams array to 0.0.\n");
-                    testStatus = false;
-                }
-            }
-        }
-    }
-    psFree(tmpModel);
-
-    printf("Testing pmModelAlloc(PS_MODEL_SERSIC)...\n");
-    tmpModel = pmModelAlloc(PS_MODEL_SERSIC);
-    if (tmpModel == NULL) {
-        printf("TEST ERROR: pmModelAlloc(PS_MODEL_SERSIC) returned a NULL pmModel\n");
-        testStatus = false;
-    } else {
-        if ((tmpModel->params->n != 8) || (tmpModel->dparams->n != 8)) {
-            printf("TEST ERROR: pmModelAlloc(PS_MODEL_SERSIC) allocated an incorrect number of params (%ld, %ld)\n",
-                   tmpModel->params->n, tmpModel->dparams->n);
-            testStatus = false;
-        } else {
-            for (psS32 i = 0 ; i < 8 ; i++) {
-                if ((tmpModel->params->data.F32[i] != 0.0) ||
-                        (tmpModel->dparams->data.F32[i] != 0.0)) {
-                    printf("TEST ERROR: pmModelAlloc(PS_MODEL_SERSIC) did not ininitialize the params/dparams array to 0.0.\n");
-                    testStatus = false;
-                }
-            }
-        }
-    }
-    psFree(tmpModel);
-
-    printf("Testing pmModelAlloc(PS_MODEL_SERSIC_CORE)...\n");
-    tmpModel = pmModelAlloc(PS_MODEL_SERSIC_CORE);
-    if (tmpModel == NULL) {
-        printf("TEST ERROR: pmModelAlloc(PS_MODEL_SERSIC_CORE) returned a NULL pmModel\n");
-        testStatus = false;
-    } else {
-        if ((tmpModel->params->n != 12) || (tmpModel->dparams->n != 12)) {
-            printf("TEST ERROR: pmModelAlloc(PS_MODEL_SERSIC_CORE) allocated an incorrect number of params (%ld, %ld)\n",
-                   tmpModel->params->n, tmpModel->dparams->n);
-            testStatus = false;
-        } else {
-            for (psS32 i = 0 ; i < 12 ; i++) {
-                if ((tmpModel->params->data.F32[i] != 0.0) ||
-                        (tmpModel->dparams->data.F32[i] != 0.0)) {
-                    printf("TEST ERROR: pmModelAlloc(PS_MODEL_SERSIC_CORE) did not ininitialize the params/dparams array to 0.0.\n");
-                    testStatus = false;
-                }
-            }
-        }
-    }
-    psFree(tmpModel);
+
+            /* XXX: Should we test that the members were set correctly?
+                        if ((tmpModel->params->n != 7) || (tmpModel->dparams->n != 7)) {
+                            printf("TEST ERROR: pmModelAlloc(PS_MODEL_GAUSS) allocated an incorrect number of params (%ld, %ld)\n",
+                                   tmpModel->params->n, tmpModel->dparams->n);
+                            testStatus = false;
+                        } else {
+                            for (psS32 i = 0 ; i < 7 ; i++) {
+                                if ((tmpModel->params->data.F32[i] != 0.0) ||
+                                        (tmpModel->dparams->data.F32[i] != 0.0)) {
+                                    printf("TEST ERROR: pmModelAlloc(PS_MODEL_GAUSS) did not ininitialize the params/dparams array to 0.0.\n");
+                                    testStatus = false;
+                                }
+                            }
+                        }
+            */
+        }
+        psFree(tmpModel);
+        i++;
+    }
+
 
     pmSource *tmpSource = pmSourceAlloc();
@@ -775,4 +735,70 @@
 }
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 #define TST04_NUM_ROWS 100
 #define TST04_NUM_COLS 100
@@ -792,111 +818,35 @@
     bool testStatus = true;
     psImage *imgData = psImageAlloc(TST04_NUM_COLS, TST04_NUM_ROWS, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < imgData->numRows; i++) {
-        for (psS32 j = 0 ; j < imgData->numCols; j++) {
-            imgData->data.F32[i][j] = TST04_SKY;
-        }
-    }
+    psImageInit(imgData, TST04_SKY);
+    psImage *imgMask = psImageAlloc(TST04_NUM_COLS, TST04_NUM_ROWS, PS_TYPE_U8);
+    psImageInit(imgMask, 0);
+    //    psImage *imgMaskS8 = psImageAlloc(TST04_NUM_COLS, TST04_NUM_ROWS, PS_TYPE_S8);
+    //    psImageInit(imgMaskS8, 0);
     psImage *imgDataF64 = psImageAlloc(TST04_NUM_COLS, TST04_NUM_ROWS, PS_TYPE_F64);
-    for (psS32 i = 0 ; i < imgDataF64->numRows; i++) {
-        for (psS32 j = 0 ; j < imgDataF64->numCols; j++) {
-            imgDataF64->data.F64[i][j] = 0.0;
-        }
-    }
-    pmSource *rc = NULL;
+    psImageInit(imgDataF64, 0.0);
     pmPeak *tmpPeak = pmPeakAlloc((psF32) (TST04_NUM_ROWS / 2),
                                   (psF32) (TST04_NUM_COLS / 2),
                                   200.0,
                                   PM_PEAK_LONE);
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceLocalSky with NULL psImage.  Should generate error and return NULL.\n");
-    rc = pmSourceLocalSky(NULL, tmpPeak, PS_STAT_SAMPLE_MEAN, 10.0, 20.0);
-    if (rc != NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a non-NULL pmSource.\n");
-        psFree(rc);
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceLocalSky with wrong-type psImage.  Should generate error and return NULL.\n");
-    rc = pmSourceLocalSky(imgDataF64, tmpPeak, PS_STAT_SAMPLE_MEAN, 10.0, 20.0);
-    if (rc != NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a non-NULL pmSource.\n");
-        psFree(rc);
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceLocalSky with NULL pmPeak.  Should generate error and return NULL.\n");
-    rc = pmSourceLocalSky(imgData, NULL, PS_STAT_SAMPLE_MEAN, 10.0, 20.0);
-    if (rc != NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a non-NULL pmSource.\n");
-        psFree(rc);
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceLocalSky with innerRadius<0.0.  Should generate error and return NULL.\n");
-    rc = pmSourceLocalSky(imgData, tmpPeak, PS_STAT_SAMPLE_MEAN, -10.0, 20.0);
-    if (rc != NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a non-NULL pmSource.\n");
-        psFree(rc);
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceLocalSky with innerRadius>outerRadius.  Should generate error and return NULL.\n");
-    rc = pmSourceLocalSky(imgData, tmpPeak, PS_STAT_SAMPLE_MEAN, 10.0, 5.0);
-    if (rc != NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a non-NULL pmSource.\n");
-        psFree(rc);
-        testStatus = false;
-    }
-    /* XXX: This is commented out since the EAM modification no longer generated NULLS for these tests.
-        printf("----------------------------------------------------------------------------------\n");
-        printf("Calling pmSourceLocalSky with subImage startRow < 0.  Should generate error and return NULL.\n");
-        tmpPeak->x = (psF32) (TST04_NUM_ROWS / 2);
-        tmpPeak->y = (psF32) (TST04_NUM_COLS / 2);
-        tmpPeak->y = 1;
-        rc = pmSourceLocalSky(imgData, tmpPeak, PS_STAT_SAMPLE_MEAN, 10.0, 20.0);
-        if (rc != NULL) {
-            printf("TEST ERROR: pmSourceLocalSky() returned a non-NULL pmSource.\n");
-            psFree(rc);
-            testStatus = false;
-        }
-        printf("----------------------------------------------------------------------------------\n");
-        printf("Calling pmSourceLocalSky with subImage endRow > numRows.  Should generate error and return NULL.\n");
-        tmpPeak->x = (psF32) (TST04_NUM_ROWS / 2);
-        tmpPeak->y = (psF32) (TST04_NUM_COLS / 2);
-        tmpPeak->y = TST04_NUM_ROWS;
-        rc = pmSourceLocalSky(imgData, tmpPeak, PS_STAT_SAMPLE_MEAN, 10.0, 20.0);
-        if (rc != NULL) {
-            printf("TEST ERROR: pmSourceLocalSky() returned a non-NULL pmSource.\n");
-            psFree(rc);
-            testStatus = false;
-        }
-        printf("----------------------------------------------------------------------------------\n");
-        printf("Calling pmSourceLocalSky with subImage startCol < 0.  Should generate error and return NULL.\n");
-        tmpPeak->x = (psF32) (TST04_NUM_ROWS / 2);
-        tmpPeak->y = (psF32) (TST04_NUM_COLS / 2);
-        tmpPeak->x = 1;
-        rc = pmSourceLocalSky(imgData, tmpPeak, PS_STAT_SAMPLE_MEAN, 10.0, 20.0);
-        if (rc != NULL) {
-            printf("TEST ERROR: pmSourceLocalSky() returned a non-NULL pmSource.\n");
-            psFree(rc);
-            testStatus = false;
-        }
-        printf("----------------------------------------------------------------------------------\n");
-        printf("Calling pmSourceLocalSky with subImage endCol > numCols.  Should generate error and return NULL.\n");
-        tmpPeak->x = (psF32) (TST04_NUM_ROWS / 2);
-        tmpPeak->y = (psF32) (TST04_NUM_COLS / 2);
-        tmpPeak->x = TST04_NUM_COLS;
-        rc = pmSourceLocalSky(imgData, tmpPeak, PS_STAT_SAMPLE_MEAN, (psF32) TST04_INNER_RADIUS, (psF32) TST04_OUTER_RADIUS);
-        if (rc != NULL) {
-            printf("TEST ERROR: pmSourceLocalSky() returned a non-NULL pmSource.\n");
-            psFree(rc);
-            testStatus = false;
-        }
-    */
+    pmSource *tmpSource = pmSourceAlloc();
+    tmpSource->pixels = imgData;
+    tmpSource->mask = imgMask;
+    tmpSource->peak = tmpPeak;
+
+    printf("----------------------------------------------------------------------------------\n");
+    printf("Calling pmSourceLocalSky with NULL tmpSource.  Should generate error and return FALSE.\n");
+    bool rc = pmSourceLocalSky(NULL, PS_STAT_SAMPLE_MEAN, 10.0);
+    if (rc != false) {
+        printf("TEST ERROR: pmSourceLocalSky() returned a non-FALSE pmSource.\n");
+        testStatus = false;
+    }
+
+    printf("----------------------------------------------------------------------------------\n");
+    printf("Calling pmSourceLocalSky with Radius<0.0.  Should generate error and return FALSE.\n");
+    rc = pmSourceLocalSky(tmpSource, PS_STAT_SAMPLE_MEAN, -10.0);
+    if (rc != false) {
+        printf("TEST ERROR: pmSourceLocalSky() returned a non-FALSE pmSource.\n");
+        testStatus = false;
+    }
 
     //
@@ -908,294 +858,57 @@
     tmpPeak->x = (psF32) (TST04_NUM_ROWS / 2);
     tmpPeak->y = (psF32) (TST04_NUM_COLS / 2);
-    rc = pmSourceLocalSky(imgData,
-                          tmpPeak,
-                          PS_STAT_SAMPLE_MEAN,
-                          (psF32) TST04_INNER_RADIUS,
-                          (psF32) TST04_OUTER_RADIUS);
-
-    if (rc == NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource.\n");
+    rc = pmSourceLocalSky(tmpSource, PS_STAT_SAMPLE_MEAN, 10.0);
+
+    if (rc == false) {
+        printf("TEST ERROR: pmSourceLocalSky() returned a FALSE pmSource.\n");
         testStatus = false;
     } else {
-        if (rc->peak == NULL) {
+        if (tmpSource->peak == NULL) {
             printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource->peak.\n");
             testStatus = false;
         } else {
-            if (rc->peak->x != tmpPeak->x) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->peak->x was %d, should have been %d.\n", rc->peak->x, tmpPeak->x);
+            if (tmpSource->peak->x != tmpPeak->x) {
+                printf("TEST ERROR: pmSourceLocalSky() pmSource->peak->x was %d, should have been %d.\n",
+                       tmpSource->peak->x, tmpPeak->x);
                 testStatus = false;
             }
 
-            if (rc->peak->y != tmpPeak->y) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->peak->y was %d, should have been %d.\n", rc->peak->y, tmpPeak->y);
+            if (tmpSource->peak->y != tmpPeak->y) {
+                printf("TEST ERROR: pmSourceLocalSky() pmSource->peak->y was %d, should have been %d.\n",
+                       tmpSource->peak->y, tmpPeak->y);
                 testStatus = false;
             }
 
-            if (rc->peak->counts != tmpPeak->counts) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->peak->counts was %f, should have been %f.\n", rc->peak->counts, tmpPeak->counts);
+            if (tmpSource->peak->counts != tmpPeak->counts) {
+                printf("TEST ERROR: pmSourceLocalSky() pmSource->peak->counts was %f, should have been %f.\n",
+                       tmpSource->peak->counts, tmpPeak->counts);
                 testStatus = false;
             }
 
-            if (rc->peak->class != tmpPeak->class) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->peak->class was %d, should have been %d.\n", rc->peak->class, tmpPeak->class);
+            if (tmpSource->peak->class != tmpPeak->class) {
+                printf("TEST ERROR: pmSourceLocalSky() pmSource->peak->class was %d, should have been %d.\n",
+                       tmpSource->peak->class, tmpPeak->class);
                 testStatus = false;
             }
         }
 
-        if (rc->pixels == NULL) {
-            printf("TEST ERROR: pmSourceLocalSky() pmSource->pixels was NULL.\n");
+        if (tmpSource->moments == NULL) {
+            printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource->moments.\n");
             testStatus = false;
         } else {
-            if (rc->pixels->numRows != (2 * TST04_OUTER_RADIUS)) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->pixels->numRows was %d, should have been %d.\n",
-                       rc->pixels->numRows, (2 * TST04_OUTER_RADIUS));
+            if (tmpSource->moments->Sky != TST04_SKY) {
+                printf("TEST ERROR: pmSourceLocalSky() pmSource->moments->Sky was %f, should have been %f.\n", tmpSource->moments->Sky, TST04_SKY);
                 testStatus = false;
             }
-
-            if (rc->pixels->numCols != (2 * TST04_OUTER_RADIUS)) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->pixels->numCols was %d, should have been %d.\n",
-                       rc->pixels->numCols, (2 * TST04_OUTER_RADIUS));
-                testStatus = false;
-            }
-
-            if (rc->pixels->col0 != (tmpPeak->x - TST04_OUTER_RADIUS)) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->pixels->col0 was %d, should have been %d.\n",
-                       rc->pixels->col0, (tmpPeak->x - TST04_OUTER_RADIUS));
-                testStatus = false;
-            }
-
-            if (rc->pixels->row0 != (tmpPeak->y - TST04_OUTER_RADIUS)) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->pixels->row0 was %d, should have been %d.\n",
-                       rc->pixels->row0, (tmpPeak->y - TST04_OUTER_RADIUS));
-                testStatus = false;
-            }
-
-            if (rc->pixels->type.type != PS_TYPE_F32) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->pixels->type was %d, should have been %d.\n",
-                       rc->pixels->type.type, PS_TYPE_F32);
-                testStatus = false;
-            }
-
-            // XXX: Test the rc->pixels-> row/col offsets.
-            // XXX: Test that the pixels corresponds to the source image pixels.
-        }
-
-        if (rc->mask == NULL) {
-            printf("TEST ERROR: pmSourceLocalSky() pmSource->mask was NULL.\n");
-            testStatus = false;
-        } else {
-            if (rc->mask->numRows != (2 * TST04_OUTER_RADIUS)) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->mask->numRows was %d, should have been %d.\n",
-                       rc->mask->numRows, (2 * TST04_OUTER_RADIUS));
-                testStatus = false;
-            }
-
-            if (rc->mask->numCols != (2 * TST04_OUTER_RADIUS)) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->mask->numCols was %d, should have been %d.\n",
-                       rc->mask->numCols, (2 * TST04_OUTER_RADIUS));
-                testStatus = false;
-            }
-
-            if (rc->mask->type.type != PS_TYPE_U8) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->mask->type was %d, should have been %d.\n",
-                       rc->mask->type.type, PS_TYPE_U8);
-                testStatus = false;
-            }
-
-            // XXX: Test the rc->mask-> row/col offsets.
-            // XXX: Test that the correct pixels were masked, not merely the number of masked pixels.
-            psS32 unmaskedPixels = 0;
-            psS32 maskedPixels = 0;
-
-            for (psS32 row = 0 ; row < rc->mask->numRows; row++) {
-                for (psS32 col = 0 ; col < rc->mask->numCols; col++) {
-                    if (rc->mask->data.U8[row][col] == 0) {
-                        unmaskedPixels++;
-                    } else {
-                        maskedPixels++;
-                    }
-                }
-            }
-            if (maskedPixels != PS_SQR(2*(TST04_INNER_RADIUS-1))) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->mask had %d masked pixels, should have been %d.\n",
-                       maskedPixels, PS_SQR(2*(TST04_INNER_RADIUS-1)));
-                testStatus = false;
-            }
-            if (unmaskedPixels != (PS_SQR(2*TST04_OUTER_RADIUS) - PS_SQR(2*(TST04_INNER_RADIUS-1)))) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->mask had %d masked pixels, should have been %d.\n",
-                       unmaskedPixels, (PS_SQR(2*TST04_OUTER_RADIUS) - PS_SQR(2*(TST04_INNER_RADIUS-1))));
-                testStatus = false;
-            }
-        }
-
-        if (rc->moments == NULL) {
-            printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource->moments.\n");
-            testStatus = false;
-        } else {
-            if (rc->moments->Sky != TST04_SKY) {
-                printf("TEST ERROR: pmSourceLocalSky() pmSource->moments->Sky was %f, should have been %f.\n", rc->moments->Sky, TST04_SKY);
-                testStatus = false;
-            }
-        }
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    psFree(rc);
-    psFree(imgData);
-    psFree(imgDataF64);
-    return(testStatus);
-}
-
-#define TST06_NUM_ROWS 100
-#define TST06_NUM_COLS 100
-#define TST06_SKY 20.0
-#define TST06_INNER_RADIUS 3
-#define TST06_OUTER_RADIUS 5
-/******************************************************************************
-test06(): We first test pmSourceLocalSky() with various NULL and unallowable
-input parameters.
- 
-XXX: Should we produce tests with boundary numbers for the inner/outer radius?
- 
-XXX: Call this with varying sizes for the image.
- *****************************************************************************/
-int test06( void )
-{
-    bool testStatus = true;
-    pmSource *tmpSource = NULL;
-    bool rc = false;
-    // Create the image used in this test.
-    psImage *imgData = psImageAlloc(TST06_NUM_COLS, TST06_NUM_ROWS, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < imgData->numRows; i++) {
-        for (psS32 j = 0 ; j < imgData->numCols; j++) {
-            imgData->data.F32[i][j] = TST06_SKY;
-        }
-    }
-    psImage *imgDataF64 = psImageAlloc(TST06_NUM_COLS, TST06_NUM_ROWS, PS_TYPE_F64);
-    for (psS32 i = 0 ; i < imgDataF64->numRows; i++) {
-        for (psS32 j = 0 ; j < imgDataF64->numCols; j++) {
-            imgDataF64->data.F64[i][j] = 0.0;
-        }
-    }
-
-    //
-    // Create a pmPeak with the center pixel set to the peak.
-    //
-    pmPeak *tmpPeak = pmPeakAlloc((psF32) (TST06_NUM_ROWS / 2),
-                                  (psF32) (TST06_NUM_COLS / 2),
-                                  200.0,
-                                  PM_PEAK_LONE);
-
-    tmpSource = pmSourceLocalSky(imgData,
-                                 tmpPeak,
-                                 PS_STAT_SAMPLE_MEAN,
-                                 (psF32) TST06_INNER_RADIUS,
-                                 (psF32) TST06_OUTER_RADIUS);
-    if (tmpSource == NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource.\n");
-        psFree(imgData);
-        psFree(imgDataF64);
-        psFree(tmpSource);
-        testStatus = false;
-        return(testStatus);
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceSetPixelsCircle with NULL pmSource.  Should generate error and return NULL.\n");
-    rc = pmSourceSetPixelsCircle(NULL, imgData, 10.0);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceSetPixelsCircle() returned TRUE.\n");
-        testStatus = false;
-    }
-    // XXX: test with pmSource->peaks NULL
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceSetPixelsCircle with NULL psImage.  Should generate error and return NULL.\n");
-    rc = pmSourceSetPixelsCircle(tmpSource, NULL, 10.0);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceSetPixelsCircle() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceSetPixelsCircle with wrong type psImage.  Should generate error and return NULL.\n");
-    rc = pmSourceSetPixelsCircle(tmpSource, imgDataF64, 10.0);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceSetPixelsCircle() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceSetPixelsCircle with radius < 0.0.  Should generate error and return NULL.\n");
-    rc = pmSourceSetPixelsCircle(tmpSource, imgData, -10.0);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceSetPixelsCircle() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    /* XXX: Commented away since the EAM mods no longer produced errors.
-        printf("----------------------------------------------------------------------------------\n");
-        printf("Calling pmSourceSetPixelsCircle with subImage startCol < 0.  Should generate error and return NULL.\n");
-        tmpSource->peak->x = (psF32) (TST06_NUM_ROWS / 2);
-        tmpSource->peak->y = (psF32) (TST06_NUM_COLS / 2);
-        tmpSource->peak->x = 1;
-        rc = pmSourceSetPixelsCircle(tmpSource, imgData, 10.0);
-        if (rc == true) {
-            printf("TEST ERROR: pmSourceSetPixelsCircle() returned TRUE.\n");
-            testStatus = false;
-        }
-
-        printf("----------------------------------------------------------------------------------\n");
-        printf("Calling pmSourceSetPixelsCircle with subImage endCol > numCols.  Should generate error and return NULL.\n");
-        tmpSource->peak->x = (psF32) (TST06_NUM_ROWS / 2);
-        tmpSource->peak->y = (psF32) (TST06_NUM_COLS / 2);
-        tmpSource->peak->x = TST06_NUM_COLS;
-        rc = pmSourceSetPixelsCircle(tmpSource, imgData, 10.0);
-        if (rc == true) {
-            printf("TEST ERROR: pmSourceSetPixelsCircle() returned TRUE.\n");
-            testStatus = false;
-        }
-
-        printf("----------------------------------------------------------------------------------\n");
-        printf("Calling pmSourceSetPixelsCircle with subImage startRow < 0.  Should generate error and return NULL.\n");
-        tmpSource->peak->x = (psF32) (TST06_NUM_ROWS / 2);
-        tmpSource->peak->y = (psF32) (TST06_NUM_COLS / 2);
-        tmpSource->peak->y = 1;
-        rc = pmSourceSetPixelsCircle(tmpSource, imgData, 10.0);
-        if (rc == true) {
-            printf("TEST ERROR: pmSourceSetPixelsCircle() returned TRUE.\n");
-            testStatus = false;
-        }
-
-        printf("----------------------------------------------------------------------------------\n");
-        printf("Calling pmSourceSetPixelsCircle with subImage endRow > numRows.  Should generate error and return NULL.\n");
-        tmpSource->peak->x = (psF32) (TST06_NUM_ROWS / 2);
-        tmpSource->peak->y = (psF32) (TST06_NUM_COLS / 2);
-        tmpSource->peak->y = TST06_NUM_ROWS;
-        rc = pmSourceSetPixelsCircle(tmpSource, imgData, 10.0);
-        if (rc == true) {
-            printf("TEST ERROR: pmSourceSetPixelsCircle() returned TRUE.\n");
-            testStatus = false;
-        }
-    */
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceSetPixelsCircle with valid data.\n");
-    tmpSource->peak->x = (psF32) (TST06_NUM_ROWS / 2);
-    tmpSource->peak->y = (psF32) (TST06_NUM_COLS / 2);
-    rc = pmSourceSetPixelsCircle(tmpSource, imgData, 10.0);
-
-    if (rc == false) {
-        printf("TEST ERROR: pmSourceSetPixelsCircle() returned FALSE.\n");
-        testStatus = false;
-    } else {
-        // XXX: Test correctness of the various psSetPixelsCircle members.
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
+        }
+    }
+
+    printf("----------------------------------------------------------------------------------\n");
+    psFree(tmpSource);
     //    psFree(imgData);
     //    psFree(imgDataF64);
-    psFree(tmpSource);
+    //    psFree(imgMask);
+    //    psFree(imgMaskS8);
     return(testStatus);
-
 }
 
@@ -1218,713 +931,48 @@
 {
     bool testStatus = true;
-    pmSource *tmpSource = NULL;
-    pmSource *rc = NULL;
-    // Create the image used in this test.
-    psImage *imgData = psImageAlloc(TST05_NUM_COLS, TST05_NUM_ROWS, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < imgData->numRows; i++) {
-        for (psS32 j = 0 ; j < imgData->numCols; j++) {
-            imgData->data.F32[i][j] = TST05_SKY;
-        }
-    }
-
-    //
-    // Create a pmPeak with the center pixel set to the peak.
-    //
-    pmPeak *tmpPeak = pmPeakAlloc((psF32) (TST05_NUM_ROWS / 2),
-                                  (psF32) (TST05_NUM_COLS / 2),
-                                  200.0,
-                                  PM_PEAK_LONE);
-
-    tmpSource = pmSourceLocalSky(imgData,
-                                 tmpPeak,
-                                 PS_STAT_SAMPLE_MEAN,
-                                 (psF32) TST05_INNER_RADIUS,
-                                 (psF32) TST05_OUTER_RADIUS);
-    if (tmpSource == NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource.\n");
-        psFree(tmpSource);
-        testStatus = false;
-        return(testStatus);
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceMoments with NULL pmSource.  Should generate error and return NULL.\n");
-    rc = pmSourceMoments(NULL, 10.0);
-    if (rc != NULL) {
-        printf("TEST ERROR: pmSourceMoments() returned a non-NULL pmSource.\n");
-        testStatus = false;
-    }
-    // XXX: test with pmSource->peaks NULL
-    // XXX: test with pmSource->pixels NULL
-    // XXX: test with pmSource->mask NULL
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceMoments with radius < 0.0.  Should generate error and return NULL.\n");
-    rc = pmSourceMoments(tmpSource, -10.0);
-    if (rc != NULL) {
-        printf("TEST ERROR: pmSourceMoments() returned a non-NULL pmSource.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    psFree(tmpSource);
-    return(testStatus);
-
-}
-
-/******************************************************************************
-test07(): We first test the various psMinLM_... routines with various NULL and
-unallowable input parameters.
- 
-XXX: We don't verify the numbers.  Must do this.
- *****************************************************************************/
-int test07( void )
-{
-    bool testStatus = true;
-    psF32 rc;
-    psVector *deriv = psVectorAlloc(7, PS_TYPE_F32);
-    psVector *params = psVectorAlloc(7, PS_TYPE_F32);
-    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < 7 ; i++) {
-        deriv->data.F32[i] = 1.0;
-        params->data.F32[i] = 1.0;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_Gauss2D with NULL deriv vector.  Should not generate error.\n");
-    rc = pmMinLM_Gauss2D(NULL, params, x);
-    if (isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_Gauss2D() returned a NAN psF32.\n");
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_Gauss2D with NULL params vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_Gauss2D(deriv, NULL, x);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_Gauss2D() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_Gauss2D with NULL x vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_Gauss2D(deriv, params, NULL);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_Gauss2D() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    psFree(deriv);
-    psFree(params);
-
-    deriv = psVectorAlloc(7, PS_TYPE_F32);
-    params = psVectorAlloc(7, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < 7 ; i++) {
-        deriv->data.F32[i] = 1.0;
-        params->data.F32[i] = 1.0;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_PsuedoGauss2D with NULL deriv vector.  Should not generate error.\n");
-    rc = pmMinLM_PsuedoGauss2D(NULL, params, x);
-    if (isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_PsuedoGauss2D() returned a NAN psF32.\n");
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_PsuedoGauss2D with NULL params vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_PsuedoGauss2D(deriv, NULL, x);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_PsuedoGauss2D() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_PsuedoGauss2D with NULL x vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_PsuedoGauss2D(deriv, params, NULL);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_PsuedoGauss2D() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    psFree(deriv);
-    psFree(params);
-
-
-    deriv = psVectorAlloc(9, PS_TYPE_F32);
-    params = psVectorAlloc(9, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < 9 ; i++) {
-        deriv->data.F32[i] = 1.0;
-        params->data.F32[i] = 1.0;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_Wauss2D with NULL deriv vector.  Should not generate error.\n");
-    rc = pmMinLM_Wauss2D(NULL, params, x);
-    if (isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_Wauss2D() returned a NAN psF32.\n");
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_Wauss2D with NULL params vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_Wauss2D(deriv, NULL, x);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_Wauss2D() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_Wauss2D with NULL x vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_Wauss2D(deriv, params, NULL);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_Wauss2D() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    psFree(deriv);
-    psFree(params);
-
-    deriv = psVectorAlloc(11, PS_TYPE_F32);
-    params = psVectorAlloc(11, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < 11 ; i++) {
-        deriv->data.F32[i] = 1.0;
-        params->data.F32[i] = 1.0;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_TwistGauss2D with NULL deriv vector.  Should not generate error.\n");
-    rc = pmMinLM_TwistGauss2D(NULL, params, x);
-    if (isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_TwistGauss2D() returned a NAN psF32.\n");
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_TwistGauss2D with NULL params vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_TwistGauss2D(deriv, NULL, x);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_TwistGauss2D() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_TwistGauss2D with NULL x vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_TwistGauss2D(deriv, params, NULL);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_TwistGauss2D() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    psFree(deriv);
-    psFree(params);
-
-    deriv = psVectorAlloc(8, PS_TYPE_F32);
-    params = psVectorAlloc(8, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < 8 ; i++) {
-        deriv->data.F32[i] = 1.0;
-        params->data.F32[i] = 1.0;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_Sersic with NULL deriv vector.  Should not generate error.\n");
-    rc = pmMinLM_Sersic(NULL, params, x);
-    if (isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_Sersic() returned a NAN psF32.\n");
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_Sersic with NULL params vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_Sersic(deriv, NULL, x);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_Sersic() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_Sersic with NULL x vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_Sersic(deriv, params, NULL);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_Sersic() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    psFree(deriv);
-    psFree(params);
-
-    deriv = psVectorAlloc(12, PS_TYPE_F32);
-    params = psVectorAlloc(12, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < 12 ; i++) {
-        deriv->data.F32[i] = 1.0;
-        params->data.F32[i] = 1.0;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_SersicCore with NULL deriv vector.  Should not generate error.\n");
-    rc = pmMinLM_SersicCore(NULL, params, x);
-    if (isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_SersicCore() returned a NAN psF32.\n");
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_SersicCore with NULL params vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_SersicCore(deriv, NULL, x);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_SersicCore() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmMinLM_SersicCore with NULL x vector.  Should generate error and return NAN.\n");
-    rc = pmMinLM_SersicCore(deriv, params, NULL);
-    if (!isnan(rc)) {
-        printf("TEST ERROR: pmMinLM_SersicCore() returned a non-NAN psF32 (%f).\n", rc);
-        testStatus = false;
-    }
-
-    psFree(deriv);
-    psFree(params);
-    psFree(x);
-    return(testStatus);
-}
-
-/******************************************************************************
-test08(): We first test pmSourceModelGuess() with various NULL and unallowable
-input parameters.
- 
-XXX: We don't verify the numbers.
- *****************************************************************************/
-int test08( void )
-{
-    bool testStatus = true;
     psImage *imgData = psImageAlloc(TST04_NUM_COLS, TST04_NUM_ROWS, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < imgData->numRows; i++) {
-        for (psS32 j = 0 ; j < imgData->numCols; j++) {
-            imgData->data.F32[i][j] = TST04_SKY;
-        }
-    }
-    pmSource *mySrc = NULL;
-    bool rc = false;
+    psImageInit(imgData, TST04_SKY);
+    psImage *imgMask = psImageAlloc(TST04_NUM_COLS, TST04_NUM_ROWS, PS_TYPE_U8);
+    psImageInit(imgMask, 0);
     pmPeak *tmpPeak = pmPeakAlloc((psF32) (TST04_NUM_ROWS / 2),
                                   (psF32) (TST04_NUM_COLS / 2),
                                   200.0,
                                   PM_PEAK_LONE);
-
-    printf("Calling pmSourceLocalSky with valid data.\n");
-    tmpPeak->x = (psF32) (TST04_NUM_ROWS / 2);
-    tmpPeak->y = (psF32) (TST04_NUM_COLS / 2);
-    mySrc = pmSourceLocalSky(imgData,
-                             tmpPeak,
-                             PS_STAT_SAMPLE_MEAN,
-                             (psF32) TST04_INNER_RADIUS,
-                             (psF32) TST04_OUTER_RADIUS);
-
-    if (mySrc == NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceModelGuess with NULL pmSource.  Should generate error, return FALSE.\n");
-    rc = pmSourceModelGuess(NULL, imgData, PS_MODEL_GAUSS);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceModelGuess() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceModelGuess with NULL psImage.  Should generate error, return FALSE.\n");
-    rc = pmSourceModelGuess(mySrc, NULL, PS_MODEL_GAUSS);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceModelGuess() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceModelGuess with bad model type.  Should generate error, return FALSE.\n");
-    rc = pmSourceModelGuess(mySrc, imgData, PS_MODEL_UNDEFINED);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceModelGuess() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceModelGuess with PS_MODEL_GAUSS\n");
-    rc = pmSourceModelGuess(mySrc, imgData, PS_MODEL_GAUSS);
-    if (rc != true) {
-        printf("TEST ERROR: pmSourceModelGuess() returned FALSE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceModelGuess with PS_MODEL_PGAUSS\n");
-    rc = pmSourceModelGuess(mySrc, imgData, PS_MODEL_PGAUSS);
-    if (rc != true) {
-        printf("TEST ERROR: pmSourceModelGuess() returned FALSE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceModelGuess with PS_MODEL_TWIST_GAUSS\n");
-    rc = pmSourceModelGuess(mySrc, imgData, PS_MODEL_TWIST_GAUSS);
-    if (rc != true) {
-        printf("TEST ERROR: pmSourceModelGuess() returned FALSE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceModelGuess with PS_MODEL_WAUSS\n");
-    rc = pmSourceModelGuess(mySrc, imgData, PS_MODEL_WAUSS);
-    if (rc != true) {
-        printf("TEST ERROR: pmSourceModelGuess() returned FALSE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceModelGuess with PS_MODEL_SERSIC\n");
-    rc = pmSourceModelGuess(mySrc, imgData, PS_MODEL_SERSIC);
-    if (rc != true) {
-        printf("TEST ERROR: pmSourceModelGuess() returned FALSE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceModelGuess with PS_MODEL_SERSIC_CORE\n");
-    rc = pmSourceModelGuess(mySrc, imgData, PS_MODEL_SERSIC_CORE);
-    if (rc != true) {
-        printf("TEST ERROR: pmSourceModelGuess() returned FALSE.\n");
-        testStatus = false;
-    }
-
-    psFree(mySrc);
-    //    psFree(tmpPeak);
-    psFree(imgData);
+    pmSource *tmpSource = pmSourceAlloc();
+    tmpSource->pixels = imgData;
+    tmpSource->mask = imgMask;
+    tmpSource->peak = tmpPeak;
+    psBool rc = pmSourceLocalSky(tmpSource, PS_STAT_SAMPLE_MEAN, 10.0);
+
+    if (rc == false) {
+        printf("TEST ERROR: pmSourceLocalSky() returned a FALSE pmSource.\n");
+        testStatus = false;
+    }
+
+    printf("----------------------------------------------------------------------------------\n");
+    printf("Calling pmSourceMoments with NULL pmSource.  Should generate error and return FALSE.\n");
+    rc = pmSourceMoments(NULL, 10.0);
+    if (rc != false) {
+        printf("TEST ERROR: pmSourceMoments() returned TRUE.\n");
+        testStatus = false;
+    }
+    // XXX: test with pmSource->peaks NULL
+    // XXX: test with pmSource->pixels NULL
+    // XXX: test with pmSource->mask NULL
+
+    printf("----------------------------------------------------------------------------------\n");
+    printf("Calling pmSourceMoments with radius < 0.0.  Should generate error and return FALSE.\n");
+    rc = pmSourceMoments(tmpSource, -10.0);
+    if (rc != false) {
+        printf("TEST ERROR: pmSourceMoments() returned TRUE.\n");
+        testStatus = false;
+    }
+
+    printf("----------------------------------------------------------------------------------\n");
+    psFree(tmpSource);
     return(testStatus);
+
 }
 
-
-#define TST09_NUM_ROWS 70
-#define TST09_NUM_COLS 70
-#define TST09_SKY 5.0
-#define TST09_INNER_RADIUS 3
-#define TST09_OUTER_RADIUS 10
-#define LEVEL (TST09_SKY + 10.0)
-/******************************************************************************
-test09(): We first test pmSourceContour() with various NULL and unallowable
-input parameters.
- 
-XXX: We don't verify the numbers.
- *****************************************************************************/
-int test09( void )
-{
-    bool testStatus = true;
-    psImage *imgData = psImageAlloc(TST09_NUM_COLS, TST09_NUM_ROWS, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < imgData->numRows; i++) {
-        for (psS32 j = 0 ; j < imgData->numCols; j++) {
-            imgData->data.F32[i][j] = TST09_SKY;
-        }
-    }
-    pmSource *mySrc = NULL;
-    psArray *rc = NULL;
-
-    pmPeak *tmpPeak = pmPeakAlloc((psF32) (TST09_NUM_ROWS / 2),
-                                  (psF32) (TST09_NUM_COLS / 2),
-                                  200.0,
-                                  PM_PEAK_LONE);
-
-    printf("Calling pmSourceLocalSky with valid data.\n");
-    tmpPeak->x = (psF32) (TST09_NUM_ROWS / 2);
-    tmpPeak->y = (psF32) (TST09_NUM_COLS / 2);
-    mySrc = pmSourceLocalSky(imgData,
-                             tmpPeak,
-                             PS_STAT_SAMPLE_MEAN,
-                             (psF32) TST09_INNER_RADIUS,
-                             (psF32) TST09_OUTER_RADIUS);
-
-    if (mySrc == NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource.\n");
-        testStatus = false;
-    }
-
-    bool rcBool = pmSourceModelGuess(mySrc, imgData, PS_MODEL_GAUSS);
-    if (rcBool != true) {
-        printf("TEST ERROR: pmSourceModelGuess() returned FALSE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceContour with NULL pmSource .  Should generate error, return NULL.\n");
-    rc = pmSourceContour(NULL, imgData, LEVEL, PS_CONTOUR_CRUDE);
-    if (rc != NULL) {
-        printf("TEST ERROR: pmSourceContour() returned non-NULL.\n");
-        testStatus = false;
-        psFree(rc);
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceContour with NULL psImage .  Should generate error, return NULL.\n");
-    rc = pmSourceContour(mySrc, NULL, LEVEL, PS_CONTOUR_CRUDE);
-    if (rc != NULL) {
-        printf("TEST ERROR: pmSourceContour() returned non-NULL.\n");
-        testStatus = false;
-        psFree(rc);
-    }
-
-    //
-    // XXX: pmSourceContour() has a problem with contour tops/bottoms.
-    // Must correct this.
-    //
-    if (0) {
-        printf("----------------------------------------------------------------------------------\n");
-        printf("Calling pmSourceContour with acceptable data.\n");
-        printf("NOTE: must figure out the parameters for this test to be meaningful.\n");
-        mySrc->modelPSF->params->data.F32[0] = TST09_SKY;
-        mySrc->modelPSF->params->data.F32[1] = 15.0;
-        mySrc->modelPSF->params->data.F32[2] = (psF32) (TST09_NUM_ROWS / 2);
-        mySrc->modelPSF->params->data.F32[3] = (psF32) (TST09_NUM_COLS / 2);
-        mySrc->modelPSF->params->data.F32[4] = 2.0;
-        mySrc->modelPSF->params->data.F32[5] = 2.0;
-        mySrc->modelPSF->params->data.F32[6] = 2.0;
-        rc = pmSourceContour(mySrc, imgData, LEVEL, PS_CONTOUR_CRUDE);
-        if (rc == NULL) {
-            printf("TEST ERROR: pmSourceContour() returned NULL.\n");
-            testStatus = false;
-        } else {
-            psFree(rc);
-        }
-    }
-
-    psFree(mySrc);
-    psFree(imgData);
-    // XXX: This psFree() causes an error.  Why?
-    //psFree(tmpPeak);
-    return(testStatus);
-}
-
-#define TST15_NUM_ROWS 100
-#define TST15_NUM_COLS 100
-#define TST15_SKY 10.0
-#define TST15_INNER_RADIUS 3
-#define TST15_OUTER_RADIUS 5
-/******************************************************************************
-test15(): We first test pmSourceAddModel() with various NULL and unallowable
-input parameters.
- 
-XXX: We don't verify the numbers.
- *****************************************************************************/
-int test15( void )
-{
-    bool testStatus = true;
-    psImage *imgData = psImageAlloc(TST15_NUM_COLS, TST15_NUM_ROWS, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < imgData->numRows; i++) {
-        for (psS32 j = 0 ; j < imgData->numCols; j++) {
-            imgData->data.F32[i][j] = TST15_SKY;
-        }
-    }
-    pmSource *mySrc = NULL;
-    psBool rc = false;
-
-    pmPeak *tmpPeak = pmPeakAlloc((psF32) (TST15_NUM_ROWS / 2),
-                                  (psF32) (TST15_NUM_COLS / 2),
-                                  200.0,
-                                  PM_PEAK_LONE);
-
-    printf("Calling pmSourceLocalSky with valid data.\n");
-    tmpPeak->x = (psF32) (TST15_NUM_ROWS / 2);
-    tmpPeak->y = (psF32) (TST15_NUM_COLS / 2);
-    mySrc = pmSourceLocalSky(imgData,
-                             tmpPeak,
-                             PS_STAT_SAMPLE_MEAN,
-                             (psF32) TST15_INNER_RADIUS,
-                             (psF32) TST15_OUTER_RADIUS);
-
-    if (mySrc == NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource.\n");
-        testStatus = false;
-    }
-
-    mySrc->modelPSF = pmModelAlloc(PS_MODEL_GAUSS);
-    mySrc->modelPSF->params->data.F32[0] = 5.0;
-    mySrc->modelPSF->params->data.F32[1] = 70.0;
-    mySrc->modelPSF->params->data.F32[2] = (psF32) (TST15_NUM_ROWS / 2);
-    mySrc->modelPSF->params->data.F32[3] = (psF32) (TST15_NUM_COLS / 2);
-    mySrc->modelPSF->params->data.F32[4] = 1.0;
-    mySrc->modelPSF->params->data.F32[5] = 1.0;
-    mySrc->modelPSF->params->data.F32[6] = 2.0;
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceAddModel with NULL psImage.  Should generate error, return FALSE.\n");
-    rc = pmSourceAddModel(NULL, mySrc, true);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceAddModel() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceAddModel with NULL psSrc.  Should generate error, return FALSE.\n");
-    rc = pmSourceAddModel(imgData, NULL, true);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceAddModel() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceAddModel with acceptable data.\n");
-    rc = pmSourceAddModel(imgData, mySrc, true);
-    if (rc != true) {
-        printf("TEST ERROR: pmSourceAddModel() returned FALSE.\n");
-        testStatus = false;
-    }
-
-    psFree(mySrc);
-    psFree(imgData);
-    return(testStatus);
-}
-
-#define TST16_NUM_ROWS 100
-#define TST16_NUM_COLS 100
-#define TST16_SKY 10.0
-#define TST16_INNER_RADIUS 3
-#define TST16_OUTER_RADIUS 5
-/******************************************************************************
-test16(): We first test pmSourceSubModel() with various NULL and unallowable
-input parameters.
- 
-XXX: We don't verify the numbers.
- *****************************************************************************/
-int test16( void )
-{
-    bool testStatus = true;
-    psImage *imgData = psImageAlloc(TST16_NUM_COLS, TST16_NUM_ROWS, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < imgData->numRows; i++) {
-        for (psS32 j = 0 ; j < imgData->numCols; j++) {
-            imgData->data.F32[i][j] = TST16_SKY;
-        }
-    }
-    pmSource *mySrc = NULL;
-    psBool rc = false;
-
-    pmPeak *tmpPeak = pmPeakAlloc((psF32) (TST16_NUM_ROWS / 2),
-                                  (psF32) (TST16_NUM_COLS / 2),
-                                  200.0,
-                                  PM_PEAK_LONE);
-
-    printf("Calling pmSourceLocalSky with valid data.\n");
-    tmpPeak->x = (psF32) (TST16_NUM_ROWS / 2);
-    tmpPeak->y = (psF32) (TST16_NUM_COLS / 2);
-    mySrc = pmSourceLocalSky(imgData,
-                             tmpPeak,
-                             PS_STAT_SAMPLE_MEAN,
-                             (psF32) TST16_INNER_RADIUS,
-                             (psF32) TST16_OUTER_RADIUS);
-
-    if (mySrc == NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource.\n");
-        testStatus = false;
-    }
-
-    mySrc->modelPSF = pmModelAlloc(PS_MODEL_GAUSS);
-    mySrc->modelPSF->params->data.F32[0] = 5.0;
-    mySrc->modelPSF->params->data.F32[1] = 70.0;
-    mySrc->modelPSF->params->data.F32[2] = (psF32) (TST16_NUM_ROWS / 2);
-    mySrc->modelPSF->params->data.F32[3] = (psF32) (TST16_NUM_COLS / 2);
-    mySrc->modelPSF->params->data.F32[4] = 1.0;
-    mySrc->modelPSF->params->data.F32[5] = 1.0;
-    mySrc->modelPSF->params->data.F32[6] = 2.0;
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceSubModel with NULL psImage.  Should generate error, return FALSE.\n");
-    rc = pmSourceSubModel(NULL, mySrc, true);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceSubModel() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceSubModel with NULL psSrc.  Should generate error, return FALSE.\n");
-    rc = pmSourceSubModel(imgData, NULL, true);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceSubModel() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceSubModel with acceptable data.\n");
-    rc = pmSourceSubModel(imgData, mySrc, true);
-    if (rc != true) {
-        printf("TEST ERROR: pmSourceSubModel() returned FALSE.\n");
-        testStatus = false;
-    }
-
-    psFree(mySrc);
-    psFree(imgData);
-    return(testStatus);
-}
-
-#define TST20_NUM_ROWS 100
-#define TST20_NUM_COLS 100
-#define TST20_SKY 10.0
-#define TST20_INNER_RADIUS 3
-#define TST20_OUTER_RADIUS 5
-/******************************************************************************
-test20(): We first test pmSourceSubModel() with various NULL and unallowable
-input parameters.
- 
-XXX: We don't verify the numbers.
- *****************************************************************************/
-int test20( void )
-{
-    bool testStatus = true;
-    psImage *imgData = psImageAlloc(TST20_NUM_COLS, TST20_NUM_ROWS, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < imgData->numRows; i++) {
-        for (psS32 j = 0 ; j < imgData->numCols; j++) {
-            imgData->data.F32[i][j] = TST20_SKY;
-        }
-    }
-    pmSource *mySrc = NULL;
-    psBool rc = false;
-
-    pmPeak *tmpPeak = pmPeakAlloc((psF32) (TST20_NUM_ROWS / 2),
-                                  (psF32) (TST20_NUM_COLS / 2),
-                                  200.0,
-                                  PM_PEAK_LONE);
-
-    printf("Calling pmSourceLocalSky with valid data.\n");
-    tmpPeak->x = (psF32) (TST20_NUM_ROWS / 2);
-    tmpPeak->y = (psF32) (TST20_NUM_COLS / 2);
-    mySrc = pmSourceLocalSky(imgData,
-                             tmpPeak,
-                             PS_STAT_SAMPLE_MEAN,
-                             (psF32) TST20_INNER_RADIUS,
-                             (psF32) TST20_OUTER_RADIUS);
-
-    if (mySrc == NULL) {
-        printf("TEST ERROR: pmSourceLocalSky() returned a NULL pmSource.\n");
-        testStatus = false;
-    }
-
-    mySrc->modelPSF = pmModelAlloc(PS_MODEL_GAUSS);
-
-
-    mySrc->modelPSF->params->data.F32[0] = 5.0;
-    mySrc->modelPSF->params->data.F32[1] = 70.0;
-    mySrc->modelPSF->params->data.F32[2] = (psF32) (TST20_NUM_ROWS / 2);
-    mySrc->modelPSF->params->data.F32[3] = (psF32) (TST20_NUM_COLS / 2);
-    mySrc->modelPSF->params->data.F32[4] = 1.0;
-    mySrc->modelPSF->params->data.F32[5] = 1.0;
-    mySrc->modelPSF->params->data.F32[6] = 2.0;
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceFitModel with NULL psImage.  Should generate error, return FALSE.\n");
-    rc = pmSourceFitModel(mySrc, NULL);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceFitModel() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceFitModel with NULL pmSource.  Should generate error, return FALSE.\n");
-    rc = pmSourceFitModel(NULL, imgData);
-    if (rc == true) {
-        printf("TEST ERROR: pmSourceFitModel() returned TRUE.\n");
-        testStatus = false;
-    }
-
-    printf("----------------------------------------------------------------------------------\n");
-    printf("Calling pmSourceFitModel with acceptable data.\n");
-    rc = pmSourceFitModel(mySrc, imgData);
-    printf("pmSourceFitModel returned %d\n", rc);
-
-    // XXX: Memory leaks are not being tested
-    psVector *junk = psVectorAlloc(10, PS_TYPE_F32);
-    junk->data.F32[0] = 0.0;
-
-    psFree(mySrc);
-    psFree(imgData);
-    return(testStatus);
-}
-
-
 // this code will