Index: /trunk/psphot/Makefile
===================================================================
--- /trunk/psphot/Makefile	(revision 4976)
+++ /trunk/psphot/Makefile	(revision 4977)
@@ -32,16 +32,27 @@
 $(SRC)/psphotSubtractPSF.$(ARCH).o \
 $(SRC)/psphotSortBySN.$(ARCH).o    \
-$(SRC)/pspsf.$(ARCH).o             \
+$(SRC)/psphotDefinePixels.$(ARCH).o    \
 $(SRC)/psLibUtils.$(ARCH).o	   \
+$(SRC)/psLine.$(ARCH).o		   \
+$(SRC)/psMinimize.$(ARCH).o	   \
+$(SRC)/psPolynomials.$(ARCH).o	   \
+$(SRC)/psEllipse.$(ARCH).o         \
 $(SRC)/psModulesUtils.$(ARCH).o	   \
-$(SRC)/pmSourceUtils.$(ARCH).o	   \
 $(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
-$(SRC)/psLine.$(ARCH).o		   \
-$(SRC)/psPolynomials.$(ARCH).o	   \
-$(SRC)/psMinimize.$(ARCH).o	   \
+$(SRC)/pmPSF.$(ARCH).o             \
+$(SRC)/pmPSFtry.$(ARCH).o          \
 $(SRC)/psImageData.$(ARCH).o	   \
-$(SRC)/psEllipse.$(ARCH).o         \
-$(SRC)/pmModelInit.$(ARCH).o       \
+$(SRC)/pmModelGroup.$(ARCH).o       \
 $(SRC)/pmObjects_EAM.$(ARCH).o
+
+MODELS = \
+$(SRC)/models/pmModel_GAUSS.c      \
+$(SRC)/models/pmModel_PGAUSS.c     \
+$(SRC)/models/pmModel_QGAUSS.c     \
+$(SRC)/models/pmModel_SGAUSS.c      
+
+$(SRC)/pmModelGroup.$(ARCH).o: $(MODELS)
+
+# deprecated
 
 FITSOURCE = \
Index: /trunk/psphot/doc/notes.txt
===================================================================
--- /trunk/psphot/doc/notes.txt	(revision 4976)
+++ /trunk/psphot/doc/notes.txt	(revision 4977)
@@ -47,6 +47,8 @@
   - fixed 'center' option in pmSourceAdd/Sub
 
+  lingering concerns:
+  - Polynomials and related functions are STILL messed up.  this is a
+    major priority!!!
   - pmSourceFitModel calculates sqrt(var), psMinimize calculates sq(sqrt(var))
-  - pmSourceContour is using the model flux, not the image flux
   - does psImageAlloc zero the image?  if so, drop init routine
     from pmSourcePSFClump
@@ -60,5 +62,5 @@
     the psImageMaskRegion, etc functions.
   - pmSourceSetPixelsCircle uses a hard-coded mask value
-  - 
+  - pmSourceContour is using the model flux, not the image flux
 
 2005.09.05
Index: /trunk/psphot/doc/psphot.txt
===================================================================
--- /trunk/psphot/doc/psphot.txt	(revision 4976)
+++ /trunk/psphot/doc/psphot.txt	(revision 4977)
@@ -1,31 +1,112 @@
+
+TBD priorities:
+
+- test psImageSmooth / remove psImageSmooth_EAM
+- unify pmSourceLocalSky / pmSourceLocalSky_EAM
+- unify pmSourceMoments_EAM / pmSourceMoments_EAM
+- merge with image data hierarchy stuff
+- better name for psPSF_Test
 
 Defined APIs:
 
-pmMoments	    *pmMomentsAlloc();
-pmModel  	    *pmModelAlloc();
-pmSource 	    *pmSourceAlloc();
-psVector 	    *pmFindVectorPeaks()
-psArray  	    *pmFindImagePeaks()
-psList   	    *pmCullPeaks()
-pmSource 	    *pmSourceLocalSky()
-bool 		     pmSourceMoments()
-pmPSFClump 	     pmSourcePSFClump()
-bool 		     pmSourceRoughClass()
-bool 		     pmSourceSetPixelsCircle()
-pmModel 	    *pmSourceModelGuess()
-psArray 	    *pmSourceContour()
-bool 		     pmSourceFitModel()
-bool 		     pmSourceAddModel()
-bool 		     pmSourceSubModel()
-int                  pmModelParameterCount ();
-char                *pmModelGetType ();
-pmModelType          pmModelSetType ();
-pmModelFunc          pmModelFunc_GetFunction ();
-pmModelFlux          pmModelFlux_GetFunction ();
-pmModelRadius        pmModelRadius_GetFunction ();
-pmModelLimits        pmModelLimits_GetFunction ();
-pmModelGuessFunc     pmModelGuessFunc_GetFunction ();
-pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction ();
-pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction ();
+top-level psphot functions : these do not need to be in psModules
+  psMetadata  *psphotArguments ();
+  psImageData *psphotSetup ();
+  psStats     *psphotImageStats ();
+  psArray     *pmPeaksSigmaLimit ();
+  psArray     *psphotSourceStats ();
+  pmPSF       *psphotChoosePSF ();
+  bool         psphotApplyPSF ();
+  bool         psphotFitGalaxies ();
+  void         psphotOutput ();
+
+psphot-specific functions
+  bool         psphotMarkPSF ();
+  bool         psphotSubtractPSF ();
+  int          psphotSortBySN ();
+  int          psphotSaveImage ();
+
+** add the following entries to the psLibSDRS
+
+psLib extra utilities
+  bool         psTimerStart ();		- already in SDRS
+  void         psTimerFree (); 		- already in SDRS
+  bool         psTimerClear ();		- add to psLibSDRS
+  psF64        psTimerMark (); 		- already in SDRS
+  psS32        psLogArguments ();	- now: psArgumentVerbosity in SDRS
+  psS32        psTraceArguments ();	- now: psArgumentVerbosity in SDRS
+  int          psArgumentGet ();	- already in SDRS
+  int          psArgumentRemove ();	- already in SDRS
+  psVector    *psVectorCreate ();	- add to psLibSDRS
+  int          psImageCountPixelMask ();- add to psLibSDRS
+  psVector    *psGetRowVectorFromImage(); - add to psLibSDRS?
+  bool         p_psVectorPrintRow ();     - add to psLibSDRS?
+
+// basic image functions
+  bool         psImageInit ();		- already in SDRS
+  void         psImageSmooth_EAM ();    - test psLib version and drop
+
+// psLine functions
+  psLine      *psLineAlloc ();		- add to psLibSDRS
+  bool         psLineInit (); 		- add to psLibSDRS
+  bool         psLineAdd ();  		- add to psLibSDRS
+
+** compare with current implementation / update psLib
+   we need to define a mechanism to carry the beta values and the
+   parameter limits
+
+// minimize 
+  psBool       p_psMinLM_GuessABP_EAM ();
+  psBool       psMinimizeLMChi2_EAM();
+  psF64        p_psMinLM_dLinear ();  -- not included in pslib.h files
+
+** already in psLib : some need to be cleaned up by MHPCC
+
+// polynomial functions
+  psF64           Polynomial2DEval();
+  psImage        *psBuildSums2D();
+  psPolynomial2D *VectorFitPolynomial2DOrd_EAM();
+  psPolynomial2D *Polynomial2DAlloc();
+  void            psPolynomial2DDump ();
+  psPolynomial2D *RobustFit2D_nomask();
+  psPolynomial2D *RobustFit2D();
+  psVector       *Polynomial2DEvalVector();
+
+  psVector       *psBuildSums1D();
+  void            psPolynomial1DDump ();
+  psF64           Polynomial1DEval_EAM();
+  psVector       *Polynomial1DEvalVector_EAM();
+  psPolynomial1D *Polynomial1DAlloc();
+  psPolynomial1D *VectorFitPolynomial1DOrd_EAM();
+
+** make sure these are all in the psModules SDRS:
+
+pmObjects.h:
+  pmMoments           *pmMomentsAlloc();
+  pmModel             *pmModelAlloc();
+  pmSource            *pmSourceAlloc();
+  psVector            *pmFindVectorPeaks()
+  psArray             *pmFindImagePeaks()
+  psList              *pmCullPeaks()
+  pmSource            *pmSourceLocalSky()
+  bool                 pmSourceMoments()
+  pmPSFClump           pmSourcePSFClump()
+  bool                 pmSourceRoughClass()
+  bool                 pmSourceSetPixelsCircle()
+  pmModel             *pmSourceModelGuess()
+  psArray             *pmSourceContour()
+  bool                 pmSourceFitModel()
+  bool                 pmSourceAddModel()
+  bool                 pmSourceSubModel()
+  int                  pmModelParameterCount ();
+  char                *pmModelGetType ();
+  pmModelType          pmModelSetType ();
+  pmModelFunc          pmModelFunc_GetFunction ();
+  pmModelFlux          pmModelFlux_GetFunction ();
+  pmModelRadius        pmModelRadius_GetFunction ();
+  pmModelLimits        pmModelLimits_GetFunction ();
+  pmModelGuessFunc     pmModelGuessFunc_GetFunction ();
+  pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction ();
+  pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction ();
 
 model function abstractions:
@@ -37,2 +118,43 @@
   typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
   typedef bool (*pmModelFitStatusFunc)(pmModel *model);
+
+conversions between elliptical shape representations
+  EllipseAxes  EllipseMomentsToAxes ();
+  EllipseShape EllipseAxesToShape ();
+  EllipseAxes  EllipseShapeToAxes ();
+
+output functions
+  bool         pmSourcesWriteText ();  	-- be careful about psImageData...
+  bool         pmSourcesWriteOBJ ();   	-- be careful about psImageData...
+  bool         pmSourcesWriteCMP ();   	-- be careful about psImageData...
+  bool         pmSourcesWriteCMF ();   	-- be careful about psImageData...
+  bool         pmSourcesWriteSX ();    	-- be careful about psImageData...
+  int          pmSourcesDophotType ();
+  bool         pmPeaksWriteText ();
+  bool         pmMomentsWriteText ();
+  bool         pmModelWritePSFs ();
+  bool         pmModelWriteFLTs ();
+  bool         pmModelWriteNULLs ();
+
+// psf utilities
+  pmPSF       *pmPSFAlloc ();
+  pmPSF_Test  *pmPSF_TestAlloc ();
+  pmPSF_Test  *pmPSF_TestModel ();
+  bool         pmPSFFromModels ();
+  pmModel     *pmModelFromPSF ();
+  bool         pmSourcePhotometry ();
+  bool         pmPSFMetricModel ();
+
+// psModule extra utilities
+  bool         pmSourceFitModel_EAM();  -- requires updated version of psMinimize...
+  bool         pmSourceLocalSky_EAM (); -- compare with pmSourceLocalSky & choose
+  bool         pmSourceMoments_EAM();   -- compare with pmSourceMoments & choose
+  bool         pmSourceDefinePixels();  -- be careful about psImageData...
+  bool         pmModelFitStatus ();
+  int          pmSourceDophotType ();
+  psF32        pmConfigLookupF32 ();	-- merge with Paul's image container hierarchy
+
+** this needs to be confronted with the phase2 / image container hierarchy
+
+// psImageData functions
+  psImageData *psImageDataAlloc ();
Index: /trunk/psphot/src/models/pmModel_GAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_GAUSS.c	(revision 4976)
+++ /trunk/psphot/src/models/pmModel_GAUSS.c	(revision 4977)
@@ -124,5 +124,5 @@
     for (int i = 4; i < 7; i++) {
 	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = Polynomial2DEval (poly, out[2], out[3]);
+	out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]);
     }
     return(true);
