Index: /branches/eam-psphot-branch/psModules/src/pmObjects.c
===================================================================
--- /branches/eam-psphot-branch/psModules/src/pmObjects.c	(revision 3789)
+++ /branches/eam-psphot-branch/psModules/src/pmObjects.c	(revision 3790)
@@ -5,6 +5,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-04-01 20:47:40 $
+ *  @version $Revision: 1.12.4.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-04-29 09:20:47 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -20,5 +20,5 @@
 /******************************************************************************
 pmPeakAlloc(): Allocate the psPeak data structure and set appropriate members.
- *****************************************************************************/
+*****************************************************************************/
 psPeak *pmPeakAlloc(psS32 x,
                     psS32 y,
@@ -38,5 +38,5 @@
 pmMomentsAlloc(): Allocate the psMoments structure and initialize the members
 to zero.
- *****************************************************************************/
+*****************************************************************************/
 psMoments *pmMomentsAlloc()
 {
@@ -64,5 +64,6 @@
 pmModelAlloc(): Allocate the psModel structure, along with its parameters,
 and initialize the type member.  Initialize the params to 0.0.
- *****************************************************************************/
+XXX EAM: changing params and dparams to psVector 
+*****************************************************************************/
 psModel *pmModelAlloc(psModelType type)
 {
@@ -73,32 +74,26 @@
     switch (type) {
     case PS_MODEL_GAUSS:
-        tmp->Nparams = 7;
-        tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
-        tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
+        tmp->params  = psVectorAlloc(7, PS_TYPE_F32);
+        tmp->dparams = psVectorAlloc(7, PS_TYPE_F32);
         break;
     case PS_MODEL_PGAUSS:
-        tmp->Nparams = 7;
-        tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
-        tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32));
+        tmp->params  = psVectorAlloc(7, PS_TYPE_F32);
+        tmp->dparams = psVectorAlloc(7, PS_TYPE_F32);
         break;
     case PS_MODEL_TWIST_GAUSS:
-        tmp->Nparams = 11;
-        tmp->params = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32));
-        tmp->dparams = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32));
+        tmp->params  = psVectorAlloc(11, PS_TYPE_F32);
+        tmp->dparams = psVectorAlloc(11, PS_TYPE_F32);
         break;
     case PS_MODEL_WAUSS:
-        tmp->Nparams = 9;
-        tmp->params = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32));
-        tmp->dparams = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32));
+        tmp->params  = psVectorAlloc(9, PS_TYPE_F32);
+        tmp->dparams = psVectorAlloc(9, PS_TYPE_F32);
         break;
     case PS_MODEL_SERSIC:
-        tmp->Nparams = 8;
-        tmp->params = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32));
-        tmp->dparams = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32));
+        tmp->params  = psVectorAlloc(8, PS_TYPE_F32);
+        tmp->dparams = psVectorAlloc(8, PS_TYPE_F32);
         break;
     case PS_MODEL_SERSIC_CORE:
-        tmp->Nparams = 12;
-        tmp->params = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32));
-        tmp->dparams = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32));
+        tmp->params  = psVectorAlloc(12, PS_TYPE_F32);
+        tmp->dparams = psVectorAlloc(12, PS_TYPE_F32);
         break;
     default:
@@ -107,7 +102,7 @@
     }
 
-    for (psS32 i = 0 ; i < tmp->Nparams ; i++) {
-        tmp->params[i] = 0.0;
-        tmp->dparams[i] = 0.0;
+    for (psS32 i = 0; i < tmp->params->n; i++) {
+        tmp->params->data.F32[i] = 0.0;
+        tmp->dparams->data.F32[i] = 0.0;
     }
 
@@ -119,5 +114,5 @@
 XXX: We don't free pixels and mask since that caused a memory error.
 We might need to increase the reference counter and decrease it here.
- *****************************************************************************/
+*****************************************************************************/
 static void p_psSourceFree(psSource *tmp)
 {
@@ -132,5 +127,5 @@
 pmSourceAlloc(): Allocate the psSource structure and initialize its members
 to NULL.
- *****************************************************************************/
+*****************************************************************************/
 psSource *pmSourceAlloc()
 {
@@ -141,4 +136,5 @@
     tmp->moments = NULL;
     tmp->models = NULL;
+    tmp->type = 0;
     p_psMemSetDeallocator(tmp, (psFreeFcn) p_psSourceFree);
 
@@ -156,5 +152,5 @@
 size of the output vector, then to set the values of the output vector.
 Depending upon actual use, this may need to be optimized.
- *****************************************************************************/
+*****************************************************************************/
 psVector *pmFindVectorPeaks(const psVector *vector,
                             psF32 threshold)
@@ -250,5 +246,5 @@
  
 XXX: Is there a better way to do this?
- *****************************************************************************/
+*****************************************************************************/
 psVector *p_psGetRowVectorFromImage(psImage *image,
                                     psU32 row)
@@ -265,22 +261,23 @@
 
 /******************************************************************************
-MyListAddPeak(): A private function which allocates a psList, if the list
+MyListAddPeak(): A private function which allocates a psArray, if the list
 argument is NULL, otherwise it adds the peak to that list.
- 
-XXX: Switch row, col args?
- *****************************************************************************/
-psList *MyListAddPeak(psList *list,
-                      psS32 row,
-                      psS32 col,
-                      psF32 counts,
-                      psPeakType type)
-{
-    psPeak *tmpPeak = pmPeakAlloc(row, col, counts, type);
+XXX EAM : changed the output to psArray 
+XXX EAM : Switched row, col args 
+XXX EAM : NOTE: this was changed in the call, so the new code is consistent
+*****************************************************************************/
+psArray *MyListAddPeak(psArray *list,
+                       psS32 row,
+                       psS32 col,
+                       psF32 counts,
+                       psPeakType type)
+{
+    psPeak *tmpPeak = pmPeakAlloc(col, row, counts, type);
 
     if (list == NULL) {
-        list = psListAlloc(tmpPeak);
-    } else {
-        psListAdd(list, PS_LIST_HEAD, tmpPeak);
-    }
+        list = psArrayAlloc(100);
+        list->n = 0;
+    }
+    psArrayAdd(list, 100, tmpPeak);
 
     return(list);
@@ -289,5 +286,5 @@
 /******************************************************************************
 pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
-above the given threshold.  Returns a psList containing location (x/y value)
+above the given threshold.  Returns a psArray containing location (x/y value)
 of all peaks.
  
@@ -299,9 +296,9 @@
 XXX: This does not work if image has either a single row, or a single column.
  
-XXX: In the output psList elements, should we use the image row/column offsets?
+XXX: In the output psArray elements, should we use the image row/column offsets?
      Currently, we do not.
- *****************************************************************************/
-psList *pmFindImagePeaks(const psImage *image,
-                         psF32 threshold)
+*****************************************************************************/
+psArray *pmFindImagePeaks(const psImage *image,
+                          psF32 threshold)
 {
     PS_IMAGE_CHECK_NULL(image, NULL);
@@ -313,5 +310,5 @@
     psU32 col = 0;
     psU32 row = 0;
-    psList *list = NULL;
+    psArray *list = NULL;
 
     //
@@ -321,4 +318,5 @@
     tmpRow = p_psGetRowVectorFromImage((psImage *) image, row);
     psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);
+
     for (psU32 i = 0 ; i < row1->n ; i++ ) {
         col = row1->data.U32[i];
@@ -359,4 +357,5 @@
         }
     }
+
     //
     // Exit if this image has a single row.
@@ -370,5 +369,5 @@
     //
     for (row = 1 ; row < (image->numRows - 1) ; row++) {
-        tmpRow = p_psGetRowVectorFromImage((psImage *) image, 0);
+        tmpRow = p_psGetRowVectorFromImage((psImage *) image, row);
         row1 = pmFindVectorPeaks(tmpRow, threshold);
 
@@ -481,5 +480,4 @@
         }
     }
-
     return(list);
 }
@@ -503,5 +501,5 @@
 
 /******************************************************************************
-psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psList that have
+psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have
 a peak value above the given maximum, or fall outside the valid region.
  
@@ -509,5 +507,7 @@
  
 XXX: warning message if valid is NULL?
- *****************************************************************************/
+ 
+XXX: changed API to create a NEW output psArray (should change name as well)
+*****************************************************************************/
 psList *pmCullPeaks(psList *peaks,
                     psF32 maxValue,
@@ -536,4 +536,28 @@
 }
 
+// XXX EAM: I changed this to return a new, subset array
+//          rather than alter the existing one
+psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid)
+{
+    PS_PTR_CHECK_NULL(peaks, NULL);
+
+    psArray *output = psArrayAlloc (200);
+    output->n = 0;
+
+    psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n);
+
+    for (int i = 0; i < peaks->n; i++) {
+        psPeak *tmpPeak = (psPeak *) peaks->data[i];
+        if (tmpPeak->counts > maxValue)
+            continue;
+        if (valid != NULL) {
+            if (IsItInThisRegion(valid, tmpPeak->x, tmpPeak->y))
+                continue;
+        }
+        psArrayAdd (output, 200, tmpPeak);
+    }
+    return(output);
+}
+
 /******************************************************************************
 psSource *pmSourceLocalSky(image, peak, innerRadius, outerRadius): this
@@ -572,5 +596,5 @@
 XXX: Don't use separate structs for the subimage and mask.  Use the source->
      members.
- *****************************************************************************/
+*****************************************************************************/
 psSource *pmSourceLocalSky(const psImage *image,
                            const psPeak *peak,
@@ -593,10 +617,13 @@
     // these variables should be renamed for clarity (imageCenterRow, etc).
     //
+    // peak->x,y is guaranteed to be on image
     psS32 SubImageCenterRow = peak->y;
     psS32 SubImageCenterCol = peak->x;
-    psS32 SubImageStartRow = SubImageCenterRow - outerRadiusS32;
-    psS32 SubImageEndRow = SubImageCenterRow + outerRadiusS32;
-    psS32 SubImageStartCol = SubImageCenterCol - outerRadiusS32;
-    psS32 SubImageEndCol = SubImageCenterCol + outerRadiusS32;
+
+    // XXX EAM : I added this code to stay on the image. So did George
+    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - outerRadiusS32);
+    psS32 SubImageEndRow    = PS_MIN (image->numRows - 1, SubImageCenterRow + outerRadiusS32);
+    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - outerRadiusS32);
+    psS32 SubImageEndCol    = PS_MIN (image->numCols - 1, SubImageCenterCol + outerRadiusS32);
     // AnulusWidth == number of pixels width in the annulus.  We add one since
     // the pixels at the inner AND outher radius are included.
