Index: /branches/eam_rel9_p0/psModules/src/objects/Makefile.am
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/Makefile.am	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/Makefile.am	(revision 5958)
@@ -15,5 +15,5 @@
 	models/pmModel_QGAUSS.c \
 	models/pmModel_SGAUSS.c 
-    
+
 psmoduleincludedir = $(includedir)
 psmoduleinclude_HEADERS = \
Index: /branches/eam_rel9_p0/psModules/src/objects/models/pmModel_PGAUSS.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 5958)
@@ -26,13 +26,13 @@
 
     if (deriv != NULL) {
-        // note difference from a pure gaussian: q = PAR[1]*r
+        psF32 *dPAR = deriv->data.F32;
         psF32 q = PAR[1]*r*r*t;
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
-        deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
+        dPAR[0] = +1.0;
+        dPAR[1] = +r;
+        dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
+        dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
+        dPAR[4] = -2.0*q*px*X;
+        dPAR[5] = -2.0*q*py*Y;
+        dPAR[6] = -q*X*Y;
     }
     return(f);
@@ -47,5 +47,5 @@
 
     beta_lim[0][0].data.F32[0] = 1000;
-    beta_lim[0][0].data.F32[1] = 10000;
+    beta_lim[0][0].data.F32[1] = 3e6;
     beta_lim[0][0].data.F32[2] = 5;
     beta_lim[0][0].data.F32[3] = 5;
@@ -63,5 +63,5 @@
 
     params_max[0][0].data.F32[0] = 1e5;
-    params_max[0][0].data.F32[1] = 1e6;
+    params_max[0][0].data.F32[1] = 1e8;
     params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
     params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
@@ -130,4 +130,5 @@
     params[5] = 1.2 / moments->Sy;
     params[6] = 0.0;
+
     return(true);
 }
Index: /branches/eam_rel9_p0/psModules/src/objects/models/pmModel_QGAUSS.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 5958)
@@ -31,16 +31,18 @@
 
     if (deriv != NULL) {
+        psF32 *dPAR = deriv->data.F32;
+
         // note difference from a pure gaussian: q = params->data.F32[1]*r
         psF32 t = PAR[1]*r*r;
         psF32 q = t*(PAR[7] + 2.25*pow(z, 1.25));
 
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
-        deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-        deriv->data.F32[7] = -t*z;
+        dPAR[0] = +1.0;
+        dPAR[1] = +r;
+        dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y);
+        dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
+        dPAR[4] = -2.0*q*px*X;
+        dPAR[5] = -2.0*q*py*Y;
+        dPAR[6] = -q*X*Y;
+        dPAR[7] = -t*z;
     }
     return(f);
@@ -55,5 +57,5 @@
 
     beta_lim[0][0].data.F32[0] = 1000;
-    beta_lim[0][0].data.F32[1] = 10000;
+    beta_lim[0][0].data.F32[1] = 3e6;
     beta_lim[0][0].data.F32[2] = 5;
     beta_lim[0][0].data.F32[3] = 5;
@@ -73,5 +75,5 @@
 
     params_max[0][0].data.F32[0] = 1e5;
-    params_max[0][0].data.F32[1] = 1e6;
+    params_max[0][0].data.F32[1] = 1e8;
     params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
     params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
@@ -198,6 +200,6 @@
     status &= ((dPAR[1]/PAR[1]) < 0.5);
 
-    if (status)
-        return true;
-    return false;
-}
+    if (!status)
+        return false;
+    return true;
+}
Index: /branches/eam_rel9_p0/psModules/src/objects/pmModelGroup.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/pmModelGroup.c	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/pmModelGroup.c	(revision 5958)
@@ -10,14 +10,59 @@
 #include "models/pmModel_SGAUSS.c"
 