Index: /trunk/psphot/src/models/pmModel_PGAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_PGAUSS.c	(revision 4976)
+++ /trunk/psphot/src/models/pmModel_PGAUSS.c	(revision 4977)
@@ -140,5 +140,5 @@
     for (int i = 4; i < 7; i++) {
 	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = Polynomial2DEval (poly, out[2], out[3]);
+	out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]);
     }
     return(true);
Index: /trunk/psphot/src/models/pmModel_QGAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_QGAUSS.c	(revision 4976)
+++ /trunk/psphot/src/models/pmModel_QGAUSS.c	(revision 4977)
@@ -165,5 +165,5 @@
     for (int i = 4; i < 8; i++) {
 	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = Polynomial2DEval (poly, out[2], out[3]);
+	out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]);
     }
     return(true);
Index: /trunk/psphot/src/models/pmModel_SGAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_SGAUSS.c	(revision 4976)
+++ /trunk/psphot/src/models/pmModel_SGAUSS.c	(revision 4977)
@@ -300,5 +300,5 @@
     for (int i = 4; i < 9; i++) {
 	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = Polynomial2DEval (poly, out[2], out[3]);
+	out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]);
     }
     return(true);
Index: /trunk/psphot/src/pmModelGroup.c
===================================================================
--- /trunk/psphot/src/pmModelGroup.c	(revision 4977)
+++ /trunk/psphot/src/pmModelGroup.c	(revision 4977)
@@ -0,0 +1,109 @@
+# include "pmModelGroup.h"
+
+double hypot(double x, double y);
+double sqrt (double x);
+
+# 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/psphot/src/pmModelGroup.h
===================================================================
--- /trunk/psphot/src/pmModelGroup.h	(revision 4977)
+++ /trunk/psphot/src/pmModelGroup.h	(revision 4977)
@@ -0,0 +1,104 @@
+
+# ifndef PM_MODEL_GROUP_H 
+# define PM_MODEL_GROUP_H 
+
+# include "pmObjects_EAM.h"
+# include "pmPSF.h"
+
+/**
+ * 
+ *  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.
+ * 
+ */
+
+/**
+ * 
+ *  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          pmModelFunc_GetFunction (pmModelType type);
+pmModelFlux          pmModelFlux_GetFunction (pmModelType type);
+pmModelRadius        pmModelRadius_GetFunction (pmModelType type);
+pmModelLimits        pmModelLimits_GetFunction (pmModelType type);
+pmModelGuessFunc     pmModelGuessFunc_GetFunction (pmModelType type);
+pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction (pmModelType type);
+pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type);
+
+// pmModelGroup utility functions
+int                  pmModelParameterCount (pmModelType type);
+char                *pmModelGetType (pmModelType type);
+pmModelType          pmModelSetType (char *name);
+
+// structure to carry model group functions
+typedef struct {
+    char *name;
+    int nParams;
+    pmModelFunc          modelFunc;
+    pmModelFlux          modelFlux;
+    pmModelRadius        modelRadius;
+    pmModelLimits        modelLimits;
+    pmModelGuessFunc     modelGuessFunc;
+    pmModelFromPSFFunc   modelFromPSFFunc;
+    pmModelFitStatusFunc modelFitStatusFunc;
+} pmModelGroup;
+
+# endif
Index: unk/psphot/src/pmModelInit.c
===================================================================
--- /trunk/psphot/src/pmModelInit.c	(revision 4976)
+++ 	(revision )
@@ -1,111 +1,0 @@
-# include <string.h>
-# include <math.h>
-# include "pmObjects_EAM.h"
-
-double hypot(double x, double y);
-double sqrt (double x);
-
-# 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/psphot/src/pmObjects_EAM.c
===================================================================
--- /trunk/psphot/src/pmObjects_EAM.c	(revision 4976)
+++ /trunk/psphot/src/pmObjects_EAM.c	(revision 4977)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-07 03:33:01 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-08 04:32:40 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -19,4 +19,6 @@
 #include "psConstants.h"
 #include "pmObjects_EAM.h"
+#include "pmModelGroup.h"
+
 
 /******************************************************************************
@@ -653,10 +655,9 @@
 
     //
-    // Loop through the subimage mask, initialize mask to 0x01 (invalid pixel)
-    // XXX EAM : use PSPHOT_MASK_INVALID?
+    // Loop through the subimage mask & initialize PSPHOT_MASK_INVALID
     //
     for (psS32 row = 0 ; row < subImageMask->numRows; row++) {
         for (psS32 col = 0 ; col < subImageMask->numCols; col++) {
-            subImageMask->data.U8[row][col] = 0x01;
+            subImageMask->data.U8[row][col] = PSPHOT_MASK_INVALID;
         }
     }
@@ -674,5 +675,5 @@
     for (psS32 row = StartRow; row < EndRow; row++) {
         for (psS32 col = StartCol; col < EndCol; col++) {
-	    subImageMask->data.U8[row][col] = 0;
+	    subImageMask->data.U8[row][col] = PSPHOT_MASK_CLEAR;
         }
     }
@@ -1116,6 +1117,5 @@
 
 	// XXX EAM : can we use the value of SATURATE if mask is NULL?
-# define MASK_SATURATED 0x02
-	int Nsatpix = psImageCountPixelMask (tmpSrc->mask, MASK_SATURATED);
+	int Nsatpix = psImageCountPixelMask (tmpSrc->mask, PSPHOT_MASK_SATURATED);
 
         // saturated star (size consistent with PSF or larger)
@@ -1615,5 +1615,5 @@
     }       
 
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
 
     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
@@ -1621,5 +1621,5 @@
     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]);
+	dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
     }
  
@@ -1640,9 +1640,9 @@
     // XXX EAM get the Gauss-Newton distance for fixed model parameters
     if (paramMask != NULL) {
-	psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
+	psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
 	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];
+	    dparams->data.F32[i] = delta->data.F32[i];
 	}
 	psFree (delta);
Index: /trunk/psphot/src/pmObjects_EAM.h
===================================================================
--- /trunk/psphot/src/pmObjects_EAM.h	(revision 4976)
+++ /trunk/psphot/src/pmObjects_EAM.h	(revision 4977)
@@ -5,6 +5,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-07 03:33:01 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-08 04:32:40 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -22,22 +22,18 @@
 #include<math.h>
 #include "pslib.h"
-
-// XXX EAM : this function should be added to psLib
-int          psImageCountPixelMask (psImage *mask, psU8 value);
-
-// XXX EAM : bug : missing from psLib/*.h
-bool psMinimizeGaussNewtonDelta (psVector *delta,
-                                 const psVector *params,
-                                 const psVector *paramMask,
-                                 const psArray  *x,
-                                 const psVector *y,
-                                 const psVector *yErr,
-                                 psMinimizeLMChi2Func func);
-
-// XXX EAM : switch to psPoly when ready
-psF64	        Polynomial2DEval(const psPolynomial2D* myPoly, psF64 x, psF64 y);
-
-// XXX EAM : psEllipse needs to be be included in psLib
-# include "psEllipse.h"
+#include "psLibUtils.h"
+// #include "pmModelGroup.h"
+
+// defined mask values for mask image pixels.
+// 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
@@ -91,6 +87,4 @@
 } pmMoments;
 
-typedef int pmModelType;
-
 /** pmPSFClump data structure
  *  
@@ -105,4 +99,6 @@
     float dY;
 } pmPSFClump;
+
+typedef int pmModelType;
 
 /** pmModel data structure
@@ -169,18 +165,4 @@
 pmSource;
 
-/** pmPSF data structure
- *  
- *  
- *  
- */
-typedef struct
-{
-    pmModelType type;			///< PSF Model in use
-    psArray *params;			///< Model parameters (psPolynomial2D)
-    float chisq;			///< PSF goodness statistic
-    int nPSFstars;			///< number of stars used to measure PSF
-}
-pmPSF;
-
 pmPeak *pmPeakAlloc(
     int x,				///< Row-coordinate in image space
@@ -316,100 +298,3 @@
 );
 