@@ -609,9 +636,12 @@
     //    printf("pmSourceLocalSky(): AnulusWidth is %d\n", AnulusWidth);
 
-    if (SubImageStartRow < 0) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
-                SubImageStartRow);
-        return(NULL);
-    }
+    // XXX EAM : these tests should not be needed: we can never hit this error because of above
+    # if (1)
+
+        if (SubImageStartRow < 0) {
+            psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
+                    SubImageStartRow);
+            return(NULL);
+        }
     if (SubImageEndRow >= image->numRows) {
         psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
@@ -629,4 +659,5 @@
         return(NULL);
     }
+    # endif
 
     //
@@ -654,4 +685,5 @@
     //
     // Loop through the subimage, mask off pixels in the inner square.
+    // XXX this uses a static mask value of 1
     //
     for (psS32 row = AnulusWidth; row <= (subImageMask->numRows - AnulusWidth) - 1; row++) {
@@ -701,5 +733,5 @@
  
 XXX: macro this for performance.
-  *****************************************************************************/
+*****************************************************************************/
 bool CheckRadius(psPeak *peak,
                  psF32 radius,
@@ -719,5 +751,6 @@
  
 XXX: macro this for performance.
-  *****************************************************************************/
+XXX: this is rather inefficient - at least compute and compare against radius^2
+*****************************************************************************/
 bool CheckRadius2(psF32 xCenter,
                   psF32 yCenter,
@@ -726,4 +759,5 @@
                   psF32 y)
 {
+    /// XXX EAM should compare with hypot (x,y) for speed
     if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) {
         return(true);
@@ -746,5 +780,5 @@
  
 XXX: mask values?
- *****************************************************************************/
+*****************************************************************************/
 psSource *pmSourceMoments(psSource *source,
                           psF32 radius)
@@ -779,8 +813,14 @@
     psF32 Y2 = 0.0;
     psF32 XY = 0.0;
-    //
+    psF32 x = 0;
+    psF32 y = 0;
+    //
+    // XXX why do I get different results for these two methods of finding Sx?
+    // XXX Sx, Sy would be better measured if we clip pixels close to sky
+    // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
     // We loop through all pixels in this subimage (source->pixels), and for each
     // pixel that is not masked, AND within the radius of the peak pixel, we
-    // proceed with the moments calculation.
+    // proceed with the moments calculation.  need to do two loops for a
+    // numerically stable result.  first loop: get the sums.
     //
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -800,7 +840,8 @@
                     X1+= xDiff * pDiff;
                     Y1+= yDiff * pDiff;
+                    XY+= xDiff * yDiff * pDiff;
+
                     X2+= PS_SQR(xDiff) * pDiff;
                     Y2+= PS_SQR(yDiff) * pDiff;
-                    XY+= xDiff * yDiff * pDiff;
 
                     if (source->pixels->data.F32[row][col] > peakPixel) {
@@ -818,13 +859,84 @@
     // Sxy             = XY / Sum
     //
-    source->moments->x = X1/Sum + ((psF32) source->peak->x);
-    source->moments->y = Y1/Sum + ((psF32) source->peak->y);
-    source->moments->Sx = sqrt(X2/Sum - PS_SQR(X1/Sum));
-    source->moments->Sy = sqrt(Y2/Sum - PS_SQR(Y1/Sum));
+    x = X1/Sum;
+    y = Y1/Sum;
+    source->moments->x = x + ((psF32) source->peak->x);
+    source->moments->y = y + ((psF32) source->peak->y);
+
     source->moments->Sxy = XY/Sum;
+    source->moments->Sum = Sum;
     source->moments->Peak = peakPixel;
     source->moments->nPixels = numPixels;
 
+    // XXX EAM : these values can be negative, so we need to limit the range
+    source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
+    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
     return(source);
+
+    // XXX EAM : the following code should be the same as above, but it is not very stable: ignore it
+    # if (0)
+        //
+        // second loop: get the difference sums
+        //
+        X2 = Y2 = 0;
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+            if ((source->mask != NULL) && (source->mask->data.U8[row][col] != 0)) {
+                psS32 imgColCoord = col + source->pixels->col0;
+                psS32 imgRowCoord = row + source->pixels->row0;
+                if (CheckRadius(source->peak,
+                                radius,
+                                imgColCoord,
+                                imgRowCoord)) {
+                    psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
+                    psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
+                    psF32 pDiff = source->pixels->data.F32[row][col] - sky;
+
+                    Sum+= pDiff;
+                    X2+= PS_SQR(xDiff - x) * pDiff;
+                    Y2+= PS_SQR(yDiff - y) * pDiff;
+                }
+            }
+        }
+    }
+
+    //
+    // second moment X = sqrt (X2/Sum)
+    //
+    source->moments->Sx = (X2/Sum);
+    source->moments->Sy = (Y2/Sum);
+    return(source);
+    # endif
+}
+
+// XXX EAM : I used
+int pmComparePeakAscend (const void **a, const void **b)
+{
+    psPeak *A = *(psPeak **)a;
+    psPeak *B = *(psPeak **)b;
+
+    psF32 diff;
+
+    diff = A->counts - B->counts;
+    if (diff < FLT_EPSILON)
+        return (-1);
+    if (diff > FLT_EPSILON)
+        return (+1);
+    return (0);
+}
+
+int pmComparePeakDescend (const void **a, const void **b)
+{
+    psPeak *A = *(psPeak **)a;
+    psPeak *B = *(psPeak **)b;
+
+    psF32 diff;
+
+    diff = A->counts - B->counts;
+    if (diff < FLT_EPSILON)
+        return (+1);
+    if (diff > FLT_EPSILON)
+        return (-1);
+    return (0);
 }
 
@@ -839,73 +951,226 @@
  
 XXX: How can this function ever return FALSE?
- *****************************************************************************/
-#define SATURATE 0.0
-#define FAINT_SN_LIM 0.0
-#define PSF_SN_LIM 0.0
-#define SATURATE 0.0
-#define SATURATE 0.0
-
-bool pmSourceRoughClass(psArray *source,
-                        psMetadata *metadata)
-{
-    PS_PTR_CHECK_NULL(source, false);
+*****************************************************************************/
+
+# define NPIX 10
+# define SCALE 0.1
+
+// XXX I am ignore memory freeing issues (EAM)
+bool pmSourceRoughClass(psArray *sources, psMetadata *metadata)
+{
+    PS_PTR_CHECK_NULL(sources, false);
     PS_PTR_CHECK_NULL(metadata, false);
     psBool rc = true;
-
-    for (psS32 i = 0 ; i < source->n ; i++) {
-        psSource *tmpSrc = (psSource *) source->data[i];
-        PS_PTR_CHECK_NULL(tmpSrc->moments, false);
+    psArray *peaks  = NULL;
+    psF32 clumpX = 0.0;
+    psF32 clumpDX = 0.0;
+    psF32 clumpY = 0.0;
+    psF32 clumpDY = 0.0;
+
+    // find the sigmaX, sigmaY clump
+    {
+        psStats *stats  = NULL;
+        psImage *splane = NULL;
+        int binX, binY;
+
+        // construct a sigma-plane image
+        splane = psImageAlloc (NPIX/SCALE, NPIX/SCALE, PS_TYPE_F32);
+
+        // place the sources in the sigma-plane image (ignore 0,0 values?)
+        for (psS32 i = 0 ; i < sources->n ; i++)
+        {
+            psSource *tmpSrc = (psSource *) sources->data[i];
+            PS_PTR_CHECK_NULL(tmpSrc, false); // just skip this one?
+            PS_PTR_CHECK_NULL(tmpSrc->moments, false); // just skip this one?
+
+            // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
+            if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) {
+                continue;
+            }
+
+            // for the moment, force splane dimensions to be 10x10 image pix
+            binX = tmpSrc->moments->Sx/SCALE;
+            if (binX < 0)
+                continue;
+            if (binX >= splane->numCols)
+                continue;
+
+            binY = tmpSrc->moments->Sy/SCALE;
+            if (binY < 0)
+                continue;
+            if (binY >= splane->numRows)
+                continue;
+
+            splane->data.F32[binY][binX] += 1.0;
+        }
+
+        // find the peak in this image
+        stats = psStatsAlloc (PS_STAT_MAX);
+        stats = psImageStats (stats, splane, NULL, 0);
+        peaks = pmFindImagePeaks (splane, stats[0].max / 2);
+        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);
+
+    }
+
+    // measure statistics on Sx, Sy if Sx, Sy within range of clump
+    {
+        psPeak *clump;
+        psF32 minSx, maxSx;
+        psF32 minSy, maxSy;
+        psVector *tmpSx = NULL;
+        psVector *tmpSy = NULL;
+        psStats *stats  = NULL;
+
+        // XXX EAM : this lets us takes the single highest peak
+        psArraySort (peaks, pmComparePeakDescend);
+        clump = peaks->data[0];
+        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d\n", clump->x, clump->y);
+
+        // define section window for clump
+        minSx = clump->x * SCALE - 0.2;
+        maxSx = clump->x * SCALE + 0.2;
+        minSy = clump->y * SCALE - 0.2;
+        maxSy = clump->y * SCALE + 0.2;
+
+        tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32);
+        tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32);
+        tmpSx->n = 0;
+        tmpSy->n = 0;
+
+        // XXX clip sources based on flux?
+        // create vectors with Sx, Sy values in window
+        for (psS32 i = 0 ; i < sources->n ; i++)
+        {
+            psSource *tmpSrc = (psSource *) sources->data[i];
+
+            if (tmpSrc->moments->Sx < minSx)
+                continue;
+            if (tmpSrc->moments->Sx > maxSx)
+                continue;
+            if (tmpSrc->moments->Sy < minSy)
+                continue;
+            if (tmpSrc->moments->Sy > maxSy)
+                continue;
+            tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx;
+            tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy;
+            tmpSx->n++;
+            tmpSy->n++;
+            if (tmpSx->n == tmpSx->nalloc) {
+                psVectorRealloc (tmpSx, tmpSx->nalloc + 100);
+                psVectorRealloc (tmpSy, tmpSy->nalloc + 100);
+            }
+        }
+
+        // measures stats of Sx, Sy
+        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+
+        stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);
+        clumpX  = stats->clippedMean;
+        clumpDX = stats->clippedStdev;
+
+        stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);
+        clumpY  = stats->clippedMean;
+        clumpDY = stats->clippedStdev;
+
+        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", clumpX, clumpY);
+        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", clumpDX, clumpDY);
+        // these values should be pushed on the metadata somewhere
+    }
+
+    int Nsat   = 0;
+    int Ngal   = 0;
+    int Nfaint = 0;
+    int Nstar  = 0;
+    int Npsf   = 0;
+    int Ncr    = 0;
+    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
+    starsn->n = 0;
+
+    // XXX allow clump size to be scaled relative to sigmas?
+    // make rough IDs based on clumpX,Y,DX,DY
+    for (psS32 i = 0 ; i < sources->n ; i++) {
+
+        psSource *tmpSrc = (psSource *) sources->data[i];
+
         tmpSrc->peak->class = 0;
 
-        psF32 sigX = 0.0;
-        psF32 sigY = 0.0;
-        // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE.
-        psF32 clumpX = 0.0;
-        psF32 clumpDX = 0.0;
-        psF32 clumpY = 0.0;
-        psF32 clumpDY = 0.0;
-
+        psF32 sigX = tmpSrc->moments->Sx;
+        psF32 sigY = tmpSrc->moments->Sy;
+
+        // check return status value (do these exist?)
+        bool status;
+        psF32 RDNOISE  = psMetadataLookupF32 (&status, metadata, "RDNOISE");
+        psF32 GAIN     = psMetadataLookupF32 (&status, metadata, "GAIN");
+        psF32 SATURATE = psMetadataLookupF32 (&status, metadata, "SATURATE");
+
+        psF32 PSF_SN_LIM   = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
+        psF32 FAINT_SN_LIM = psMetadataLookupF32 (&status, metadata, "FAINT_SN_LIM");
+
+        // saturated object (star or single pixel not distinguished)
         if (tmpSrc->moments->Peak > SATURATE) {
-            tmpSrc->peak->class|= PS_SOURCE_SATURATED;
-        } else {
-            // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE.
-            psF32 gain = 0.0;
-            psF32 readNoise = 0.0;
-            psF32 S = tmpSrc->moments->Sum;
-            psF32 A = PS_PI * tmpSrc->moments->Sx * tmpSrc->moments->Sy;
-            psF32 B = tmpSrc->moments->Sky;
-            psF32 SN = (PS_SQRT_F32(gain) * S) /
-                       PS_SQRT_F32(S + (A * B) + ((A * readNoise * readNoise) / PS_SQRT_F32(gain)));
-            if (SN < FAINT_SN_LIM) {
-                tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR;
-            }
-            if (SN < PSF_SN_LIM) {
-                tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR;
-            }
-            // XXX: The SDRS is not real clear on how to calculate sigX, sigY.
-            if ((fabs(sigX - clumpX) < clumpDX) &&
-                    (fabs(sigY - clumpY) < clumpDY)) {
-                tmpSrc->peak->class|= PS_SOURCE_PSFSTAR;
-            }
-
-            if ((sigX < (clumpX - clumpDX)) &&
-                    (sigY < (clumpY - clumpDY)))
-                tmpSrc->peak->class|= PS_SOURCE_DEFECT;
-        }
-
-        if ((sigX > (clumpX + clumpDX)) &&
-                (sigY > (clumpY + clumpDY))) {
-            tmpSrc->peak->class|= PS_SOURCE_GALAXY;
-        }
-
-        if (tmpSrc->peak->class == 0) {
-            tmpSrc->peak->class|= PS_SOURCE_OTHER;
-        }
-    }
+            tmpSrc->type |= PS_SOURCE_SATURATED;
+            Nsat ++;
+            continue;
+        }
+
+        // too small to be stellar
+        if ((sigX < (clumpX - clumpDX)) || (sigY < (clumpY - clumpDY))) {
+            tmpSrc->type |= PS_SOURCE_DEFECT;
+            Ncr ++;
+            continue;
+        }
+
+        // possible galaxy
+        if ((sigX > (clumpX + clumpDX)) || (sigY > (clumpY + clumpDY))) {
+            tmpSrc->type |= PS_SOURCE_GALAXY;
+            Ngal ++;
+            continue;
+        }
+
+        // the rest are probable stellar objects
+        psF32 S  = tmpSrc->moments->Sum;
+        psF32 A  = PS_PI * sigX * sigY;
+        psF32 B  = tmpSrc->moments->Sky;
+        psF32 RT = PS_SQRT_F32(S + (A * B) + (A * PS_SQR(RDNOISE) / PS_SQRT_F32(GAIN)));
+        psF32 SN = (S * PS_SQRT_F32(GAIN) / RT);
+
+        starsn->data.F32[starsn->n] = SN;
+        starsn->n ++;
+        Nstar ++;
+
+        // faint star
+        if (SN < FAINT_SN_LIM) {
+            tmpSrc->type |= PS_SOURCE_FAINTSTAR;
+            Nfaint ++;
+            continue;
+        }
+
+        // PSF star
+        if (SN > PSF_SN_LIM) {
+            tmpSrc->type |= PS_SOURCE_PSFSTAR;
+            Npsf ++;
+            continue;
+        }
+
+        // random type of star
+        tmpSrc->type |= PS_SOURCE_OTHER;
+    }
+
+    {
+        psStats *stats  = NULL;
+        stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+        stats = psVectorStats (stats, starsn, NULL, NULL, 0);
+        psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
+    }
+
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:  %3d\n", Nstar);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:   %3d\n", Npsf);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nfaint: %3d\n", Nfaint);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ngal:   %3d\n", Ngal);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:   %3d\n", Nsat);
+    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:    %3d\n", Ncr);
 
     return(rc);
 }