-static pmModelGroup models[] = {
-                                   {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
-                                   {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
-                                   {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
-                                   {"PS_MODEL_SGAUSS",       9, pmModelFunc_SGAUSS,  pmModelFlux_SGAUSS,  pmModelRadius_SGAUSS,  pmModelLimits_SGAUSS,  pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS},
-                               };
+static pmModelGroup defaultModels[] = {
+                                          {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
+                                          {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
+                                          {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
+                                          {"PS_MODEL_SGAUSS",       9, pmModelFunc_SGAUSS,  pmModelFlux_SGAUSS,  pmModelRadius_SGAUSS,  pmModelLimits_SGAUSS,  pmModelGuess_SGAUSS, pmModelFromPSF_SGAUSS, pmModelFitStatus_SGAUSS},
+                                      };
+
+static pmModelGroup *models = NULL;
+static int Nmodels = 0;
+
+static void ModelGroupFree (pmModelGroup *modelGroup)
+{
+
+    if (modelGroup == NULL)
+        return;
+    psFree (modelGroup);
+    return;
+}
+
+pmModelGroup *pmModelGroupAlloc (int nModels)
+{
+
+    pmModelGroup *modelGroup = (pmModelGroup *) psAlloc (nModels * sizeof(pmModelGroup));
+    psMemSetDeallocator(modelGroup, (psFreeFunc) ModelGroupFree);
+    return (modelGroup);
+}
+
+void pmModelGroupAdd (pmModelGroup *model)
+{
+
+    if (models == NULL) {
+        pmModelGroupInit ();
+    }
+
+    Nmodels ++;
+    models = (pmModelGroup *) psRealloc (models, Nmodels*sizeof(pmModelGroup));
+    models[Nmodels-1] = model[0];
+    return;
+}
+
+void pmModelGroupInit (void)
+{
+
+    int Nnew = sizeof (defaultModels) / sizeof (pmModelGroup);
+
+    models = pmModelGroupAlloc (Nnew);
+    for (int i = 0; i < Nnew; i++) {
+        models[i] = defaultModels[i];
+    }
+    Nmodels = Nnew;
+    return;
+}
 
 pmModelFunc pmModelFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -29,5 +74,4 @@
 pmModelFlux pmModelFlux_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -39,5 +83,4 @@
 pmModelRadius pmModelRadius_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -49,5 +92,4 @@
 pmModelLimits pmModelLimits_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -59,5 +101,4 @@
 pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -69,5 +110,4 @@
 pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -79,5 +119,4 @@
 pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -89,5 +128,4 @@
 psS32 pmModelParameterCount (pmModelType type)
 {
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
@@ -99,6 +137,4 @@
 psS32 pmModelSetType (char *name)
 {
-
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     for (int i = 0; i < Nmodels; i++) {
         if (!strcmp(models[i].name, name)) {
@@ -111,6 +147,4 @@
 char *pmModelGetType (pmModelType type)
 {
-
-    int Nmodels = sizeof (models) / sizeof (pmModelGroup);
     if ((type < 0) || (type >= Nmodels)) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
Index: /branches/eam_rel9_p0/psModules/src/objects/pmModelGroup.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/pmModelGroup.h	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/pmModelGroup.h	(revision 5958)
@@ -9,6 +9,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-10-10 19:53:40 $
+ *  @version $Revision: 1.1.16.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-10 04:46:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -209,3 +209,12 @@
 pmModelGroup;
 
+// allocate a pmModelGroup to hold nModels entries
+pmModelGroup *pmModelGroupAlloc (int nModels);
+
+// initialize the internal (static) model group with the default models
+void pmModelGroupInit (void);
+
+// add a new model to the internal (static) model group
+void pmModelGroupAdd (pmModelGroup *model);
+
 # endif
Index: /branches/eam_rel9_p0/psModules/src/objects/pmObjects.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/pmObjects.c	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/pmObjects.c	(revision 5958)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.5.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-10 04:46:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -271,4 +271,7 @@
     tmp->chisq = 0.0;
     tmp->nIter = 0;
+    tmp->radius = 0;
+    tmp->status = PM_MODEL_UNTRIED;
+
     psS32 Nparams = pmModelParameterCount(type);
     if (Nparams == 0) {
@@ -307,7 +310,9 @@
     tmp->mask = NULL;
     tmp->moments = NULL;
+    tmp->blends = NULL;
     tmp->modelPSF = NULL;
     tmp->modelFLT = NULL;
-    tmp->type = 0;
+    tmp->type = PM_SOURCE_UNKNOWN;
+    tmp->mode = PM_SOURCE_DEFAULT;
     psMemSetDeallocator(tmp, (psFreeFunc) sourceFree);
 
@@ -827,4 +832,5 @@
     psS32 numPixels = 0;
     psF32 Sum = 0.0;
+    psF32 Var = 0.0;
     psF32 X1 = 0.0;
     psF32 Y1 = 0.0;
@@ -850,6 +856,7 @@
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
         for (psS32 col = 0; col < source->pixels->numCols ; col++) {
-            if ((source->mask != NULL) && (source->mask->data.U8[row][col]))
+            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {
                 continue;
+            }
 
             psF32 xDiff = col + source->pixels->col0 - xPeak;
@@ -858,14 +865,18 @@
             // XXX EAM : calculate xDiff, yDiff up front;
             //           radius is just a function of (xDiff, yDiff)
-            if (!VALID_RADIUS(xDiff, yDiff, R2))
+            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
                 continue;
+            }
 
             psF32 pDiff = source->pixels->data.F32[row][col] - sky;
+            psF32 wDiff = source->weight->data.F32[row][col];
 
             // XXX EAM : check for valid S/N in pixel
             // XXX EAM : should this limit be user-defined?
-            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1)
+            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {
                 continue;
-
+            }
+
+            Var += wDiff;
             Sum += pDiff;
             X1  += xDiff * pDiff;
@@ -880,7 +891,9 @@
         }
     }
+
+    // if we have less than (1/4) of the possible pixels, force a retry
     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
-    if ((numPixels < 3) || (Sum <= 0)) {
-        psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
+    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
+        psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
         return (false);
@@ -899,5 +912,5 @@
     y = Y1/Sum;
     if ((fabs(x) > radius) || (fabs(y) > radius)) {
-        psTrace (".psModules.pmSourceMoments", 5,
+        psTrace (".psModules.pmSourceMoments", 3,
                  "large centroid swing; invalid peak %d, %d\n",
                  source->peak->x, source->peak->y);
@@ -912,4 +925,5 @@
     source->moments->Sxy = XY/Sum - x*y;
     source->moments->Sum = Sum;
+    source->moments->SN  = Sum / sqrt(Var);
     source->moments->Peak = peakPixel;
     source->moments->nPixels = numPixels;
@@ -994,8 +1008,16 @@
         psImage *splane = NULL;
         int binX, binY;
+        bool status;
+
+        psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
+        if (!status)
+            SX_MAX = 10.0;
+        psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
+        if (!status)
+            SY_MAX = 10.0;
 
         // construct a sigma-plane image
         // psImageAlloc does zero the data
-        splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
+        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
         for (int i = 0; i < splane->numRows; i++)
         {
@@ -1131,5 +1153,5 @@
 XXX: How can this function ever return FALSE?
  
-XXX EAM : add the saturated mask value to metadata 
+EAM: I moved S/N calculation to pmSourceMoments, using weight image
 *****************************************************************************/
 
@@ -1146,5 +1168,8 @@
     int Ncr      = 0;
     int Nsatstar = 0;
-
+    // psRegion allArray = psRegionSet (0, 0, 0, 0);
+    psRegion inner;
+
+    // report stats on S/N values for star-like objects
     psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
     starsn->n = 0;
@@ -1152,8 +1177,5 @@
     // check return status value (do these exist?)
     bool status;
-    psF32 RDNOISE    = psMetadataLookupF32 (&status, metadata, "RDNOISE");
-    psF32 GAIN       = psMetadataLookupF32 (&status, metadata, "GAIN");
     psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
-    // psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
 
     // XXX allow clump size to be scaled relative to sigmas?
@@ -1168,25 +1190,15 @@
         psF32 sigY = tmpSrc->moments->Sy;
 
-        // calculate and save signal-to-noise estimates
-        psF32 S  = tmpSrc->moments->Sum;
-        psF32 A  = 4 * M_PI * sigX * sigY;
-        psF32 B  = tmpSrc->moments->Sky;
-        psF32 RT = sqrt(S + (A * B) + (A * PS_SQR(RDNOISE) / sqrt(GAIN)));
-        psF32 SN = (S * sqrt(GAIN) / RT);
-        tmpSrc->moments->SN = SN;
-
         // XXX EAM : can we use the value of SATURATE if mask is NULL?
-        //
-        // XXX: Must verify this region (the region argument was added to psImageCountPixelMask()
-        // after EAM wrote this code.
-        //
-        psRegion tmpRegion = psRegionSet(0, tmpSrc->mask->numCols, 0, tmpSrc->mask->numRows);
-        int Nsatpix = psImageCountPixelMask(tmpSrc->mask, tmpRegion, PSPHOT_MASK_SATURATED);
+        inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
+        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);
 
         // saturated star (size consistent with PSF or larger)
         // Nsigma should be user-configured parameter
         bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
+        big = true;
         if ((Nsatpix > 1) && big) {
-            tmpSrc->type = PM_SOURCE_SATSTAR;
+            tmpSrc->type = PM_SOURCE_STAR;
+            tmpSrc->mode = PM_SOURCE_SATSTAR;
             Nsatstar ++;
             continue;
@@ -1196,4 +1208,5 @@
         if (Nsatpix > 1) {
             tmpSrc->type = PM_SOURCE_SATURATED;
+            tmpSrc->mode = PM_SOURCE_DEFAULT;
             Nsat ++;
             continue;
@@ -1205,4 +1218,5 @@
         if ((sigX < 0.05) || (sigY < 0.05)) {
             tmpSrc->type = PM_SOURCE_DEFECT;
+            tmpSrc->mode = PM_SOURCE_DEFAULT;
             Ncr ++;
             continue;
@@ -1212,4 +1226,5 @@
         if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
             tmpSrc->type = PM_SOURCE_GALAXY;
+            tmpSrc->mode = PM_SOURCE_DEFAULT;
             Ngal ++;
             continue;
@@ -1217,12 +1232,13 @@
 
         // the rest are probable stellar objects
-        starsn->data.F32[starsn->n] = SN;
+        starsn->data.F32[starsn->n] = tmpSrc->moments->SN;
         starsn->n ++;
         Nstar ++;
 
-        // PSF star (within 1.5 sigma of clump center
+        // PSF star (within 1.5 sigma of clump center, S/N > limit)
         psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
-        if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
-            tmpSrc->type = PM_SOURCE_PSFSTAR;
+        if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {
+            tmpSrc->type = PM_SOURCE_STAR;
+            tmpSrc->mode = PM_SOURCE_PSFSTAR;
             Npsf ++;
             continue;
@@ -1230,5 +1246,6 @@
 
         // random type of star
-        tmpSrc->type = PM_SOURCE_OTHER;
+        tmpSrc->type = PM_SOURCE_STAR;
+        tmpSrc->mode = PM_SOURCE_DEFAULT;
     }
 
@@ -1551,14 +1568,7 @@
 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
-/******************************************************************************
-pmSourceFitModel(source, model): must create the appropiate arguments to the
-LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions.
- 
-XXX: should there be a mask value?
-XXX EAM : fit the specified model (not necessarily the one in source)
-*****************************************************************************/
-bool pmSourceFitModel_v5(pmSource *source,
-                         pmModel *model,
-                         const bool PSF)
+bool pmSourceFitModel (pmSource *source,
+                       pmModel *model,
+                       const bool PSF)
 {
     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
@@ -1569,11 +1579,12 @@
     PS_ASSERT_PTR_NON_NULL(source->mask, false);
     PS_ASSERT_PTR_NON_NULL(source->weight, false);
+
+    // XXX EAM : is it necessary for the mask & weight to exist?  the
+    //           tests below could be conditions (!NULL)
+
     psBool fitStatus = true;
     psBool onPic     = true;
     psBool rc        = true;
 
-    // XXX EAM : is it necessary for the mask & weight to exist?  the
-    //           tests below could be conditions (!NULL)
-
     psVector *params = model->params;
     psVector *dparams = model->dparams;
@@ -1582,45 +1593,48 @@
     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
 
-    int nParams = PSF ? params->n - 4 : params->n;
-
-    // find the number of valid pixels
-    // XXX EAM : this loop and the loop below could just be one pass
-    //           using the psArrayAdd and psVectorExtend functions
-    psS32 count = 0;
+    int nParams = PSF ? 4 : params->n;
+
+    // maximum number of valid pixels
+    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
+
+    // construct the coordinate and value entries
+    psArray *x = psArrayAlloc(nPix);
+    psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);
+    psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32);
+
+    nPix = 0;
     for (psS32 i = 0; i < source->pixels->numRows; i++) {
         for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                count++;
-            }
-        }
-    }
-    if (count <  nParams + 1) {
+            if (source->mask->data.U8[i][j]) {
+                continue;
+            }
+            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+            // Convert i/j to image space:
+            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
+            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+            x->data[nPix] = (psPtr *) coord;
+            y->data.F32[nPix] = source->pixels->data.F32[i][j];
+
+            // psMinimizeLMChi2 takes wt = 1/dY^2
+            if (source->weight->data.F32[i][j] == 0) {
+                continue;
+            }
+            yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
+            nPix++;
+        }
+    }
+    if (nPix <  nParams + 1) {
         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
         psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        model->status = PM_MODEL_BADARGS;
+        psFree (x);
+        psFree (y);
+        psFree (yErr);
         return(false);
     }
-
-    // construct the coordinate and value entries
-    psArray *x = psArrayAlloc(count);
-    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
-    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
-    psS32 tmpCnt = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-                // XXX: Convert i/j to image space:
-                // XXX EAM: coord order is (x,y) == (col,row)
-                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
-                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
-                x->data[tmpCnt] = (psPtr *) coord;
-                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
-                yErr->data.F32[tmpCnt] = sqrt (source->weight->data.F32[i][j]);
-                // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then
-                //           the minimization function calculates sq()
-                tmpCnt++;
-            }
-        }
-    }
+    x->n = nPix;
+    y->n = nPix;
+    yErr->n = nPix;
 
     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
@@ -1638,150 +1652,5 @@
     }
 
-    // XXX EAM : covar must be F64?
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
-
-    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
-    for (int i = 0; i < dparams->n; i++) {
-        if ((paramMask != NULL) && paramMask->data.U8[i])
-            continue;
-        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
-    }
-
-    // XXX EAM: we need to do something (give an error?) if rc is false
-    // XXX EAM: psMinimizeLMChi2 does not check convergence
-
-    // XXX models can go insane: reject these
-    onPic &= (params->data.F32[2] >= source->pixels->col0);
-    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
-    onPic &= (params->data.F32[3] >= source->pixels->row0);
-    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
-
-    // XXX EAM: save the resulting chisq, nDOF, nIter
-    model->chisq = myMin->value;
-    model->nIter = myMin->iter;
-    model->nDOF  = y->n - nParams;
-
-    // XXX EAM get the Gauss-Newton distance for fixed model parameters
-    if (paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
-        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
-        for (int i = 0; i < dparams->n; i++) {
-            if (!paramMask->data.U8[i])
-                continue;
-            dparams->data.F32[i] = delta->data.F64[i];
-        }
-        psFree (delta);
-    }
-
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myMin);
-    psFree(covar);
-    psFree(paramMask);
-
-    rc = (onPic && fitStatus);
-    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
-    return(rc);
-}
-
-// XXX EAM : new version with parameter range limits and weight enhancement
-bool pmSourceFitModel (pmSource *source,
-                       pmModel *model,
-                       const bool PSF)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->moments, false);
-    PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
-    PS_ASSERT_PTR_NON_NULL(source->weight, false);
-
-    // XXX EAM : is it necessary for the mask & weight to exist?  the
-    //           tests below could be conditions (!NULL)
-
-    psBool fitStatus = true;
-    psBool onPic     = true;
-    psBool rc        = true;
-    psF32  Ro, ymodel;
-
-    psVector *params = model->params;
-    psVector *dparams = model->dparams;
-    psVector *paramMask = NULL;
-
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-
-    // XXX EAM : I need to use the sky value to constrain the weight model
-    int nParams = PSF ? params->n - 4 : params->n;
-    psF32 So = params->data.F32[0];
-
-    // find the number of valid pixels
-    // XXX EAM : this loop and the loop below could just be one pass
-    //           using the psArrayAdd and psVectorExtend functions
-    psS32 count = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                count++;
-            }
-        }
-    }
-    if (count <  nParams + 1) {
-        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
-        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
-        return(false);
-    }
-
-    // construct the coordinate and value entries
-    psArray *x = psArrayAlloc(count);
-    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
-    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
-    psS32 tmpCnt = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            if (source->mask->data.U8[i][j] == 0) {
-                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-                // XXX: Convert i/j to image space:
-                // XXX EAM: coord order is (x,y) == (col,row)
-                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
-                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
-                x->data[tmpCnt] = (psPtr *) coord;
-                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
-
-                // compare observed flux to model flux to adjust weight
-                ymodel = modelFunc (NULL, model->params, coord);
-
-                // this test enhances the weight based on deviation from the model flux
-                Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
-
-                // psMinimizeLMChi2_EAM takes wt = 1/dY^2
-                if (source->weight->data.F32[i][j] == 0) {
-                    yErr->data.F32[tmpCnt] = 0.0;
-                } else {
-                    yErr->data.F32[tmpCnt] = 1.0 / (source->weight->data.F32[i][j] * Ro);
-                }
-                tmpCnt++;
-            }
-        }
-    }
-
-    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
-                            PM_SOURCE_FIT_MODEL_TOLERANCE);
-
-    // PSF model only fits first 4 parameters, FLT model fits all
-    if (PSF) {
-        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
-        for (int i = 0; i < 4; i++) {
-            paramMask->data.U8[i] = 0;
-        }
-        for (int i = 4; i < paramMask->n; i++) {
-            paramMask->data.U8[i] = 1;
-        }
-    }
-
-    // XXX EAM : I've added three types of parameter range checks
-    // XXX EAM : this requires my new psMinimization functions
+    // Set the parameter range checks
     pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
     psVector *beta_lim = NULL;