-/**
- * 
- *  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.
- * 
- */
-
-/**
- * 
- *  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          pmModelFunc_GetFunction (pmModelType type);
-pmModelFlux          pmModelFlux_GetFunction (pmModelType type);
-pmModelRadius        pmModelRadius_GetFunction (pmModelType type);
-pmModelLimits        pmModelLimits_GetFunction (pmModelType type);
-pmModelGuessFunc     pmModelGuessFunc_GetFunction (pmModelType type);
-pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction (pmModelType type);
-pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type);
-
-// pmModelGroup utility functions
-int                  pmModelParameterCount (pmModelType type);
-char                *pmModelGetType (pmModelType type);
-pmModelType          pmModelSetType (char *name);
-
-// structure to carry model group functions
-typedef struct {
-    char *name;
-    int nParams;
-    pmModelFunc          modelFunc;
-    pmModelFlux          modelFlux;
-    pmModelRadius        modelRadius;
-    pmModelLimits        modelLimits;
-    pmModelGuessFunc     modelGuessFunc;
-    pmModelFromPSFFunc   modelFromPSFFunc;
-    pmModelFitStatusFunc modelFitStatusFunc;
-} pmModelGroup;
-
 #endif
Index: /trunk/psphot/src/pmPSF.c
===================================================================
--- /trunk/psphot/src/pmPSF.c	(revision 4977)
+++ /trunk/psphot/src/pmPSF.c	(revision 4977)
@@ -0,0 +1,121 @@
+# include <pslib.h>
+# include "psLibUtils.h"
+# include "pmObjects_EAM.h"
+# include "pmPSF.h"
+# include "pmModelGroup.h"
+
+/********  pmPSF functions  **************************************************/
+
+// function to free a pmPSF structure
+static void pmPSFFree (pmPSF *psf) {
+
+  if (psf == NULL) return;
+
+  psFree (psf->params);
+  return;
+}
+
+// allocate a pmPSF
+// a PSF always has 4 parameters fewer than the equivalent model
+pmPSF *pmPSFAlloc (pmModelType type) {
+
+    int Nparams;
+
+    pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF));
+
+    psf->type = type;
+    psf->chisq = 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] = Polynomial2DAlloc_EAM(PS_POLYNOMIAL_ORD, 1, 1);
+    }
+
+    psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
+    return(psf);
+}
+
+// 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
+// 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] = VectorClipFitPolynomial2D_EAM (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);
+}
+
+// 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/psphot/src/pmPSF.h
===================================================================
--- /trunk/psphot/src/pmPSF.h	(revision 4977)
+++ /trunk/psphot/src/pmPSF.h	(revision 4977)
@@ -0,0 +1,24 @@
+
+# ifndef PM_PSF_H 
+# define PM_PSF_H
+
+/** pmPSF data structure
+ *  
+ *  
+ *  
+ */
+typedef struct
+{
+    pmModelType type;			///< PSF Model in use
+    psArray *params;			///< Model parameters (psPolynomial2D)
+    float chisq;			///< PSF goodness statistic
+    int nPSFstars;			///< number of stars used to measure PSF
+}
+pmPSF;
+
+// pmPSF utilities
+pmPSF       *pmPSFAlloc (pmModelType type);
+bool	     pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
+pmModel	    *pmModelFromPSF (pmModel *model, pmPSF *psf);
+
+# endif
Index: /trunk/psphot/src/pmPSFtry.c
===================================================================
--- /trunk/psphot/src/pmPSFtry.c	(revision 4977)
+++ /trunk/psphot/src/pmPSFtry.c	(revision 4977)
@@ -0,0 +1,295 @@
+# include <pslib.h>
+# include "psLibUtils.h"
+# include "pmObjects_EAM.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) {
+
+    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
+
+    // XXX probably need to increment ref counter
+    test->modelType   = pmModelSetType (modelName);
+    test->psf         = pmPSFAlloc (test->modelType);
+    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);
+    test->ApResid     = 0;
+    test->dApResid    = 0;
+    test->skyBias     = 0;
+
+    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->modelType); 
+	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->ApResid, try->dApResid, try->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 (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 = Polynomial1DAlloc_EAM (PS_POLYNOMIAL_ORD, 1);
+
+  // XXX EAM : this is the intended API (cycle 7? cycle 8?) 
+  poly = VectorFitPolynomial1D_EAM (poly, maskB, 1, daBin, NULL, rfBin);
+
+  // XXX EAM : replace this when the above version is implemented
+  // poly = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL);
+
+  psVector *daBinFit = Polynomial1DEvalVector_EAM (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->ApResid = poly->coeff[0];
+  try->dApResid = stats->clippedStdev;
+  try->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/psphot/src/pmPSFtry.h
===================================================================
--- /trunk/psphot/src/pmPSFtry.h	(revision 4977)
+++ /trunk/psphot/src/pmPSFtry.h	(revision 4977)
@@ -0,0 +1,34 @@
+
+# ifndef PM_PSF_TRY_H 
+# define PM_PSF_TRY_H
+
+// data to try a given PSF model type
+typedef struct {
+    pmModelType modelType;
+    pmPSF      *psf;
+    psArray    *sources;      // pointers to the original sources
+    psArray    *modelFLT;     // model fits, floating parameters 
+    psArray    *modelPSF;     // model fits, PSF parameters
+    psVector   *mask;
+    psVector   *metric;
+    psVector   *fitMag;
+    float       ApResid;
+    float       dApResid;
+    float       skyBias;
+} pmPSFtry;
+
+enum {
+    PSFTRY_MASK_CLEAR    = 0x00,
+    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,
+} pmPSFtryMaskValues;
+
+// pmPSFtry utilities
+pmPSFtry    *pmPSFtryAlloc (psArray *stars, char *modelName);
+pmPSFtry    *pmPSFtryModel (psArray *sources, char *modelName, float radius);
+bool	     pmPSFtryMetric (pmPSFtry *try, float RADIUS);
+
+# endif
Index: /trunk/psphot/src/pmSourceUtils.c
===================================================================
--- /trunk/psphot/src/pmSourceUtils.c	(revision 4976)
+++ /trunk/psphot/src/pmSourceUtils.c	(revision 4977)
@@ -1,44 +1,1 @@
 # include "psphot.h"
-
-bool pmSourceDefinePixels(pmSource *mySource, 
-			  const psImageData *imdata,
-			  psF32 x, 
-			  psF32 y,
-			  psF32 Radius)
-{
-    psRegion srcRegion;
-
-    // Grab a subimage of the original image of size (2 * outerRadius).
-    srcRegion = psRegionForSquare (x, y, Radius);
-    srcRegion = psRegionForImage (imdata->image, &srcRegion);
-
-    mySource->pixels = psImageSubset(imdata->image, srcRegion);
-    mySource->noise  = psImageSubset(imdata->noise, srcRegion);
-    mySource->mask   = psImageSubset(imdata->mask,  srcRegion);
-
-    return(mySource);
-}
-
-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/psphot/src/psLibUtils.c
===================================================================
--- /trunk/psphot/src/psLibUtils.c	(revision 4976)
+++ /trunk/psphot/src/psLibUtils.c	(revision 4977)
@@ -1,4 +1,5 @@
 # include <strings.h>  // for strncasecmp
 # include <pslib.h>
+# include "psLibUtils.h"
 
 // XXX EAM : this is NOT included in the C99 headers ??
@@ -157,51 +158,4 @@
     memcpy (tmpVector->data.F32, image->data.F32[row], image->numCols*sizeof(psF32));
     return(tmpVector);
-}
-
-// extract config informatin from config data or from image header, if specified
-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);
 }
 