-
-
 
 /******************************************************************************
@@ -919,5 +1184,5 @@
 XXX: The circle will have a diameter of (1+radius).  This is different from
      the pmSourceSetLocal() function.
- *****************************************************************************/
+*****************************************************************************/
 bool pmSourceSetPixelCircle(psSource *source,
                             const psImage *image,
@@ -927,6 +1192,6 @@
     PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
     PS_PTR_CHECK_NULL(source, false);
-    //    PS_PTR_CHECK_NULL(source->moments, false);
-    PS_PTR_CHECK_NULL(source->peak, false);
+    PS_PTR_CHECK_NULL(source->moments, false);
+    // PS_PTR_CHECK_NULL(source->peak, false);
     PS_FLOAT_COMPARE(0.0, radius, false);
 
@@ -940,15 +1205,19 @@
     psS32 SubImageCenterRow = source->peak->y;
     psS32 SubImageCenterCol = source->peak->x;
-    psS32 SubImageStartRow = SubImageCenterRow - radiusS32;
-    psS32 SubImageEndRow = SubImageCenterRow + radiusS32;
-    psS32 SubImageStartCol = SubImageCenterCol - radiusS32;
-    psS32 SubImageEndCol = SubImageCenterCol + radiusS32;
-
-    if (SubImageStartRow < 0) {
-        psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
-                SubImageStartRow);
-        return(false);
-    }
-    if (SubImageEndRow+1 >= image->numRows) {
+    // XXX EAM : for the circle to stay on the image
+    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - radiusS32);
+    psS32 SubImageEndRow    = PS_MIN (image->numRows - 1, SubImageCenterRow + radiusS32);
+    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
+    psS32 SubImageEndCol    = PS_MIN (image->numCols - 1, SubImageCenterCol + radiusS32);
+
+    // XXX EAM : this should not be needed: we can never hit this error
+    # if (1)
+
+        if (SubImageStartRow < 0) {
+            psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",
+                    SubImageStartRow);
+            return(false);
+        }
+    if (SubImageEndRow >= image->numRows) {
         psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",
                 SubImageEndRow);