@@ -1807,19 +1676,10 @@
     }
 
-    // XXX EAM: we need to do something (give an error?) if rc is false
-    // XXX EAM: psMinimizeLMChi2 does not check convergence
-
-    // XXX models can go insane: reject these
-    onPic &= (params->data.F32[2] >= source->pixels->col0);
-    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
-    onPic &= (params->data.F32[3] >= source->pixels->row0);
-    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
-
-    // XXX EAM: save the resulting chisq, nDOF, nIter
+    // save the resulting chisq, nDOF, nIter
     model->chisq = myMin->value;
     model->nIter = myMin->iter;
     model->nDOF  = y->n - nParams;
 
-    // XXX EAM get the Gauss-Newton distance for fixed model parameters
+    // get the Gauss-Newton distance for fixed model parameters
     if (paramMask != NULL) {
         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
@@ -1832,8 +1692,28 @@
     }
 
-    psFree(paramMask);
+    // set the model success or failure status
+    if (!fitStatus) {
+        model->status = PM_MODEL_NONCONVERGE;
+    } else {
+        model->status = PM_MODEL_SUCCESS;
+    }
+
+    // models can go insane: reject these
+    onPic &= (params->data.F32[2] >= source->pixels->col0);
+    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
+    onPic &= (params->data.F32[3] >= source->pixels->row0);
+    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
+    if (!onPic) {
+        model->status = PM_MODEL_OFFIMAGE;
+    }
+
+    source->mode |= PM_SOURCE_FITTED;
+
     psFree(x);
     psFree(y);