Index: /trunk/psphot/src/psLibUtils.h
===================================================================
--- /trunk/psphot/src/psLibUtils.h	(revision 4977)
+++ /trunk/psphot/src/psLibUtils.h	(revision 4977)
@@ -0,0 +1,74 @@
+
+# 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 : bug : missing from psLib/*.h
+bool psMinimizeGaussNewtonDelta (psVector *delta,
+                                 const psVector *params,
+                                 const psVector *paramMask,
+                                 const psArray  *x,
+                                 const psVector *y,
+                                 const psVector *yErr,
+                                 psMinimizeLMChi2Func func);
+
+// minimize 
+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);
+void 	     psTimerFree ();
+bool 	     psTimerClear (char *name);
+psF64 	     psTimerMark (char *name);
+psS32 	     psLogArguments (int *argc, char **argv);
+psS32 	     psTraceArguments (int *argc, char **argv);
+int	     psArgumentGet (int argc, char **argv, char *arg);
+int	     psArgumentRemove (int N, int *argc, char **argv);
+psVector    *psVectorCreate (double lower, double upper, double delta, psElemType type);
+psVector    *psGetRowVectorFromImage(psImage *image, psU32 row);
+int          psImageCountPixelMask (psImage *mask, psU8 value);
+
+// basic image functions
+bool         psImageInit (psImage *image,...);
+void	     psImageSmooth_EAM (psImage *image, float sigma, float Nsigma);
+
+// psLine functions
+psLine      *psLineAlloc (int Nline);
+bool	     psLineInit (psLine *line);
+bool	     psLineAdd (psLine *line, char *format, ...);
+
+// not included in the .h file?
+bool p_psVectorPrint (int fd, psVector *a, char *name);
+
+// not included in the .h file?
+bool p_psVectorPrintRow (int fd, psVector *a, char *name);
+
+// polynomial functions
+psPolynomial1D *Polynomial1DAlloc_EAM(psPolynomialType type, psS32 nOrder);
+void	        Polynomial1DDump_EAM(psPolynomial1D *poly);
+psF64		Polynomial1DEval_EAM(const psPolynomial1D* myPoly, psF64 x);
+psVector       *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x);
+psPolynomial1D *VectorFitPolynomial1D_EAM(psPolynomial1D* myPoly, psVector* mask, psMaskType maskValue, const psVector* y, const psVector* yErr, const psVector* x);
+
+psPolynomial2D *Polynomial2DAlloc_EAM(psPolynomialType type, psS32 nXorder, psS32 nYorder);
+void	        Polynomial2DDump_EAM(psPolynomial2D *poly);
+psF64	        Polynomial2DEval_EAM(const psPolynomial2D* myPoly, psF64 x, psF64 y);
+psVector       *Polynomial2DEvalVector_EAM(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
+psVector       *Polynomial2DEvalVectorD_EAM(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
+psPolynomial2D *VectorFitPolynomial2D_EAM(psPolynomial2D* myPoly, psVector* mask, psMaskType maskValue, const psVector* z, const psVector* zErr, const psVector* x, const psVector* y);
+psPolynomial2D *VectorClipFitPolynomial2D_EAM(psPolynomial2D* poly, psStats *stats, psVector* mask, psMaskType maskValue, const psVector* z, const psVector* zErr, const psVector* x, const psVector* y);
+
+# endif
Index: /trunk/psphot/src/psLine.c
===================================================================
--- /trunk/psphot/src/psLine.c	(revision 4976)
+++ /trunk/psphot/src/psLine.c	(revision 4977)
@@ -1,3 +1,4 @@
-# include "psphot.h"
+# include <pslib.h>
+# include "psLibUtils.h"
 
 static void psLineFree (psLine *line) {
Index: /trunk/psphot/src/psMinimize.c
===================================================================
--- /trunk/psphot/src/psMinimize.c	(revision 4976)
+++ /trunk/psphot/src/psMinimize.c	(revision 4977)
@@ -1,3 +1,4 @@
-# include "psphot.h"
+# include <pslib.h>
+# include "psLibUtils.h"
 
 // XXX EAM : this file contains my alternate implementations of psMinimizeLMChi2
Index: /trunk/psphot/src/psModulesUtils.c
===================================================================
--- /trunk/psphot/src/psModulesUtils.c	(revision 4976)
+++ /trunk/psphot/src/psModulesUtils.c	(revision 4977)
@@ -1,3 +1,50 @@
-# include "psphot.h"
+# include "psModulesUtils.h"
+
+// extract config informatin from config data or from image header, if specified
+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);
+}
 
 // XXX EAM : these are my alternate implementations of psModule functions
@@ -321,2 +368,26 @@
     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/psphot/src/psModulesUtils.h
===================================================================
--- /trunk/psphot/src/psModulesUtils.h	(revision 4977)
+++ /trunk/psphot/src/psModulesUtils.h	(revision 4977)
@@ -0,0 +1,21 @@
+
+# ifndef PS_MODULE_UTILS
+# define PS_MODULE_UTILS
+
+# include <strings.h>  // for strcasecmp
+# include <pslib.h>
+# include "psLibUtils.h"
+# include "pmObjects_EAM.h"
+# include "pmModelGroup.h"
+
+// psModule extra utilities
+bool 	     pmSourceLocalSky_EAM (pmSource *source, psStatsOptions statsOptions, psF32 Radius);
+bool 	     pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF);
+bool 	     pmSourceMoments_EAM(pmSource *source, psF32 radius);
+bool 	     pmModelFitStatus (pmModel *model);
+int	     pmSourceDophotType (pmSource *source);
+bool	     pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask);
+
+psF32        pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
+
+# endif
Index: /trunk/psphot/src/psPolynomials.c
===================================================================
--- /trunk/psphot/src/psPolynomials.c	(revision 4976)
+++ /trunk/psphot/src/psPolynomials.c	(revision 4977)
@@ -1,3 +1,4 @@
-# include "psphot.h"
+# include <pslib.h>
+# include "psLibUtils.h"
 
 // XXX EAM : this file defines alternate versions of the polynomial fitting
@@ -9,5 +10,5 @@
 
 // write out the terms of the given 1D polynomial
-void psPolynomial1DDump (psPolynomial1D *poly) {
+void Polynomial1DDump_EAM (psPolynomial1D *poly) {
 
     for (int i = 0; i < poly->n + 1; i++) {
@@ -16,50 +17,17 @@
 }    
 
-psF64 Polynomial1DEval_EAM(psF64 x, const psPolynomial1D* myPoly)
-{
-    psS32 loop_x = 0;
-    psF64 polySum = 0.0;
-    psF64 xSum = 1.0;
-
-    for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
-        if (myPoly->mask[loop_x] == 0) {
-            polySum += xSum * myPoly->coeff[loop_x];
-        }
-        xSum *= x;
-    }
-
-    return(polySum);
-}
-
-psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly,
-				       const psVector *x)
-{
-    // PS_POLY_CHECK_NULL(myPoly, NULL);
-    // PS_VECTOR_CHECK_NULL(x, NULL);
-    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
-
-    psVector *tmp;
-
-    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
-    for (psS32 i=0;i<x->n;i++) {
-        tmp->data.F64[i] = Polynomial1DEval_EAM(x->data.F64[i], myPoly);
-    }
-
-    return(tmp);
-}
-
-static void psPolynomial1DFree (psPolynomial1D *poly) {
-
-  if (poly == NULL) return;
-  psFree (poly->coeff);
-  psFree (poly->coeffErr);
-  psFree (poly->mask);
-  return;
+static void Polynomial1DFree_EAM (psPolynomial1D *poly) {
+
+    if (poly == NULL) return;
+    psFree (poly->coeff);
+    psFree (poly->coeffErr);
+    psFree (poly->mask);
+    return;
 }
 
 // XXX EAM : use Nterm = Norder + 1 definition  
 // XXX EAM : should we provide both order and nterms in struct?
-psPolynomial1D* Polynomial1DAlloc(psS32 nOrder,
-				  psPolynomialType type)
+psPolynomial1D* Polynomial1DAlloc_EAM(psPolynomialType type,
+				      psS32 nOrder)
 {
     // PS_INT_CHECK_NON_NEGATIVE(nOrder, NULL);
@@ -80,12 +48,45 @@
         newPoly->mask[i] = 0;
     }
-    psMemSetDeallocator(newPoly, (psFreeFunc) psPolynomial1DFree);
+    psMemSetDeallocator(newPoly, (psFreeFunc) Polynomial1DFree_EAM);
     return(newPoly);
 }
 
+psF64 Polynomial1DEval_EAM(const psPolynomial1D* myPoly, psF64 x)
+{
+    psS32 loop_x = 0;
+    psF64 polySum = 0.0;
+    psF64 xSum = 1.0;
+
+    for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
+        if (myPoly->mask[loop_x] == 0) {
+            polySum += xSum * myPoly->coeff[loop_x];
+        }
+        xSum *= x;
+    }
+
+    return(polySum);
+}
+
+psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly,
+				     const psVector *x)
+{
+    // PS_POLY_CHECK_NULL(myPoly, NULL);
+    // PS_VECTOR_CHECK_NULL(x, NULL);
+    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+
+    psVector *tmp;
+
+    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
+    for (psS32 i=0;i<x->n;i++) {
+        tmp->data.F64[i] = Polynomial1DEval_EAM(myPoly, x->data.F64[i]);
+    }
+
+    return(tmp);
+}
+
 // XXX EAM : my alternate BuildSums1D