@@ -960,13 +1229,16 @@
         return(false);
     }
-    if (SubImageEndCol+1 >= image->numCols) {
+    if (SubImageEndCol >= image->numCols) {
         psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",
                 SubImageEndCol);
         return(false);
     }
+    # endif
 
     // XXX: Must recycle image.
+    // XXX EAM: this message reflects a programming error we know about.
+    //          i am setting it to a trace message which we can take out
     if (source->pixels != NULL) {
-        psLogMsg(__func__, PS_LOG_WARN,
+        psTrace (".psModule.pmObjects.pmSourceSetPixelCircle", 4,
                  "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL.  Freeing and reallocating.\n");
         psFree(source->pixels);
@@ -975,6 +1247,6 @@
                                    SubImageStartCol,
                                    SubImageStartRow,
-                                   SubImageEndCol+1,
-                                   SubImageEndRow+1);
+                                   SubImageEndCol,
+                                   SubImageEndRow);
 
     // XXX: Must recycle image.
@@ -982,11 +1254,11 @@
         psFree(source->mask);
     }
-    source->mask = psImageAlloc(1 + 2 * radiusS32,
-                                1 + 2 * radiusS32,
+    source->mask = psImageAlloc(source->pixels->numCols,
+                                source->pixels->numRows,
                                 PS_TYPE_F32);
 
     //
     // Loop through the subimage mask, initialize mask to 0 or 1.
-    //
+    // XXX EAM: valid pixels should have 0, not 1
     for (psS32 row = 0 ; row < source->mask->numRows; row++) {
         for (psS32 col = 0 ; col < source->mask->numCols; col++) {
@@ -997,7 +1269,7 @@
                              (psF32) col,
                              (psF32) row)) {
+                source->mask->data.U8[row][col] = 0;
+            } else {
                 source->mask->data.U8[row][col] = 1;
-            } else {
-                source->mask->data.U8[row][col] = 0;
             }
         }
@@ -1005,5 +1277,4 @@
     return(true);
 }