+    psFree(yErr);
     psFree(myMin);
+    psFree(covar);
+    psFree(paramMask);
 
     rc = (onPic && fitStatus);
@@ -1846,5 +1726,7 @@
                              pmModel *model,
                              bool center,
-                             psS32 flag)
+                             bool sky,
+                             bool add
+                                )
 {
     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
@@ -1859,4 +1741,6 @@
     psS32 imageCol;
     psS32 imageRow;
+    psF32 skyValue = params->data.F32[0];
+    psF32 pixelValue;
 
     for (psS32 i = 0; i < image->numRows; i++) {
@@ -1864,5 +1748,5 @@
             if ((mask != NULL) && mask->data.U8[i][j])
                 continue;
-            psF32 pixelValue;
+
             // XXX: Should you be adding the pixels for the entire subImage,
             // or a radius of pixels around it?
@@ -1882,15 +1766,21 @@
             x->data.F32[0] = (float) imageCol;
             x->data.F32[1] = (float) imageRow;
-            pixelValue = modelFunc (NULL, params, x);
-            // fprintf (stderr, "%f %f  %d %d  %f\n", x->data.F32[0], x->data.F32[1], i, j, pixelValue);
-
-            if (flag == 1) {
-                pixelValue = -pixelValue;
-            }
-
-            // XXX: Must figure out how to calculate the image coordinates and
-            // how to use the boolean "center" flag.
-
-            image->data.F32[i][j]+= pixelValue;
+
+            // set the appropriate pixel value for this coordinate
+            if (sky) {
+                pixelValue = modelFunc (NULL, params, x);
+            } else {
+                pixelValue = modelFunc (NULL, params, x) - skyValue;
+            }
+
+
+            // add or subtract the value
+            if (add
+               ) {
+                image->data.F32[i][j] += pixelValue;
+            }
+            else {
+                image->data.F32[i][j] -= pixelValue;
+            }
         }
     }