-static psVector *BuildSums1D(psVector* sums, 
-			     psF64 x,
-			     psS32 nTerm)
+static psVector *BuildSums1D_EAM(psVector* sums, 
+				 psF64 x,
+				 psS32 nTerm)
 {
     psS32 nSum = 0;
@@ -109,9 +110,10 @@
 
 // XXX EAM : test version of 1d fitting
-psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly,
-					     psVector *mask,
-					     const psVector *x,
-					     const psVector *y,
-					     const psVector *yErr)
+psPolynomial1D* VectorFitPolynomial1D_EAM(psPolynomial1D* myPoly,
+					  psVector *mask,
+					  psMaskType maskValue,
+					  const psVector *y,
+					  const psVector *yErr,
+					  const psVector *x)
 {
     // I think this is 1 dimension down
@@ -130,9 +132,9 @@
     // XXX EAM : change from FILE to fd breaks this code:
     if (psTraceGetLevel (".psLib.dataManip.VectorFitPolynomial1DOrd") >= 5) {
-      FILE *f = psTraceGetDestination ();
-      fprintf (f, "VectorFitPolynomial1D()\n");
-      for (int i = 0; i < x->n; i++) {
-	fprintf (f, "(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
-      }
+	FILE *f = psTraceGetDestination ();
+	fprintf (f, "VectorFitPolynomial1D()\n");
+	for (int i = 0; i < x->n; i++) {
+	    fprintf (f, "(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
+	}
     }
 # endif
@@ -155,6 +157,6 @@
     // Build the B and A data structs.
     for (int k = 0; k < x->n; k++) {
-        if ((mask != NULL) && mask->data.U8[k]) continue;
-	xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
+        if ((mask != NULL) && (mask->data.U8[k] & maskValue)) continue;
+	xSums = BuildSums1D_EAM(xSums, x->data.F64[k], nTerm);
       
 	if (yErr == NULL) {
@@ -190,5 +192,5 @@
     } 
     else 
-    // LUD version of the fit
+	// LUD version of the fit
     {
 	psImage *ALUD = NULL;
@@ -218,86 +220,23 @@
 // ********************** 2D polynomial functions ******************
 
-// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
-psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly,
-				 const psVector *x,
-				 const psVector *y)
-
-{
-    // PS_POLY_CHECK_NULL(myPoly, NULL);
-    // PS_VECTOR_CHECK_NULL(x, NULL);
-    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
-    // PS_VECTOR_CHECK_NULL(y, NULL);
-    // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
-
-    psVector *tmp;
-    psS32 vecLen=x->n;
-
-    // Determine the length of the output vector to by the minimum of the x,y vectors
-    if (y->n < vecLen) {
-        vecLen = y->n;
-    }
-
-    // Create output vector to return
-    tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
-
-    // Evaluate the polynomial at the specified points
-    for (psS32 i=0; i<vecLen; i++) {
-        tmp->data.F32[i] = Polynomial2DEval(myPoly,x->data.F32[i],y->data.F32[i]);
-    }
-
-    // Return output vector
-    return(tmp);
-}
-
-// XXX EAM : this version uses the F64 vectors
-psVector *Polynomial2DEvalVectorD(const psPolynomial2D *myPoly,
-				  const psVector *x,
-				  const psVector *y)
-
-{
-    // PS_POLY_CHECK_NULL(myPoly, NULL);
-    // PS_VECTOR_CHECK_NULL(x, NULL);
-    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
-    // PS_VECTOR_CHECK_NULL(y, NULL);
-    // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
-
-    psVector *tmp;
-    psS32 vecLen=x->n;
-
-    // Determine the length of the output vector to by the minimum of the x,y vectors
-    if (y->n < vecLen) {
-        vecLen = y->n;
-    }
-
-    // Create output vector to return
-    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
-
-    // Evaluate the polynomial at the specified points
-    for (psS32 i=0; i<vecLen; i++) {
-        tmp->data.F64[i] = Polynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]);
-    }
-
-    // Return output vector
-    return(tmp);
-}
-
-static void psPolynomial2DFree(psPolynomial2D *poly) {
-
-  if (poly == NULL) return;
-  for (int i = 0; i < poly->nX + 1; i++) {
-    psFree (poly->coeff[i]);
-    psFree (poly->coeffErr[i]);
-    psFree (poly->mask[i]);
-  }
-  psFree (poly->coeff);
-  psFree (poly->coeffErr);
-  psFree (poly->mask);
-  return;
+static void Polynomial2DFree_EAM(psPolynomial2D *poly) {
+
+    if (poly == NULL) return;
+    for (int i = 0; i < poly->nX + 1; i++) {
+	psFree (poly->coeff[i]);
+	psFree (poly->coeffErr[i]);
+	psFree (poly->mask[i]);
+    }
+    psFree (poly->coeff);
+    psFree (poly->coeffErr);
+    psFree (poly->mask);
+    return;
 }
 
 // XXX EAM : use Nterm = Norder + 1 definition  
 // the user requests a polynomial of order Norder
-psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder,
-                                    psPolynomialType type)
+psPolynomial2D* Polynomial2DAlloc_EAM(psPolynomialType type,
+				      psS32 nXorder, 
+				      psS32 nYorder)
 {
     // PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
@@ -332,128 +271,10 @@
         }
     }
-    psMemSetDeallocator(newPoly, (psFreeFunc) psPolynomial2DFree);
+    psMemSetDeallocator(newPoly, (psFreeFunc) Polynomial2DFree_EAM);
     return(newPoly);
 }
 
-// XXX EAM : BuildSums2D in analogy with BuildSums1D
-static psImage *BuildSums2D(psImage* sums,
-			    psF64 x,      psF64 y,
-			    psS32 nXterm, psS32 nYterm)
-{
-    psS32 nXsum = 0;
-    psS32 nYsum = 0;
-    psF64 xSum = 1.0;
-    psF64 ySum = 1.0;
-
-    nXsum = 2*nXterm;
-    nYsum = 2*nYterm;
-    if (sums == NULL) {
-        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
-    }
-    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
-	psFree (sums);
-        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
-    }
-
-    ySum = 1.0;
-    for (int j = 0; j < nYsum; j++) {
-	xSum = ySum;
-	for (int i = 0; i < nXsum; i++) {
-	    sums->data.F64[i][j] = xSum;
-	    xSum *= x;
-	}
-	ySum *= y;
-    }
-    return (sums);
-}
-
-// XXX EAM : test version of 2d fitting
-psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
-					     psVector* mask,
-					     const psVector* x,
-					     const psVector* y,
-					     const psVector* z,
-					     const psVector* zErr)
-{
-    // I think this is 1 dimension down
-    psImage*  	  A = NULL;
-    psVector* 	  B = NULL;
-    psImage*   Sums = NULL;
-    psF64 wt;
-    psS32 nTerm, nXterm, nYterm;
-
-    nXterm = myPoly->nX + 1;
-    nYterm = myPoly->nY + 1;
-    nTerm = nXterm * nYterm;
-
-    A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
-    B = psVectorAlloc(nTerm, PS_TYPE_F64);
-
-    // Initialize data structures (why is this not a function!)
-    for (int i = 0; i < nTerm; i++) {
-        B->data.F64[i] = 0.0;
-        for (int j = 0; j < nTerm; j++) {
-            A->data.F64[i][j] = 0.0;
-        }
-    }
-
-    // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)
-
-    // Build the B and A data structs.
-    for (int k  = 0; k < x->n; k++) {
-      if ((mask != NULL) && mask->data.U8[k]) continue;
-	Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
-      
-	if (zErr == NULL) {
-	    wt = 1.0;
-	} else {
-	    // this should probably by zErr^2 !!
-	    // this should filter zErr == 0 values
-	    wt = 1.0 / zErr->data.F64[k];
-	}
-
-	// we could skip half of the array and assign at the end
-	// we must handle masked orders
-	for (int n = 0; n < nXterm; n++) {
-	    for (int m = 0; m < nYterm; m++) {
-		B->data.F64[n+m*nXterm] += z->data.F64[k] * Sums->data.F64[n][m] * wt;
-	    }
-	}
-
-	for (int i = 0; i < nXterm; i++) {
-	    for (int j = 0; j < nYterm; j++) {
-		for (int n = 0; n < nXterm; n++) {
-		    for (int m = 0; m < nYterm; m++) {
-			A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
-		    }
-		}
-	    }
-	}
-    }
- 
-    // does the solution in place
-    psGaussJordan (A, B);
-    
-    // XXX: How do we know if these routines were successful?
-    // ALUD = psMatrixLUD(ALUD, &outPerm, A);
-    // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-
-    for (int n = 0; n < nXterm; n++) {
-	for (int m = 0; m < nYterm; m++) {
-	    myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
-	}
-    }
-
-    psFree(A);
-    psFree(B);
-    psFree(Sums);
-
-    psTrace(".psLib.dataManip.VectorFitPolynomial2DOrd", 4,
-            "---- VectorFitPolynomial2DOrd() begin ----\n");
-    return (myPoly);
-}
-
 // write out the terms of the given 2D polynomial
-void psPolynomial2DDump (psPolynomial2D *poly) {
+void Polynomial2DDump_EAM (psPolynomial2D *poly) {
 
     for (int i = 0; i < poly->nX + 1; i++) {
@@ -464,103 +285,7 @@
 }    
 
-psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
-			    const psVector* x,
-			    const psVector* y,
-			    const psVector* z,
-			    const psVector* dz)
-{
-    psVector *X;
-    psVector *Y;
-    psVector *Z;
-    psVector *dZ;
-
-    psVector *zFit   = NULL;
-    psVector *zResid = NULL;
-    psStats  *stats  = NULL;
-
-    X  = psVectorCopy (NULL, x, PS_TYPE_F64);
-    Y  = psVectorCopy (NULL, y, PS_TYPE_F64);
-    Z  = psVectorCopy (NULL, z, PS_TYPE_F64);
-    dZ = psVectorCopy (NULL, dz, PS_TYPE_F64);
-
-    for (int N = 0; N < 3; N++) {
-	// XXX EAM : this would be better defined with an element mask
-	poly   = VectorFitPolynomial2DOrd_EAM (poly, NULL, X, Y, Z, dZ);
-	zFit   = Polynomial2DEvalVectorD (poly, x, y);
-	zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
-
-	stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-	stats  = psVectorStats (stats, zResid, NULL, NULL, 0);
-	psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g (%d pts)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues);
-
-	// re-create X, Y, Z, dZ if pts are valid
-	int n = 0;
-	for (int i = 0; i < zResid->n; i++) {
-	    if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
-	      continue;
-	    } 
-	    X->data.F64[n]  =  x->data.F64[i];
-	    Y->data.F64[n]  =  y->data.F64[i];
-	    Z->data.F64[n]  =  z->data.F64[i];
-	    dZ->data.F64[n] = dz->data.F64[i];
-	    n++;
-	}
-	X->n = n;
-	Y->n = n;
-	Z->n = n;
-	dZ->n = n;
-    }
-    return (poly);
-}
-
-// XXX EAM : be careful here with F32 vs F64 vectors 
-psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
-			    psVector* mask,
-			    const psVector* x,
-			    const psVector* y,
-			    const psVector* z,
-			    const psVector* dz)
-{
-    // PS_VECTOR_CHECK_NULL(mask, NULL);
-    // PS_VECTOR_CHECK_NULL(x, NULL);
-    // PS_VECTOR_CHECK_NULL(y, NULL);
-    // PS_VECTOR_CHECK_NULL(z, NULL);
-    // PS_VECTOR_CHECK_NULL(dz, NULL);
-
-    psVector *zFit   = NULL;
-    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
-    psStats  *stats  = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
-
-    for (int N = 0; N < 3; N++) {
-	poly   = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz);
-	zFit   = Polynomial2DEvalVectorD (poly, x, y);
-	zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
-
-	stats  = psVectorStats (stats, zResid, NULL, mask, 1);
-	psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n", 
-		 stats->sampleMean, stats->sampleStdev);
-
-	// set mask if pts are not valid
-	// we are masking out any point which is out of range 
-	// recovery is not allowed with this scheme
-	for (int i = 0; i < zResid->n; i++) {
-	  if (mask->data.U8[i]) continue;
-	  if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
-	    mask->data.U8[i] = 1;
-	    continue;
-	  }       
-	}
-	psFree (zFit);
-    }
-    psFree (zResid);
-    psFree (stats);
-    return (poly);
-}
-
-// XXX EAM : VectorFitPolynomial2DOrd and Polynomial2DEvalVector require different types (F32 vs F64)
-
-psF64 Polynomial2DEval(const psPolynomial2D* myPoly,
-		       psF64 x,
-		       psF64 y)
+psF64 Polynomial2DEval_EAM(const psPolynomial2D* myPoly,
+			   psF64 x,
+			   psF64 y)
 {
     // PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -587,2 +312,228 @@
     return(polySum);
 }