-
 
 /******************************************************************************
@@ -1022,5 +1293,5 @@
 image, not subImage coords.  Therefore, the calls to the model evaluation
 functions will be in image, not subImage coords.  Remember this.
- *****************************************************************************/
+*****************************************************************************/
 bool pmSourceModelGuess(psSource *source,
                         const psImage *image,
@@ -1038,48 +1309,35 @@
     source->models = pmModelAlloc(model);
 
+    psVector *params = source->models->params;
+
     switch (model) {
     case PS_MODEL_GAUSS:
-        source->models->params[0] = source->moments->Sky;
-        source->models->params[1] = source->peak->counts - source->moments->Sky;
-        source->models->params[2] = source->moments->x;
-        source->models->params[3] = source->moments->y;
-        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
-        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
-        source->models->params[6] = source->moments->Sxy;
+        params->data.F32[0] = source->moments->Sky;
+        params->data.F32[1] = source->peak->counts - source->moments->Sky;
+        params->data.F32[2] = source->moments->x;
+        params->data.F32[3] = source->moments->y;
+        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
+        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
+        params->data.F32[6] = source->moments->Sxy;
         return(true);
 
     case PS_MODEL_PGAUSS:
-        source->models->params[0] = source->moments->Sky;
-        source->models->params[1] = source->peak->counts - source->moments->Sky;
-        source->models->params[2] = source->moments->x;
-        source->models->params[3] = source->moments->y;
-        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
-        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
-        source->models->params[6] = source->moments->Sxy;
+        params->data.F32[0] = source->moments->Sky;
+        params->data.F32[1] = source->peak->counts - source->moments->Sky;
+        params->data.F32[2] = source->moments->x;
+        params->data.F32[3] = source->moments->y;
+        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
+        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
+        params->data.F32[6] = source->moments->Sxy;
         return(true);
 
-    case PS_MODEL_TWIST_GAUSS:
-        source->models->params[0] = source->moments->Sky;
-        source->models->params[1] = source->peak->counts - source->moments->Sky;
-        source->models->params[2] = source->moments->x;
-        source->models->params[3] = source->moments->y;
-        // XXX: What are these?
-        // source->models->params[4] = SxInner;
-        // source->models->params[5] = SyInner;
-        // source->models->params[6] = SxyInner;
-        // source->models->params[7] = SxOuter;
-        // source->models->params[8] = SyOuter;
-        // source->models->params[9] = SxyOuter;
-        // source->models->params[10] = N;
-        return(true);
-
     case PS_MODEL_WAUSS:
-        source->models->params[0] = source->moments->Sky;
-        source->models->params[1] = source->peak->counts - source->moments->Sky;
-        source->models->params[2] = source->moments->x;
-        source->models->params[3] = source->moments->y;
-        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
-        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
-        source->models->params[6] = source->moments->Sxy;
+        params->data.F32[0] = source->moments->Sky;
+        params->data.F32[1] = source->peak->counts - source->moments->Sky;
+        params->data.F32[2] = source->moments->x;
+        params->data.F32[3] = source->moments->y;
+        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
+        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
+        params->data.F32[6] = source->moments->Sxy;
         // XXX: What are these?
         // source->models->params[7] = B2;
@@ -1087,32 +1345,47 @@
         return(true);
 
+        // XXX EAM : I might drop this model (or rather, replace it)
+    case PS_MODEL_TWIST_GAUSS:
+        params->data.F32[0] = source->moments->Sky;
+        params->data.F32[1] = source->peak->counts - source->moments->Sky;
+        params->data.F32[2] = source->moments->x;
+        params->data.F32[3] = source->moments->y;
+        // XXX: What are these?
+        // params->data.F32[4] = SxInner;
+        // params->data.F32[5] = SyInner;
+        // params->data.F32[6] = SxyInner;
+        // params->data.F32[7] = SxOuter;
+        // params->data.F32[8] = SyOuter;
+        // params->data.F32[9] = SxyOuter;
+        // params->data.F32[10] = N;
+        return(true);
+
     case PS_MODEL_SERSIC:
-        source->models->params[0] = source->moments->Sky;
-        source->models->params[1] = source->peak->counts - source->moments->Sky;
-        source->models->params[2] = source->moments->x;
-        source->models->params[3] = source->moments->y;
-        source->models->params[4] = sqrt(2.0) / source->moments->Sx;
-        source->models->params[5] = sqrt(2.0) / source->moments->Sy;
-        source->models->params[6] = source->moments->Sxy;
+        params->data.F32[0] = source->moments->Sky;
+        params->data.F32[1] = source->peak->counts - source->moments->Sky;
+        params->data.F32[2] = source->moments->x;
+        params->data.F32[3] = source->moments->y;
+        params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
+        params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
+        params->data.F32[6] = source->moments->Sxy;
         // XXX: What are these?
-        //source->models->params[7] = Nexp;
+        //params->data.F32[7] = Nexp;
         return(true);
 
     case PS_MODEL_SERSIC_CORE:
-        source->models->params[0] = source->moments->Sky;
-        source->models->params[1] = source->peak->counts - source->moments->Sky;
-        source->models->params[2] = source->moments->x;
-        source->models->params[3] = source->moments->y;
+        params->data.F32[0] = source->moments->Sky;
+        params->data.F32[1] = source->peak->counts - source->moments->Sky;
+        params->data.F32[2] = source->moments->x;
+        params->data.F32[3] = source->moments->y;
         // XXX: What are these?
-        //source->models->params[4] SxInner;
-        //source->models->params[5] SyInner;
-        //source->models->params[6] SxyInner;
-        //source->models->params[7] Zd;
-        //source->models->params[8] SxOuter;
-        //source->models->params[9] SyOuter;
-        //source->models->params[10] = SxyOuter;
-        //source->models->params[11] = Nexp;
+        // params->data.F32[4] SxInner;
+        // params->data.F32[5] SyInner;
+        // params->data.F32[6] SxyInner;
+        // params->data.F32[7] Zd;
+        // params->data.F32[8] SxOuter;
+        // params->data.F32[9] SyOuter;
+        // params->data.F32[10] = SxyOuter;
+        // params->data.F32[11] = Nexp;
         return(true);
-
     default:
         psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
@@ -1141,8 +1414,8 @@
 XXX: For a while, the first psVectorAlloc() was generating a seg fault during
 testing.  Try to reproduce that and debug.
- *****************************************************************************/
-psF32 evalModel(psSource *src,
-                psU32 row,
-                psU32 col)
+*****************************************************************************/
+static psF32 evalModel(psSource *src,
+                       psU32 row,
+                       psU32 col)
 {
     PS_PTR_CHECK_NULL(src, false);
@@ -1152,8 +1425,7 @@
     // XXX: The following step will not be necessary if the models->params
     // member is a psVector.  Suggest to IfA.
-    psVector *params = psVectorAlloc(7, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < src->models->Nparams ; i++) {
-        params->data.F32[i] = src->models->params[i];
-    }
+
+    // XXX EAM: done: models->params is now a vector
+    psVector *params = src->models->params;
 
     //
@@ -1188,6 +1460,4 @@
         return(NAN);
     }
-
-    psFree(params);
     psFree(x);
     return(tmpF);
@@ -1205,10 +1475,10 @@
  
 XXX: The result is returned in image coords.
- *****************************************************************************/
-psF32 findValue(psSource *source,
-                psF32 level,
-                psU32 row,
-                psU32 col,
-                psU32 dir)
+*****************************************************************************/
+static psF32 findValue(psSource *source,
+                       psF32 level,
+                       psU32 row,
+                       psU32 col,
+                       psU32 dir)
 {
     //
@@ -1280,5 +1550,5 @@
 XXX: What is momde?
 XXX: The top, bottom of the contour is not correctly determined.
- *****************************************************************************/
+*****************************************************************************/
 psArray *pmSourceContour(psSource *source,
                          const psImage *image,
@@ -1382,7 +1652,7 @@
 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x);
 
-//XXX: What should these values be?
-#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100
-#define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0
+// XXX EAM : these are better starting values, but should be available from metadata?
+#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 20
+#define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
 /******************************************************************************
 pmSourceFitModel(source, image): must create the appropiate arguments to the
@@ -1391,5 +1661,5 @@
 XXX: should there be a mask value?
 XXX: Probably should remove the "image" argument. 
- *****************************************************************************/
+*****************************************************************************/
 bool pmSourceFitModel(psSource *source,
                       const psImage *image)
@@ -1403,7 +1673,11 @@
     PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
     psBool rc;
+
+    // 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++) {
+    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++;
@@ -1411,16 +1685,23 @@
         }
     }
+
+    // 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++) {
+    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:
-                coord->data.F32[0] = (psF32) (i + source->pixels->row0);
-                coord->data.F32[1] = (psF32) (j + source->pixels->col0);
+                // 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];
+
+                // XXX EAM : this is approximate: need to apply the gain and rdnoise
+                yErr->data.F32[tmpCnt] = sqrt(PS_MAX(1, source->pixels->data.F32[i][j]));
                 tmpCnt++;
             }