@@ -1907,8 +1797,9 @@
                       psImage *mask,
                       pmModel *model,
-                      bool center)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 0);
+                      bool center,
+                      bool sky)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);
     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
@@ -1920,8 +1811,9 @@
                       psImage *mask,
                       pmModel *model,
-                      bool center)
-{
-    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
-    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, 1);
+                      bool center,
+                      bool sky)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);
     psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
Index: /branches/eam_rel9_p0/psModules/src/objects/pmObjects.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/pmObjects.h	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/pmObjects.h	(revision 5958)
@@ -10,6 +10,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.4.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-10 04:46:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -118,10 +118,14 @@
 pmPSFClump;
 
+// type of model carried by the pmModel structure
 typedef int pmModelType;
-#define PS_MODEL_GAUSS 0
-#define PS_MODEL_PGAUSS 1
-#define PS_MODEL_QGAUSS 2
-#define PS_MODEL_SGAUSS 3
-
+
+typedef enum {
+    PM_MODEL_UNTRIED,               ///< model fit not yet attempted
+    PM_MODEL_SUCCESS,               ///< model fit succeeded
+    PM_MODEL_NONCONVERGE,           ///< model fit did not converge
+    PM_MODEL_OFFIMAGE,              ///< model fit drove out of range
+    PM_MODEL_BADARGS                ///< model fit called with invalid args
+} pmModelStatus;
 
 /** pmModel data structure
@@ -142,4 +146,5 @@
     int nDOF;    ///< number of degrees of freedom
     int nIter;    ///< number of iterations to reach min
+    int status;         ///< fit status
     float radius;   ///< fit radius actually used
 }
@@ -157,23 +162,26 @@
  */
 typedef enum {
+    PM_SOURCE_UNKNOWN,                  ///< a cosmic-ray
     PM_SOURCE_DEFECT,                   ///< a cosmic-ray
     PM_SOURCE_SATURATED,                ///< random saturated pixels
-
-    PM_SOURCE_SATSTAR,                  ///< a saturated star
-    PM_SOURCE_PSFSTAR,                  ///< a PSF star
-    PM_SOURCE_GOODSTAR,                 ///< a good-quality star
-
-    PM_SOURCE_POOR_FIT_PSF,             ///< poor quality PSF fit
-    PM_SOURCE_FAIL_FIT_PSF,             ///< failed to get a good PSF fit
-    PM_SOURCE_FAINTSTAR,                ///< below S/N cutoff
-
+    PM_SOURCE_STAR,                     ///< a good-quality star
     PM_SOURCE_GALAXY,                   ///< an extended object (galaxy)
-    PM_SOURCE_FAINT_GALAXY,             ///< a galaxy below S/N cutoff
-    PM_SOURCE_DROP_GALAXY,              ///< ?
-    PM_SOURCE_FAIL_FIT_GAL,             ///< failed on the galaxy fit
-    PM_SOURCE_POOR_FIT_GAL,             ///< poor quality galaxy fit
-
-    PM_SOURCE_OTHER,                    ///< unidentified
 } pmSourceType;