+
+// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
+psVector *Polynomial2DEvalVector_EAM(const psPolynomial2D *myPoly,
+				     const psVector *x,
+				     const psVector *y)
+
+{
+    // PS_POLY_CHECK_NULL(myPoly, NULL);
+    // PS_VECTOR_CHECK_NULL(x, NULL);
+    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    // PS_VECTOR_CHECK_NULL(y, NULL);
+    // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+
+    psVector *tmp;
+    psS32 vecLen=x->n;
+
+    // Determine the length of the output vector to by the minimum of the x,y vectors
+    if (y->n < vecLen) {
+        vecLen = y->n;
+    }
+
+    // Create output vector to return
+    tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
+
+    // Evaluate the polynomial at the specified points
+    for (psS32 i=0; i<vecLen; i++) {
+        tmp->data.F32[i] = Polynomial2DEval_EAM(myPoly, x->data.F32[i], y->data.F32[i]);
+    }
+
+    // Return output vector
+    return(tmp);
+}
+
+// XXX EAM : this version uses the F64 vectors
+psVector *Polynomial2DEvalVectorD_EAM(const psPolynomial2D *myPoly,
+				  const psVector *x,
+				  const psVector *y)
+
+{
+    // PS_POLY_CHECK_NULL(myPoly, NULL);
+    // PS_VECTOR_CHECK_NULL(x, NULL);
+    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    // PS_VECTOR_CHECK_NULL(y, NULL);
+    // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+
+    psVector *tmp;
+    psS32 vecLen=x->n;
+
+    // Determine the length of the output vector to by the minimum of the x,y vectors
+    if (y->n < vecLen) {
+        vecLen = y->n;
+    }
+
+    // Create output vector to return
+    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
+
+    // Evaluate the polynomial at the specified points
+    for (psS32 i=0; i<vecLen; i++) {
+        tmp->data.F64[i] = Polynomial2DEval_EAM(myPoly,x->data.F64[i],y->data.F64[i]);
+    }
+
+    // Return output vector
+    return(tmp);
+}
+
+// XXX EAM : BuildSums2D in analogy with BuildSums1D
+static psImage *BuildSums2D_EAM(psImage* sums,
+				psF64 x,      psF64 y,
+				psS32 nXterm, psS32 nYterm)
+{
+    psS32 nXsum = 0;
+    psS32 nYsum = 0;
+    psF64 xSum = 1.0;
+    psF64 ySum = 1.0;
+
+    nXsum = 2*nXterm;
+    nYsum = 2*nYterm;
+    if (sums == NULL) {
+        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
+    }
+    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
+	psFree (sums);
+        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
+    }
+
+    ySum = 1.0;
+    for (int j = 0; j < nYsum; j++) {
+	xSum = ySum;
+	for (int i = 0; i < nXsum; i++) {
+	    sums->data.F64[i][j] = xSum;
+	    xSum *= x;
+	}
+	ySum *= y;
+    }
+    return (sums);
+}
+
+// XXX EAM : test version of 2d fitting
+psPolynomial2D* VectorFitPolynomial2D_EAM(psPolynomial2D* myPoly,
+					     psVector* mask,
+					     psMaskType maskValue,
+					     const psVector* z,
+					     const psVector* zErr,
+					     const psVector* x,
+					     const psVector* y)
+{
+    // I think this is 1 dimension down
+    psImage*  	  A = NULL;
+    psVector* 	  B = NULL;
+    psImage*   Sums = NULL;
+    psF64 wt;
+    psS32 nTerm, nXterm, nYterm;
+
+    nXterm = myPoly->nX + 1;
+    nYterm = myPoly->nY + 1;
+    nTerm = nXterm * nYterm;
+
+    A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
+    B = psVectorAlloc(nTerm, PS_TYPE_F64);
+
+    // Initialize data structures (why is this not a function!)
+    for (int i = 0; i < nTerm; i++) {
+        B->data.F64[i] = 0.0;
+        for (int j = 0; j < nTerm; j++) {
+            A->data.F64[i][j] = 0.0;
+        }
+    }
+
+    // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)
+
+    // Build the B and A data structs.
+    for (int k  = 0; k < x->n; k++) {
+	if ((mask != NULL) && (mask->data.U8[k] & maskValue)) continue;
+	Sums = BuildSums2D_EAM(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
+      
+	if (zErr == NULL) {
+	    wt = 1.0;
+	} else {
+	    // this should probably by zErr^2 !!
+	    // this should filter zErr == 0 values
+	    wt = 1.0 / zErr->data.F64[k];
+	}
+
+	// we could skip half of the array and assign at the end
+	// we must handle masked orders
+	for (int n = 0; n < nXterm; n++) {
+	    for (int m = 0; m < nYterm; m++) {
+		B->data.F64[n+m*nXterm] += z->data.F64[k] * Sums->data.F64[n][m] * wt;
+	    }
+	}
+
+	for (int i = 0; i < nXterm; i++) {
+	    for (int j = 0; j < nYterm; j++) {
+		for (int n = 0; n < nXterm; n++) {
+		    for (int m = 0; m < nYterm; m++) {
+			A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
+		    }
+		}
+	    }
+	}
+    }
+ 
+    // does the solution in place
+    psGaussJordan (A, B);
+    
+    // XXX: How do we know if these routines were successful?
+    // ALUD = psMatrixLUD(ALUD, &outPerm, A);
+    // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
+
+    for (int n = 0; n < nXterm; n++) {
+	for (int m = 0; m < nYterm; m++) {
+	    myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
+	}
+    }
+
+    psFree(A);
+    psFree(B);
+    psFree(Sums);
+
+    psTrace(".psLib.dataManip.VectorFitPolynomial2DOrd", 4,
+            "---- VectorFitPolynomial2DOrd() begin ----\n");
+    return (myPoly);
+}
+
+// XXX EAM : be careful here with F32 vs F64 vectors 
+psPolynomial2D* VectorClipFitPolynomial2D_EAM(psPolynomial2D* poly,
+					      psStats *stats,
+					      psVector* mask,
+					      psMaskType maskValue,
+					      const psVector* z,
+					      const psVector* dz,
+					      const psVector* x,
+					      const psVector* y)
+{
+    psVector *zFit   = NULL;
+    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
+
+    // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats:
+    stats->options |= (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+
+    for (int N = 0; N < 3; N++) {
+	poly   = VectorFitPolynomial2D_EAM (poly, mask, maskValue, z, dz, x, y);
+	zFit   = Polynomial2DEvalVectorD_EAM (poly, x, y);
+	zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
+
+	stats  = psVectorStats (stats, zResid, NULL, mask, maskValue);
+	psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n", 
+		 stats->sampleMean, stats->sampleStdev);
+
+	// set mask if pts are not valid
+	// we are masking out any point which is out of range 
+	// recovery is not allowed with this scheme
+	for (int i = 0; i < zResid->n; i++) {
+	    if (mask->data.U8[i]) continue;
+	    if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
+		mask->data.U8[i] &= 0x01;
+		continue;
+	    }       
+	}
+	psFree (zFit);
+    }
+    psFree (zResid);
+    return (poly);
+}
+
+// XXX EAM : VectorFitPolynomial2D and Polynomial2DEvalVector require different types (F32 vs F64)
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 4976)
+++ /trunk/psphot/src/psphot.h	(revision 4977)
@@ -3,5 +3,10 @@
 # include <unistd.h>   // for unlink
 # include <pslib.h>
+# include "psLibUtils.h"
 # include "pmObjects_EAM.h"
+# include "psModulesUtils.h"
+# include "pmPSF.h"
+# include "pmPSFtry.h"
+# include "pmModelGroup.h"
 
 typedef struct {
@@ -11,34 +16,4 @@
     psMetadata *header;
 } psImageData;
-
-// data to test a given PSF model type
-typedef struct {
-    pmModelType modelType;
-    pmPSF      *psf;
-    psArray    *sources;      // pointers to the original sources
-    psArray    *modelFLT;     // model fits, floating parameters 
-    psArray    *modelPSF;     // model fits, PSF parameters
-    psVector   *mask;
-    psVector   *metric;
-    psVector   *fitMag;
-    float       ApResid;
-    float       dApResid;
-    float       skyBias;
-} pmPSF_Test;
-
-// structure to carry a dynamic string
-typedef struct {
-    int NLINE;
-    int Nline;
-    char *line;
-} psLine;
-
-enum {
-    PSPHOT_MASK_INVALID   = 0x01,
-    PSPHOT_MASK_SATURATED = 0x02,
-    PSPHOT_MASK_MARKED    = 0x08,
-} MaskValues;
-
-# define psMemCopy(A)(psMemIncrRefCounter((A)))
 
 // top-level psphot functions