@@ -1431,31 +1712,24 @@
                             PM_SOURCE_FIT_MODEL_TOLERANCE);
 
-    psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32);
+    psVector *params = source->models->params;
 
     switch (source->models->type) {
     case PS_MODEL_GAUSS:
-
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
-                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec);
+        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Gauss2D);
         break;
     case PS_MODEL_PGAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
-                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec);
+        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_PsuedoGauss2D);
         break;
     case PS_MODEL_TWIST_GAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
-                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec);
+        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Wauss2D);
         break;
     case PS_MODEL_WAUSS:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
-                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec);
+        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_TwistGauss2D);
         break;
     case PS_MODEL_SERSIC:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
-                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec);
+        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_Sersic);
         break;
     case PS_MODEL_SERSIC_CORE:
-        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y,
-                              NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec);
+        rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, yErr, pmMinLM_SersicCore);
         break;
     default:
@@ -1463,9 +1737,13 @@
         rc = false;
     }
+    // XXX EAM: we need to do something (give an error?) if rc is false
+    // XXX EAM: save the resulting chisq, nDOF, nIter
+    source->models->chisq = myMin->value;
+    source->models->nDOF  = y->n - params->n;
+    source->models->nIter = myMin->iter;
 
     psFree(x);
     psFree(y);
     psFree(myMin);
-    psFree(params);
     return(rc);
 }
@@ -1484,12 +1762,9 @@
     PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);
 
-    psVector *params = psVectorAlloc(src->models->Nparams, PS_TYPE_F32);
     psVector *x = psVectorAlloc(2, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < src->models->Nparams ; i++) {
-        params->data.F32[i] = src->models->params[i];
-    }
-
-    for (psS32 i = 0 ; i < src->pixels->numRows ; i++) {
-        for (psS32 j = 0 ; j < src->pixels->numCols ; j++) {
+    psVector *params = src->models->params;
+
+    for (psS32 i = 0; i < src->pixels->numRows; i++) {
+        for (psS32 j = 0; j < src->pixels->numCols; j++) {
             psF32 pixelValue;
             // XXX: Should you be adding the pixels for the entire subImage,
@@ -1525,5 +1800,4 @@
                 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
                 psFree(x);
-                psFree(params);
                 return(false);
             }
@@ -1539,5 +1813,4 @@
     }
     psFree(x);
-    psFree(params);
     return(true);
 }
@@ -1574,10 +1847,14 @@
 
 
-/******************************************************************************
-pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x,
-y" coordinate pair.  This function computes the gaussian, specified by the
-parameters in "params" at that x,y point and returns the value.  The
-derivatives are also caculated and returned in the "deriv" argument.
- 
+/**
+   all of these object representation functions have the same form : func(*deriv, *params, *x)
+ 
+   the argument "x" contains a single "x,y" coordinate pair.  The function computes the object
+   model, based on the parameters in "params" at the x,y point specified by *x, and returns the value.  
+   The derivatives are also caculated and returned in the "deriv" argument.  parameter error checking is 
+   skipped because speed is most important.
+**/
+
+/******************************************************************************
     params->data.F32[0] = So;
     params->data.F32[1] = Zo;
@@ -1587,39 +1864,19 @@
     params->data.F32[5] = sqrt(2.0) / SigmaY;
     params->data.F32[6] = Sxy;
- 
-XXX: Consider getting rid of the parameter checks since this might consume
-a significant fraction of this function CPU time.
- 
-XXX: I added the following.  Must conform with IfA.  If deriv==NULL, then
-we simply don't perform the derivative calculations.
- 
-XXX: It is not clear whether the subImage coords, or the image coords should
-be used when calling these functions.  I don't think it really matters, as
-long as we are consistent.
- *****************************************************************************/
-psF32 pmMinLM_Gauss2D(psVector *deriv,
+*****************************************************************************/
+psF64 pmMinLM_Gauss2D(psVector *deriv,
                       psVector *params,
                       psVector *x)
 {
-    PS_VECTOR_CHECK_NULL(params, NAN);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(params, 7, NAN);
-    PS_VECTOR_CHECK_NULL(x, NAN);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
-
-    psF32 X = x->data.F32[0] - params->data.F32[2];
-    psF32 Y = x->data.F32[1] - params->data.F32[3];
+    psF32 X  = x->data.F32[0] - params->data.F32[2];
+    psF32 Y  = x->data.F32[1] - params->data.F32[3];
     psF32 px = params->data.F32[4]*X;
     psF32 py = params->data.F32[5]*Y;
-    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 r = exp(-z);
-    psF32 f = params->data.F32[1]*r + params->data.F32[0];
+    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
+    psF32 r  = exp(-z);
+    psF32 q  = params->data.F32[1]*r;
+    psF32 f  = q + params->data.F32[0];
 
     if (deriv != NULL) {
-        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
-        PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);
-
-        psF32 q = params->data.F32[1]*r;
         deriv->data.F32[0] = +1.0;
         deriv->data.F32[1] = +r;
@@ -1630,49 +1887,6 @@
         deriv->data.F32[6] = -q*X*Y;
     }
-
     return(f);
 }
-
-/******************************************************************************
-p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above
-function in a form that is usable in the LM minimization routines.
- *****************************************************************************/
-psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv,
-                                psVector *params,
-                                psArray *x)
-{
-    PS_IMAGE_CHECK_NULL(deriv, NULL);
-    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
-    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
-    PS_VECTOR_CHECK_NULL(params, NULL);
-    PS_VECTOR_CHECK_EMPTY(params, NULL);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
-    PS_PTR_CHECK_NULL(x, NULL);
-    if (deriv->numRows != x->n) {
-        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
-    }
-    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
-    // XXX: use static memory here.
-    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
-
-    for (psS32 i = 0 ; i < x->n ; i++) {
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            tmpRow->data.F32[j] = deriv->data.F32[i][j];
-        }
-
-        psVector *tmpVec2 = (psVector *) x->data[i];
-        tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow,
-                                              params,
-                                              tmpVec2);
-
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            deriv->data.F32[i][j] = tmpRow->data.F32[j];
-        }
-    }
-
-    psFree(tmpRow);
-    return(tmpVec);
-}
-
 
 /******************************************************************************
@@ -1684,88 +1898,31 @@
     params->data.F32[5] = sqrt(2) / SigmaY;
     params->data.F32[6] = Sxy;
- *****************************************************************************/
-psF32 pmMinLM_PsuedoGauss2D(psVector *deriv,
+*****************************************************************************/
+psF64 pmMinLM_PsuedoGauss2D(psVector *deriv,
                             psVector *params,
                             psVector *x)
 {
-    PS_VECTOR_CHECK_NULL(params, NAN);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(params, 7, NAN);
-    PS_VECTOR_CHECK_NULL(x, NAN);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
-
-    psF32 X = x->data.F32[0] - params->data.F32[2];
-    psF32 Y = x->data.F32[1] - params->data.F32[3];
+    psF32 X  = x->data.F32[0] - params->data.F32[2];
+    psF32 Y  = x->data.F32[1] - params->data.F32[3];
     psF32 px = params->data.F32[4]*X;
     psF32 py = params->data.F32[5]*Y;
-    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 t = 1 + z + 0.5*z*z;
-    psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
-    psF32 f = params->data.F32[1]*r + params->data.F32[0];
+    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
+    psF32 t  = 1 + z + 0.5*z*z;
+    psF32 r  = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
+    psF32 f  = params->data.F32[1]*r + params->data.F32[0];
 
     if (deriv != NULL) {
-        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
-        PS_VECTOR_CHECK_SIZE(deriv, 7, NAN);
-
-        //
         // note difference from a pure gaussian: q = params->data.F32[1]*r
-        //
         psF32 q = params->data.F32[1]*r*r*t;
         deriv->data.F32[0] = +1.0;
         deriv->data.F32[1] = +r;
         deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
-        deriv->data.F32[3] = q *
-                             (2.0*py*params->data.F32[5] + params->data.F32[6]*X);
+        deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
         deriv->data.F32[4] = -2.0*q*px*X;
         deriv->data.F32[5] = -2.0*q*py*Y;
         deriv->data.F32[6] = -q*X*Y;
     }
-
     return(f);
 }
-
-/******************************************************************************
-p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the
-above function in a form that is usable in the LM minimization routines.
- *****************************************************************************/
-psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv,
-                                      psVector *params,
-                                      psArray *x)
-{
-    PS_IMAGE_CHECK_NULL(deriv, NULL);
-    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
-    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
-    PS_VECTOR_CHECK_NULL(params, NULL);
-    PS_VECTOR_CHECK_EMPTY(params, NULL);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
-    PS_PTR_CHECK_NULL(x, NULL);
-    if (deriv->numRows != x->n) {
-        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
-    }
-    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
-    // XXX: use static memory here.
-    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
-
-    for (psS32 i = 0 ; i < x->n ; i++) {
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            tmpRow->data.F32[j] = deriv->data.F32[i][j];
-        }
-
-        tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow,
-                              params,
-                              (psVector *) x->data[i]);
-
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            deriv->data.F32[i][j] = tmpRow->data.F32[j];
-        }
-    }
-
-    psFree(tmpRow);
-    return(tmpVec);
-}
-
-
-
 
 /******************************************************************************
@@ -1779,16 +1936,9 @@
     params->data.F32[7] = B2;
     params->data.F32[8] = B3;
- *****************************************************************************/
-psF32 pmMinLM_Wauss2D(psVector *deriv,
+*****************************************************************************/
+psF64 pmMinLM_Wauss2D(psVector *deriv,
                       psVector *params,
                       psVector *x)
 {
-    PS_VECTOR_CHECK_NULL(params, NAN);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(params, 9, NAN);
-    PS_VECTOR_CHECK_NULL(x, NAN);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
-
     psF32 X = x->data.F32[0] - params->data.F32[2];
     psF32 Y = x->data.F32[1] - params->data.F32[2];
@@ -1801,10 +1951,5 @@
 
     if (deriv != NULL) {
-        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
-        PS_VECTOR_CHECK_SIZE(deriv, 9, NAN);
-        //
         // note difference from gaussian: q = params->data.F32[1]*r
-        //
-
         psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0));
         deriv->data.F32[0] = +1.0;
@@ -1817,56 +1962,8 @@
         deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t;
         deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0;
-        //
         // The values of 100 dampen the swing of params->data.F32[7,8] */
-        //
-    }
-
+    }
     return(f);
 }