+
+typedef enum {
+    PM_SOURCE_DEFAULT    = 0x0000, ///<
+    PM_SOURCE_PSFMODEL   = 0x0001, ///<
+    PM_SOURCE_FLTMODEL   = 0x0002, ///<
+    PM_SOURCE_SUBTRACTED = 0x0004, ///<
+    PM_SOURCE_FITTED     = 0x0008, ///<
+    PM_SOURCE_FAIL       = 0x0010, ///<
+    PM_SOURCE_POOR       = 0x0020, ///<
+    PM_SOURCE_PAIR       = 0x0040, ///<
+    PM_SOURCE_PSFSTAR    = 0x0080, ///<
+    PM_SOURCE_SATSTAR    = 0x0100, ///<
+    PM_SOURCE_BLEND      = 0x0200, ///<
+    PM_SOURCE_LINEAR     = 0x0400, ///<
+    PM_SOURCE_TEMPSUB    = 0x0800, ///< XXX get me a better name!
+} pmSourceMode;
 
 /** pmSource data structure
@@ -194,4 +202,8 @@
     pmModel *modelFLT;   ///< FLT (floating) Model fit (parameters and type).
     pmSourceType type;   ///< Best identification of object.
+    pmSourceMode mode;   ///< Best identification of object.
+    psArray *blends;
+    float apMag;
+    float fitMag;
 }
 pmSource;
@@ -507,5 +519,6 @@
     psImage *mask,   ///< The image pixel mask (valid == 0)
     pmModel *model,   ///< The input pmModel
-    bool center    ///< A boolean flag that determines whether pixels are centered
+    bool center,    ///< A boolean flag that determines whether pixels are centered
+    bool sky        ///< A boolean flag that determines if the sky is subtracted
 );
 
@@ -525,5 +538,6 @@
     psImage *mask,   ///< The image pixel mask (valid == 0)
     pmModel *model,   ///< The input pmModel
-    bool center    ///< A boolean flag that determines whether pixels are centered
+    bool center,    ///< A boolean flag that determines whether pixels are centered
+    bool sky        ///< A boolean flag that determines if the sky is subtracted
 );
 
Index: /branches/eam_rel9_p0/psModules/src/objects/pmPSF.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/pmPSF.c	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/pmPSF.c	(revision 5958)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-10 04:46:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -69,6 +69,6 @@
  X-center
  Y-center
- ???: Sky background value?
- ???: Zo?
+ Sky background value
+ Object Normalization
  *****************************************************************************/
 pmPSF *pmPSFAlloc (pmModelType type)