@@ -56,16 +31,9 @@
 bool         psphotSubtractPSF (pmSource *source);
 int 	     psphotSortBySN (const void **a, const void **b);
+int          psphotSaveImage (psMetadata *header, psImage *image, char *filename);
+bool 	     psphotDefinePixels (pmSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius);
 
 // psImageData functions
 psImageData *psImageDataAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header);
-
-// psf utilities
-pmPSF       *pmPSFAlloc (pmModelType type);
-pmPSF_Test  *pmPSF_TestAlloc (psArray *stars, char *modelName);
-pmPSF_Test  *pmPSF_TestModel (psArray *sources, char *modelName, float radius);
-bool	     pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
-pmModel	    *pmModelFromPSF (pmModel *model, pmPSF *psf);
-bool	     pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask);
-bool	     pmPSFMetricModel (pmPSF_Test *test, float RADIUS);
 
 // output functions
@@ -81,74 +49,2 @@
 bool 	     pmModelWriteFLTs (psArray *sources, char *filename);
 bool 	     pmModelWriteNULLs (psArray *sources, char *filename);
-
-// psModule extra utilities
-bool 	     pmSourceDefinePixels(pmSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius);
-bool 	     pmSourceLocalSky_EAM (pmSource *source, psStatsOptions statsOptions, psF32 Radius);
-bool 	     pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF);
-bool 	     pmSourceMoments_EAM(pmSource *source, psF32 radius);
-bool 	     pmModelFitStatus (pmModel *model);
-int	     pmSourceDophotType (pmSource *source);
-
-// minimize 
-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);
-void 	     psTimerFree ();
-bool 	     psTimerClear (char *name);
-psF64 	     psTimerMark (char *name);
-psS32 	     psLogArguments (int *argc, char **argv);
-psS32 	     psTraceArguments (int *argc, char **argv);
-int	     psArgumentGet (int argc, char **argv, char *arg);
-int	     psArgumentRemove (int N, int *argc, char **argv);
-psF32        pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
-psVector    *psVectorCreate (double lower, double upper, double delta, psElemType type);
-psVector    *psGetRowVectorFromImage(psImage *image, psU32 row);
-int          psImageCountPixelMask (psImage *mask, psU8 value);
-
-// basic image functions
-bool         psImageInit (psImage *image,...);
-void	     psImageSmooth_EAM (psImage *image, float sigma, float Nsigma);
-int          psphotSaveImage (psMetadata *header, psImage *image, char *filename);
-
-// psLine functions
-psLine      *psLineAlloc (int Nline);
-bool	     psLineInit (psLine *line);
-bool	     psLineAdd (psLine *line, char *format, ...);
-
-// polynomial functions
-psF64	        Polynomial2DEval(const psPolynomial2D* myPoly, psF64 x, psF64 y);
-psImage	       *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
-psPolynomial2D *VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
-psPolynomial2D *Polynomial2DAlloc(psS32 nXorder, psS32 nYorder, psPolynomialType type);
-void	        psPolynomial2DDump (psPolynomial2D *poly);
-psPolynomial2D *RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
-psPolynomial2D *RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
-psVector       *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
-
-psVector       *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm);
-void	        psPolynomial1DDump (psPolynomial1D *poly);
-psF64		Polynomial1DEval_EAM(psF64 x, const psPolynomial1D* myPoly);
-psVector       *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x);
-psPolynomial1D *Polynomial1DAlloc(psS32 nOrder, psPolynomialType type);
-psPolynomial1D *VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr);
-
-// XXX EAM : bug : missing from psLib/*.h
-bool psMinimizeGaussNewtonDelta (psVector *delta,
-                                 const psVector *params,
-                                 const psVector *paramMask,
-                                 const psArray  *x,
-                                 const psVector *y,
-                                 const psVector *yErr,
-                                 psMinimizeLMChi2Func func);
-
-// not included in the .h file?
-bool p_psVectorPrint (
-    int fd,                            ///< output file descriptor
-    psVector *a,                       ///< vector to print
-    char *name                         ///< name of vector (for title)
-);
-
-bool p_psVectorPrintRow (int fd, psVector *a, char *name);
Index: /trunk/psphot/src/psphotApplyPSF.c
===================================================================
--- /trunk/psphot/src/psphotApplyPSF.c	(revision 4976)
+++ /trunk/psphot/src/psphotApplyPSF.c	(revision 4977)
@@ -62,7 +62,9 @@
 	if (model->radius > OUTER_RADIUS) {
 	  // (re)-allocate image, noise, mask arrays for each peak (square of radius OUTER)
-	  pmSourceDefinePixels (source, imdata, x, y, model->radius);
+	    psphotDefinePixels (source, imdata, x, y, model->radius);
 	}
 
+	if (i > 66) psTraceSetLevel (".psLib.dataManip.psMinimizeLMChi2", 5);
+	    
 	// fit PSF model (set/unset the pixel mask)
 	psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
Index: /trunk/psphot/src/psphotArguments.c
===================================================================
--- /trunk/psphot/src/psphotArguments.c	(revision 4976)
+++ /trunk/psphot/src/psphotArguments.c	(revision 4977)
@@ -9,5 +9,5 @@
 
   // basic pslib options
-  fprintf (stderr, "starting...\n");
+  fprintf (stderr, "starting... %s\n", psLibVersion());
   psLogSetFormat ("M");
   psLogArguments (argc, argv);
Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 4976)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 4977)
@@ -1,5 +1,5 @@
 # include "psphot.h"
 
-// test PSF models and select best option
+// try PSF models and select best option
 
 pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *skystats) 
@@ -7,5 +7,5 @@
     bool       	    status;
     char       	   *modelName;
-    pmPSF_Test 	   *test  = NULL;
+    pmPSFtry 	   *try   = NULL;
     psArray    	   *stars = NULL;
     psMetadataItem *item  = NULL;
@@ -34,14 +34,14 @@
     psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); 
 
-    // set up an array to store the test results
-    psArray *tests = psArrayAlloc (list->n);
+    // set up an array to store the results
+    psArray *models = psArrayAlloc (list->n);
 
-    // test each model option listed in config
-    for (int i = 0; i < tests->n; i++) { 
+    // try each model option listed in config
+    for (int i = 0; i < models->n; i++) { 
 
 	item = psListGetAndIncrement (iter);
 	modelName = item->data.V;
 
-	tests->data[i] = pmPSF_TestModel (stars, modelName, RADIUS);
+	models->data[i] = pmPSFtryModel (stars, modelName, RADIUS);
 	psFree (modelName);
 	psFree (item);
@@ -51,12 +51,12 @@
     psFree (stars);
 
-    // select the best of the test models
+    // select the best of the models
     // here we are using the clippedStdev on the metric as the indicator
-    test = tests->data[0];
+    try = models->data[0];
     int bestN = 0;
-    float bestM = test->dApResid;
-    for (int i = 1; i < tests->n; i++) {
-	test = tests->data[i];
-	float M = test->dApResid;
+    float bestM = try->dApResid;
+    for (int i = 1; i < models->n; i++) {
+	try = models->data[i];
+	float M = try->dApResid;
 	if (M < bestM) {
 	    bestM = M;
@@ -65,29 +65,12 @@
     }
 
-    // keep only the selected test:
-    test = tests->data[bestN];
-    modelName = pmModelGetType (test->modelType);
-    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, test->ApResid, test->dApResid);
+    // keep only the selected model:
+    try = models->data[bestN];
+    modelName = pmModelGetType (try->modelType);
+    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, try->ApResid, try->dApResid);
 
-    // XXX I am having trouble tracking down the double frees vs leaks
-    # if (0)
-    // we are keeping the modelPSF fits from the PSF test, 
-    // but we replace these later on (apply_psf_model)
-    // set source->model based on best psf model fit
-    for (int i = 0; i < test->sources->n; i++) {
-	pmSource *source = test->sources->data[i];
-	// drop masked sources from PSFSTAR list
-	if (test->mask->data.U8[i]) {
-	  source->type = PS_SOURCE_OTHER;
-	  source->modelPSF = NULL;
-	} else {
-	  // source->modelPSF = psMemCopy(test->modelPSF->data[i]);
-	  source->modelPSF = test->modelPSF->data[i];
-	}
-    }
-    # endif
-
-    pmPSF *psf = psMemCopy(test->psf);
-    psFree (tests);
+    // keep only the pmPSF resulting from this analysis
+    pmPSF *psf = psMemCopy(try->psf);
+    psFree (models);
     return (psf);
 }
Index: /trunk/psphot/src/psphotDefinePixels.c
===================================================================
--- /trunk/psphot/src/psphotDefinePixels.c	(revision 4977)
+++ /trunk/psphot/src/psphotDefinePixels.c	(revision 4977)
@@ -0,0 +1,21 @@
+# include "psphot.h"
+
+bool psphotDefinePixels(pmSource *mySource, 
+			const psImageData *imdata,
+			psF32 x, 
+			psF32 y,
+			psF32 Radius)
+{
+    psRegion srcRegion;
+
+    // Grab a subimage of the original image of size (2 * outerRadius).
+    srcRegion = psRegionForSquare (x, y, Radius);
+    srcRegion = psRegionForImage (imdata->image, &srcRegion);
+
+    mySource->pixels = psImageSubset(imdata->image, srcRegion);
+    mySource->noise  = psImageSubset(imdata->noise, srcRegion);
+    mySource->mask   = psImageSubset(imdata->mask,  srcRegion);
+
+    return(mySource);
+}
+
Index: /trunk/psphot/src/psphotFitGalaxies.c
===================================================================
--- /trunk/psphot/src/psphotFitGalaxies.c	(revision 4976)
+++ /trunk/psphot/src/psphotFitGalaxies.c	(revision 4977)
@@ -72,5 +72,5 @@
 	if (model->radius > OUTER_RADIUS) {
 	  // (re)-allocate image, noise, mask arrays for each peak (square of radius OUTER)
-	  pmSourceDefinePixels (source, imdata, x, y, model->radius);
+	  psphotDefinePixels (source, imdata, x, y, model->radius);
 	}
 
Index: /trunk/psphot/src/psphotSourceStats.c
===================================================================
--- /trunk/psphot/src/psphotSourceStats.c	(revision 4976)
+++ /trunk/psphot/src/psphotSourceStats.c	(revision 4977)
@@ -24,5 +24,5 @@
 	// allocate image, noise, mask arrays for each peak (square of radius OUTER)
 	// XXX EAM : this has 12 leaks (v.5)
-	pmSourceDefinePixels (source, imdata, source->peak->x, source->peak->y, OUTER);
+	psphotDefinePixels (source, imdata, source->peak->x, source->peak->y, OUTER);
 
 	// measure a local sky value
Index: /trunk/psphot/src/pspsf.c
===================================================================
--- /trunk/psphot/src/pspsf.c	(revision 4976)
+++ /trunk/psphot/src/pspsf.c	(revision 4977)
@@ -1,39 +1,12 @@
 # include "psphot.h"
 
-static void pmPSFFree (pmPSF *psf) {
-
-  if (psf == NULL) return;
-
-  psFree (psf->params);
-  return;
-}
-
-// a PSF always has 4 parameters fewer than the equivalent model
-pmPSF *pmPSFAlloc (pmModelType type) {
-
-    int Nparams;
-
-    pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF));
-
-    psf->type = type;
-    psf->chisq = 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++) {
-	// we need a way to set the allowed range of orders
-	psf->params->data[i] = Polynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
-    }
-
-    psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
-    return(psf);
-}
-
-static void pmPSF_TestFree (pmPSF_Test *test) {
+// ********  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;
@@ -49,7 +22,8 @@
 }
 