-
-/******************************************************************************
-p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above
-function in a form that is usable in the LM minimization routines.
- *****************************************************************************/
-psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv,
-                                psVector *params,
-                                psArray *x)
-{
-    PS_IMAGE_CHECK_NULL(deriv, NULL);
-    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
-    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
-    PS_VECTOR_CHECK_NULL(params, NULL);
-    PS_VECTOR_CHECK_EMPTY(params, NULL);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
-    PS_PTR_CHECK_NULL(x, NULL);
-    if (deriv->numRows != x->n) {
-        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
-    }
-    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
-    // XXX: use static memory here.
-    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
-
-    for (psS32 i = 0 ; i < x->n ; i++) {
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            tmpRow->data.F32[j] = deriv->data.F32[i][j];
-        }
-
-        tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow,
-                                              params,
-                                              (psVector *) x->data[i]);
-
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            deriv->data.F32[i][j] = tmpRow->data.F32[j];
-        }
-    }
-
-    psFree(tmpRow);
-    return(tmpVec);
-}
-
-
-
-
-
 
 // XXX: What should these be?
@@ -1885,16 +1982,9 @@
     params->data.F32[9] = SxyOuter;
     params->data.F32[10] = N;
- *****************************************************************************/
-psF32 pmMinLM_TwistGauss2D(psVector *deriv,
+*****************************************************************************/
+psF64 pmMinLM_TwistGauss2D(psVector *deriv,
                            psVector *params,
                            psVector *x)
 {
-    PS_VECTOR_CHECK_NULL(params, NAN);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(params, 11, NAN);
-    PS_VECTOR_CHECK_NULL(x, NAN);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
-
     psF32 X = x->data.F32[0] - params->data.F32[2];
     psF32 Y = x->data.F32[1] - params->data.F32[3];
@@ -1907,11 +1997,7 @@
     psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10]));
 
-
     psF32 f = params->data.F32[5]*r + params->data.F32[6];
 
     if (deriv != NULL) {
-        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
-        PS_VECTOR_CHECK_SIZE(deriv, 11, NAN);
-
         psF32 q1 = params->data.F32[5]*PS_SQR(r);
         psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0));
@@ -1921,16 +2007,12 @@
         deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X);
 
-        //
         // These fudge factors impede the growth of params->data.F32[4] beyond
         // params->data.F32[7].
-        //
         psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]);
         psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
         deriv->data.F32[4] = -2.0*q1*px1*X*f2;
 
-        //
         // These fudge factors impede the growth of params->data.F32[5] beyond
         // params->data.F32[8].
-        //
         f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]);
         f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
@@ -1945,46 +2027,4 @@
     return(f);
 }
-
-/******************************************************************************
-p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above
-function in a form that is usable in the LM minimization routines.
- *****************************************************************************/
-psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv,
-                                     psVector *params,
-                                     psArray *x)
-{
-    PS_IMAGE_CHECK_NULL(deriv, NULL);
-    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
-    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
-    PS_VECTOR_CHECK_NULL(params, NULL);
-    PS_VECTOR_CHECK_EMPTY(params, NULL);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
-    PS_PTR_CHECK_NULL(x, NULL);
-    if (deriv->numRows != x->n) {
-        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
-    }
-    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
-    // XXX: use static memory here.
-    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
-
-    for (psS32 i = 0 ; i < x->n ; i++) {
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            tmpRow->data.F32[j] = deriv->data.F32[i][j];
-        }
-
-        tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow,
-                              params,
-                              (psVector *) x->data[i]);
-
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            deriv->data.F32[i][j] = tmpRow->data.F32[j];
-        }
-    }
-
-    psFree(tmpRow);
-    return(tmpVec);
-}
-
-
 
 /******************************************************************************
@@ -1998,62 +2038,11 @@
     params->data.F32[6] = Sxy;
     params->data.F32[7] = Nexp;
- *****************************************************************************/
-psF32 pmMinLM_Sersic(psVector *deriv,
+*****************************************************************************/
+psF64 pmMinLM_Sersic(psVector *deriv,
                      psVector *params,
                      psVector *x)
 {
-    PS_VECTOR_CHECK_NULL(params, NAN);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(params, 8, NAN);
-    PS_VECTOR_CHECK_NULL(x, NAN);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
-
-    if (deriv != NULL) {
-        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
-        PS_VECTOR_CHECK_SIZE(deriv, 8, NAN);
-    }
-
     psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
     return(0.0);
-}
-/******************************************************************************
-p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above
-function in a form that is usable in the LM minimization routines.
- *****************************************************************************/
-psVector *p_pmMinLM_Sersic_Vec(psImage *deriv,
-                               psVector *params,
-                               psArray *x)
-{
-    PS_IMAGE_CHECK_NULL(deriv, NULL);
-    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
-    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
-    PS_VECTOR_CHECK_NULL(params, NULL);
-    PS_VECTOR_CHECK_EMPTY(params, NULL);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
-    PS_PTR_CHECK_NULL(x, NULL);
-    if (deriv->numRows != x->n) {
-        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
-    }
-    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
-    // XXX: use static memory here.
-    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
-
-    for (psS32 i = 0 ; i < x->n ; i++) {
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            tmpRow->data.F32[j] = deriv->data.F32[i][j];
-        }
-
-        tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow,
-                                             params,
-                                             (psVector *) x->data[i]);
-
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            deriv->data.F32[i][j] = tmpRow->data.F32[j];
-        }
-    }
-
-    psFree(tmpRow);
-    return(tmpVec);
 }
 