@@ -83,4 +83,24 @@
     psf->dApResid = 0.0;
     psf->skyBias  = 0.0;
+
+    // the ApTrend components are (x, y, rflux)
+    psf->ApTrend = psPolynomial3DAlloc (2, 2, 1, PS_POLYNOMIAL_ORD);
+
+    // we do not allow any rflux vs X,Y cross-terms
+    psf->ApTrend->mask[0][1][1] = 1;  // rflux x^0 y^1
+    psf->ApTrend->mask[0][2][1] = 1;  // rflux x^0 y^2
+    psf->ApTrend->mask[1][0][1] = 1;  // rflux x^1 y^0
+    psf->ApTrend->mask[1][1][1] = 1;  // rflux x^1 y^1
+    psf->ApTrend->mask[1][2][1] = 1;  // rflux x^1 y^2
+    psf->ApTrend->mask[2][0][1] = 1;  // rflux x^2 y^0
+    psf->ApTrend->mask[2][1][1] = 1;  // rflux x^2 y^1
+    psf->ApTrend->mask[2][2][1] = 1;  // rflux x^2 y^2
+
+    // only 2nd order terms, no such combinations
+    psf->ApTrend->mask[2][2][0] = 1;  // x^2 y^2
+    psf->ApTrend->mask[2][1][0] = 1;  // x^2 y^1
+    psf->ApTrend->mask[1][2][0] = 1;  // x^1 y^2
+
+    // total free parameters = 18 - 11 = 7
 
     Nparams = pmModelParameterCount (type);
Index: /branches/eam_rel9_p0/psModules/src/objects/pmPSF.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/pmPSF.h	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/pmPSF.h	(revision 5958)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-10-10 19:53:40 $
+ *  @version $Revision: 1.1.18.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-10 04:46:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -33,9 +33,11 @@
     pmModelType type;   ///< PSF Model in use
     psArray *params;   ///< Model parameters (psPolynomial2D)
+    psPolynomial3D *ApTrend;  ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst)
+    float ApResid;   ///< ???
+    float dApResid;   ///< ???
+    float skyBias;   ///< ???
     float chisq;   ///< PSF goodness statistic
-    float ApResid;                      ///< ???
-    float dApResid;                     ///< ???
-    float skyBias;                      ///< ???
     int nPSFstars;   ///< number of stars used to measure PSF
+    int nApResid;   ///< number of stars used to measure ApResid
 }
 pmPSF;
Index: /branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.c	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.c	(revision 5958)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-10 04:46:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -119,7 +119,4 @@
     psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
 