-pmPSF_Test *pmPSF_TestAlloc (psArray *sources, char *modelName) {
-
-    pmPSF_Test *test = (pmPSF_Test *) psAlloc(sizeof(pmPSF_Test));
+// allocate a pmPSFtry based on the desired sources and the model (identified by name)
+pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName) {
+
+    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
 
     // XXX probably need to increment ref counter
@@ -74,14 +48,23 @@
     }	
 
-    psMemSetDeallocator(test, (psFreeFunc) pmPSF_TestFree);
+    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
     return (test);
 }
 
-// test->mask values indicate reason source was rejected:
+// 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:
 // 1: outlier in psf polynomial fit
 // 2: flt model failed to converge
 // 3: psf model failed to converge
 // 4: invalid source photometry
-pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS) 
+// XXX EAM : use an enum for these values?
+
+pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS) 
 {
     bool status;
@@ -93,12 +76,12 @@
     int Npsf = 0;
 
-    pmPSF_Test *test = pmPSF_TestAlloc (sources, modelName);
+    pmPSFtry *try = pmPSFtryAlloc (sources, modelName);
 
     // stage 1:  fit an independent model (freeModel) to all sources
     psTimerStart ("fit");
-    for (int i = 0; i < test->sources->n; i++) {
-
-	pmSource *source = test->sources->data[i];
-	pmModel  *model  = pmSourceModelGuess (source, test->modelType); 
+    for (int i = 0; i < try->sources->n; i++) {
+
+	pmSource *source = try->sources->data[i];
+	pmModel  *model  = pmSourceModelGuess (source, try->modelType); 
 	x = source->peak->x;
 	y = source->peak->y;
@@ -112,31 +95,31 @@
 	// exclude the poor fits
 	if (!status) {
-	  test->mask->data.U8[i] = 2;
+	  try->mask->data.U8[i] = 2;
 	  psFree (model);
 	  continue;
 	}
-	test->modelFLT->data[i] = model;
+	try->modelFLT->data[i] = model;
 	Nflt ++;
     }
-    psLogMsg ("psphot.psftest", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
-    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
+    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 (test->modelFLT, "modelsFLT.dat");
+    // DumpModelFits (try->modelFLT, "modelsFLT.dat");
 
     // stage 2: construct a psf (pmPSF) from this collection of model fits
-    pmPSFFromModels (test->psf, test->modelFLT, test->mask);
+    pmPSFFromModels (try->psf, try->modelFLT, try->mask);
 
     // stage 3: refit with fixed shape parameters
     psTimerStart ("fit");
-    for (int i = 0; i < test->sources->n; i++) {
+    for (int i = 0; i < try->sources->n; i++) {
 	// masked for: bad model fit, outlier in parameters
-	if (test->mask->data.U8[i]) continue;
-
-	pmSource *source = test->sources->data[i];
-	pmModel  *modelFLT = test->modelFLT->data[i];
+	if (try->mask->data.U8[i]) 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, test->psf); 
+	pmModel *modelPSF = pmModelFromPSF (modelFLT, try->psf); 
 	x = source->peak->x;
 	y = source->peak->y;
@@ -147,5 +130,5 @@
 	// skip poor fits
 	if (!status) {
-	    test->mask->data.U8[i] = 3;
+	    try->mask->data.U8[i] = 3;
 	    psFree (modelPSF);
 	    goto next_source;
@@ -153,5 +136,5 @@
 
 	// otherwise, save the resulting model
-	test->modelPSF->data[i] = modelPSF;
+	try->modelPSF->data[i] = modelPSF;
 
 	// XXX : use a different aperture radius from the fit radius?
@@ -159,10 +142,10 @@
 	// XXX : pass 'source' as input?
 	if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
-	    test->mask->data.U8[i] = 4;
+	    try->mask->data.U8[i] = 4;
 	    goto next_source;
 	}	    
 
-	test->metric->data.F64[i] = obsMag - fitMag;
-	test->fitMag->data.F64[i] = fitMag;
+	try->metric->data.F64[i] = obsMag - fitMag;
+	try->fitMag->data.F64[i] = fitMag;
 	Npsf ++;
 
@@ -171,86 +154,25 @@
 
     }
-    psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
-    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
+    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 (test->modelPSF, "modelsPSF.dat");
-
-    // XXX this function wants aperture radius from pmSourcePhotometry
-    pmPSFMetricModel (test, RADIUS);
-    psLogMsg ("psphot.pspsf", 3, "test model %s, ap-fit: %f +/- %f, sky bias: %f\n", 
-	      modelName, test->ApResid, test->dApResid, test->skyBias);
-
-    return (test);
-}
-
-// input: an array of pmModels, pre-allocated psf
-// some of the array entries may be NULL, 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
-	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?
-
-    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
-	}
-
-	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 (x);
-    psFree (y);
-    psFree (z);
-    psFree (dz);
-    return (true);
-}
-
-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);
-
-    // do we need a different model for the PSF vs FLT version?
-    pmModel *modelPSF = pmModelAlloc (psf->type);
-
-    // set model parameters for this source based on PSF information
-    modelFromPSFFunc (modelPSF, modelFLT, psf);
-
-    return (modelPSF);
-}
-
-bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS) {
+    // 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->ApResid, try->dApResid, try->skyBias);
+
+    return (try);
+}
+
+
+bool pmPSFtryMetric (pmPSFtry *try, float RADIUS) {
 
   float dBin;
   int   nKeep, nSkip;
 
-  // the measured (aperture - fit) magnitudes (dA == test->metric)
+  // 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
@@ -262,13 +184,13 @@
 
   // rflux = ten(0.4*fitMag);
-  psVector *rflux = psVectorAlloc (test->sources->n, PS_TYPE_F64);
-  for (int i = 0; i < test->sources->n; i++) {
-    if (test->mask->data.U8[i]) continue;
-    rflux->data.F64[i] = pow(10.0, 0.4*test->fitMag->data.F64[i]);
+  psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
+  for (int i = 0; i < try->sources->n; i++) {
+    if (try->mask->data.U8[i]) 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, test->mask, 0xff);
+  psVectorStats (stats, rflux, NULL, try->mask, 0xff);
   
   // build binned versions of rflux, metric
@@ -284,14 +206,14 @@
   for (int i = 0; i < daBin->n; i++) {
 
-    psVector *tmp = psVectorAlloc (test->sources->n, PS_TYPE_F64);
+    psVector *tmp = psVectorAlloc (try->sources->n, PS_TYPE_F64);
     tmp->n = 0;
 
     // accumulate data within bin range
-    for (int j = 0; j < test->sources->n; j++) {
+    for (int j = 0; j < try->sources->n; j++) {
       // masked for: bad model fit, outlier in parameters
-      if (test->mask->data.U8[j]) continue;
+      if (try->mask->data.U8[j]) continue;
     
       // skip points with extreme dA values
-      if (fabs(test->metric->data.F64[j]) > 0.5) continue;
+      if (fabs(try->metric->data.F64[j]) > 0.5) continue;
 
       // skip points outside of this bin
@@ -299,5 +221,5 @@
       if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
 
-      tmp->data.F64[tmp->n] = test->metric->data.F64[j];
+      tmp->data.F64[tmp->n] = try->metric->data.F64[j];
       tmp->n ++;
     }
@@ -342,8 +264,13 @@
 
   // linear fit to rfBin, daBin
-  psPolynomial1D *poly = Polynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+  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 = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL);
 
-  psVector *daBinFit = Polynomial1DEvalVector_EAM (poly, rfBin);
+  psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
   psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
 
@@ -351,7 +278,7 @@
   stats = psVectorStats (stats, daResid, NULL, maskB, 1);
 
-  test->ApResid = poly->coeff[0];
-  test->dApResid = stats->clippedStdev;
-  test->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+  try->ApResid = poly->coeff[0];
+  try->dApResid = stats->clippedStdev;
+  try->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
 
   psFree (rflux);