@@ -2072,72 +2061,10 @@
     params->data.F32[10] = SxyOuter;
     params->data.F32[11] = Nexp;
- *****************************************************************************/
-psF32 pmMinLM_SersicCore(psVector *deriv,
+*****************************************************************************/
+psF64 pmMinLM_SersicCore(psVector *deriv,
                          psVector *params,
                          psVector *x)
 {
-    PS_VECTOR_CHECK_NULL(params, NAN);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(params, 12, NAN);
-    PS_VECTOR_CHECK_NULL(x, NAN);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN);
-    PS_VECTOR_CHECK_SIZE(x, 2, NAN);
-
-    if (deriv != NULL) {
-        PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN);
-        PS_VECTOR_CHECK_SIZE(deriv, 12, NAN);
-    }
-
     psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet.");
     return(0.0);
 }
-/******************************************************************************
-p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above
-function in a form that is usable in the LM minimization routines.
- *****************************************************************************/
-psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv,
-                                   psVector *params,
-                                   psArray *x)
-{
-    PS_IMAGE_CHECK_NULL(deriv, NULL);
-    PS_IMAGE_CHECK_EMPTY(deriv, NULL);
-    PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL);
-    PS_VECTOR_CHECK_NULL(params, NULL);
-    PS_VECTOR_CHECK_EMPTY(params, NULL);
-    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL);
-    PS_PTR_CHECK_NULL(x, NULL);
-    if (deriv->numRows != x->n) {
-        psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x.");
-    }
-    psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32);
-    // XXX: use static memory here.
-    psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32);
-
-    for (psS32 i = 0 ; i < x->n ; i++) {
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            tmpRow->data.F32[j] = deriv->data.F32[i][j];
-        }
-
-        tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow,
-                              params,
-                              (psVector *) x->data[i]);
-
-        for (psS32 j = 0 ; j < tmpRow->n ; j++) {
-            deriv->data.F32[i][j] = tmpRow->data.F32[j];
-        }
-    }
-
-    psFree(tmpRow);
-    return(tmpVec);
-}
-
-
-
-/******************************************************************************
- *****************************************************************************/
-psF32 pmMinLM_PsuedoSersic(psVector *deriv,
-                           psVector *params,
-                           psVector *x)
-{
-    return(0.0);
-}
Index: /branches/eam-psphot-branch/psModules/src/pmObjects.h
===================================================================
--- /branches/eam-psphot-branch/psModules/src/pmObjects.h	(revision 3789)
+++ /branches/eam-psphot-branch/psModules/src/pmObjects.h	(revision 3790)
@@ -5,6 +5,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.7.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-04-06 19:34:06 $
+ *  @version $Revision: 1.7.2.1.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-04-29 09:20:47 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -54,5 +54,5 @@
 
 typedef enum {
-    PS_MODEL_GAUSS,  ///< Regular 2-D Gaussian
+    PS_MODEL_GAUSS,   ///< Regular 2-D Gaussian
     PS_MODEL_PGAUSS,  ///< Psuedo 2-D Gaussian
     PS_MODEL_TWIST_GAUSS, ///< 2-D Twisted Gaussian
@@ -67,8 +67,9 @@
 {
     psModelType type;       ///< Model to be used.
-    psS32 Nparams;          ///< Number of parameters.
-    psF32 *params;          ///< Paramater values.
-    psF32 *dparams;         ///< Parameter errors.
+    psVector *params;       ///< Paramater values.
+    psVector *dparams;      ///< Parameter errors.
     psF32 chisq;            ///< Fit chi-squared.
+    psS32 nDOF;             ///< number of degrees of freedom
+    psS32 nIter;            ///< number of iterations to reach min
 }
 psModel;
@@ -121,7 +122,7 @@
 value) of all peaks.
  *****************************************************************************/
-psList *pmFindImagePeaks(const psImage *image, ///< The input image where peaks will be found (psF32)
-                         psF32 threshold ///< Threshold above which to find a peak
-                        );
+psArray *pmFindImagePeaks(const psImage *image, ///< The input image where peaks will be found (psF32)
+                          psF32 threshold ///< Threshold above which to find a peak
+                         );
 
 /******************************************************************************
@@ -209,6 +210,7 @@
 /******************************************************************************
 XXX: Why only *x argument?
- *****************************************************************************/
-psF32 pmMinLM_Gauss2D(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
+XXX EAM: psMinimizeLMChi2Func returns psF64, not psF32
+ *****************************************************************************/
+psF64 pmMinLM_Gauss2D(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
                       psVector *params,  ///< A psVector which holds the parameters of this function
                       psVector *x  ///< A psVector which holds the row/col coordinate
@@ -218,5 +220,5 @@
 XXX: Why only *x argument?
  *****************************************************************************/
-psF32 pmMinLM_PsuedoGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
+psF64 pmMinLM_PsuedoGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
                             psVector *params, ///< A psVector which holds the parameters of this function
                             psVector *x  ///< A psVector which holds the row/col coordinate
@@ -225,5 +227,5 @@
 /******************************************************************************
  *****************************************************************************/
-psF32 pmMinLM_Wauss2D(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
+psF64 pmMinLM_Wauss2D(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
                       psVector *params,  ///< A psVector which holds the parameters of this function
                       psVector *x  ///< A psVector which holds the row/col coordinate
@@ -232,5 +234,5 @@
 /******************************************************************************
  *****************************************************************************/
-psF32 pmMinLM_TwistGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
+psF64 pmMinLM_TwistGauss2D(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
                            psVector *params, ///< A psVector which holds the parameters of this function
                            psVector *x  ///< A psVector which holds the row/col coordinate
@@ -239,5 +241,5 @@
 /******************************************************************************
  *****************************************************************************/
-psF32 pmMinLM_Sersic(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
+psF64 pmMinLM_Sersic(psVector *deriv,  ///< A possibly-NULL structure for the output derivatives
                      psVector *params,  ///< A psVector which holds the parameters of this function
                      psVector *x  ///< A psVector which holds the row/col coordinate
@@ -246,5 +248,5 @@
 /******************************************************************************
  *****************************************************************************/
-psF32 pmMinLM_SersicCore(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
+psF64 pmMinLM_SersicCore(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
                          psVector *params, ///< A psVector which holds the parameters of this function
                          psVector *x  ///< A psVector which holds the row/col coordinate
@@ -253,5 +255,5 @@
 /******************************************************************************
  *****************************************************************************/
-psF32 pmMinLM_PsuedoSersic(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
+psF64 pmMinLM_PsuedoSersic(psVector *deriv, ///< A possibly-NULL structure for the output derivatives
                            psVector *params, ///< A psVector which holds the parameters of this function
                            psVector *x  ///< A psVector which holds the row/col coordinate