-    // make this optional?
-    // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat");
-
     // stage 2: construct a psf (pmPSF) from this collection of model fits
     pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask);
@@ -169,4 +166,6 @@
 
     }
+    psfTry->psf->nPSFstars = Npsf;
+
     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);
@@ -177,20 +176,11 @@
     // XXX this function wants aperture radius for pmSourcePhotometry
     pmPSFtryMetric (psfTry, RADIUS);
-    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
-              modelName,
-              psfTry->psf->ApResid,
-              psfTry->psf->dApResid,
-              psfTry->psf->skyBias);
-
+    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n",
+              modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias);
     return (psfTry);
 }
 
-
 bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
 {
-
-    float dBin;
-    int   nKeep, nSkip;
-
     // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
@@ -210,106 +200,24 @@
     }
 
-    // find min and max of (1/flux):
-    psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-    psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL);
-
-    // build binned versions of rflux, metric
-    dBin = (stats->max - stats->min) / 10.0;
-    psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64);
-    psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
-    psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
-    psFree (stats);
-
-    psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
-
-    // group data in daBin bins, measure lower 50% mean
-    for (int i = 0; i < daBin->n; i++) {
-
-        psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
-        tmp->n = 0;
-
-        // accumulate data within bin range
-        for (int j = 0; j < psfTry->sources->n; j++) {
-            // masked for: bad model fit, outlier in parameters
-            if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL)
-                continue;
-
-            // skip points with extreme dA values
-            if (fabs(psfTry->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] = psfTry->metric->data.F64[j];
-            tmp->n ++;
-        }
-
-        // is this a valid point?
-        maskB->data.U8[i] = 0;
-        if (tmp->n < 2) {
-            maskB->data.U8[i] = 1;
-            psFree (tmp);
-            continue;
-        }
-
-        // dA values are contaminated with low outliers
-        // measure statistics only on upper 50% of points
-        // this would be easier if we could sort in reverse:
-        //
-        // psVectorSort (tmp, tmp);
-        // tmp->n = 0.5*tmp->n;
-        // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-        // psVectorStats (stats, tmp, NULL, NULL, 0);
-        // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
-
-        psVectorSort (tmp, tmp);
-        nKeep = 0.5*tmp->n;
-        nSkip = tmp->n - nKeep;
-
-        psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
-        for (int j = 0; j < tmp2->n; j++) {
-            tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
-        }
-
-        stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-        psVectorStats (stats, tmp2, NULL, NULL, 0);
-        psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
-
-        daBin->data.F64[i] = stats->sampleMedian;
-
-        psFree (stats);
-        psFree (tmp);
-        psFree (tmp2);
-    }
-
-    // linear fit to rfBin, daBin
-    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
-
-    // XXX EAM : this is the intended API (cycle 7? cycle 8?)
-    poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
-
-    // XXX EAM : replace this when the above version is implemented
-    // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL);
-
-    psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
-    psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
-
-    stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
-    stats = psVectorStats (stats, daResid, NULL, maskB, 1);
-
-    psfTry->psf->ApResid = poly->coeff[0];
-    psfTry->psf->dApResid = stats->clippedStdev;
-    psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+    // use 3hi/1lo sigma clipping on the rflux vs metric fit
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    stats->min = 1.0;
+    stats->max = 3.0;
+    stats->clipIter = 3;
+
+    // fit ApResid only to rflux, ignore x,y variations for now
+    // linear clipped fit of ApResid to rflux
+    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+    poly = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, rflux);
+    psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+
+    psfTry->psf->ApTrend->coeff[0][0][0] = poly->coeff[0];
+    psfTry->psf->ApTrend->coeff[0][0][1] = 0;
+
+    psfTry->psf->skyBias  = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+    psfTry->psf->ApResid  = poly->coeff[0];
+    psfTry->psf->dApResid = stats->sampleStdev;
 
     psFree (rflux);
-    psFree (rfBin);
-    psFree (daBin);
-    psFree (maskB);
-    psFree (daBinFit);
-    psFree (daResid);
     psFree (poly);
     psFree (stats);
@@ -317,2 +225,9 @@
     return true;
 }
+
+/*
+  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
+  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
+  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
+*/
+
Index: /branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.h	(revision 5957)
+++ /branches/eam_rel9_p0/psModules/src/objects/pmPSFtry.h	(revision 5958)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 20:32:44 $
+ *  @version $Revision: 1.2.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-10 04:46:03 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -109,3 +109,15 @@
 );
 
+/** pmPSFtryMetric_Alt()
+ *
+ * This function is used to measure the PSF model metric for the set of
+ * results contained in the pmPSFtry structure (alternative implementation).
+ *
+ */
+bool pmPSFtryMetric_Alt(
+    pmPSFtry *try
+    ,                      ///< Add comment.
+    float RADIUS                        ///< Add comment.
+);
+
 # endif