Index: trunk/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- trunk/psLib/src/math/psMinimizePolyFit.c	(revision 6185)
+++ trunk/psLib/src/math/psMinimizePolyFit.c	(revision 6186)
@@ -10,16 +10,10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-23 20:44:29 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-23 22:25:31 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  *
- *  XXX: must follow coding name standards on local functions.
- *  XXX: put local functions in front.
  *  XXX: For clip-fit functions, what should we do if the mask is NULL?
- *  XXX: Generate macros for
- * PS_ASSERT_VECTOR_TYPES_MATCH()
- * PS_ASSERT_VECTOR_SIZES_MATCH()
- *  They should have better error messages.
  *
  */
@@ -42,4 +36,16 @@
 /*****************************************************************************/
 
+#define PS_VECTOR_GEN_CHEBY_INDEX(VEC, SIZE, TYPE) \
+VEC = psVectorAlloc(SIZE, TYPE); \
+if (TYPE == PS_TYPE_F64) { \
+    for (psS32 i = 0 ; i < SIZE ; i++) { \
+        VEC->data.F64[i] = ((2.0 / ((psF64) (SIZE - 1))) * ((psF64) i)) - 1.0; \
+    }\
+} else if (TYPE == PS_TYPE_F32){ \
+    for (psS32 i = 0 ; i < SIZE ; i++) { \
+        VEC->data.F32[i] = ((2.0 / ((psF32) (SIZE - 1))) * ((psF32) i)) - 1.0; \
+    }\
+}\
+
 /*****************************************************************************/
 /* TYPE DEFINITIONS                                                          */
@@ -58,34 +64,4 @@
 /*****************************************************************************/
 
-/******************************************************************************
- ******************************************************************************
- Analytical 1-D fitting routines.
- ******************************************************************************
- *****************************************************************************/
-// XXX: Make this a general type conversion macro, or function
-#define PS_VECTOR_GEN_F64_FROM_F32(VECF64, VECF32) \
-VECF64 = psVectorAlloc(VECF32->n, PS_TYPE_F64); \
-for (psS32 i = 0 ; i < VECF32->n ; i++) { \
-    VECF64->data.F64[i] = (psF64) VECF32->data.F32[i]; \
-} \
-
-// XXX: Consolidate these
-#define PS_VECTOR_GEN_CHEBY_INDEX(VECF64, SIZE) \
-VECF64 = psVectorAlloc(SIZE, PS_TYPE_F64); \
-for (psS32 i = 0 ; i < SIZE ; i++) { \
-    VECF64->data.F64[i] = ((2.0 / ((psF64) (SIZE - 1))) * ((psF64) i)) - 1.0; \
-}\
-
-#define PS_VECTOR_GEN_CHEBY_INDEX_F64(VECF64, SIZE) \
-VECF64 = psVectorAlloc(SIZE, PS_TYPE_F64); \
-for (psS32 i = 0 ; i < SIZE ; i++) { \
-    VECF64->data.F64[i] = ((2.0 / ((psF64) (SIZE - 1))) * ((psF64) i)) - 1.0; \
-}\
-
-#define PS_VECTOR_GEN_CHEBY_INDEX_F32(VECF32, SIZE) \
-VECF32 = psVectorAlloc(SIZE, PS_TYPE_F32); \
-for (psS32 i = 0 ; i < SIZE ; i++) { \
-    VECF32->data.F32[i] = ((2.0 / ((psF32) (SIZE - 1))) * ((psF32) i)) - 1.0; \
-}\
 /******************************************************************************
 BuildSums1D(sums, x, polyOrder, sums): this routine calculates the powers of
@@ -113,5 +89,5 @@
 
     xSum = 1.0;
-    for (int i = 0; i < nSum; i++) {
+    for (psS32 i = 0; i < nSum; i++) {
         sums->data.F64[i] = xSum;
         xSum *= x;
@@ -148,21 +124,11 @@
 
     xSum = 1.0;
-    for (int i = 0; i < nXsum; i++) {
+    for (psS32 i = 0; i < nXsum; i++) {
         ySum = xSum;
-        for (int j = 0; j < nYsum; j++) {
+        for (psS32 j = 0; j < nYsum; j++) {
             sums->data.F64[i][j] = ySum;
             ySum *= y;
         }
         xSum *= x;
-    }
-
-    if (0) {
-        printf("--------------------- BuildSums2D(%.2f %.2f) ---------------------\n", x, y);
-        for (int i = 0; i < nXsum; i++) {
-            for (int j = 0; j < nYsum; j++) {
-                printf("(%.2f) ", sums->data.F64[i][j]);
-            }
-            printf("\n");
-        }
     }
 
@@ -196,7 +162,7 @@
     if (sums == NULL) {
         sums = (psF64 ***) psAlloc (nXsum*sizeof(psF64));
-        for (int i = 0; i < nXsum; i++) {
+        for (psS32 i = 0; i < nXsum; i++) {
             sums[i] = (psF64 **) psAlloc (nYsum*sizeof(psF64));
-            for (int j = 0; j < nYsum; j++) {
+            for (psS32 j = 0; j < nYsum; j++) {
                 sums[i][j] = (psF64 *) psAlloc (nZsum*sizeof(psF64));
             }
@@ -207,9 +173,9 @@
     if (1) {
         zSum = 1.0;
-        for (int k = 0; k < nZsum; k++) {
+        for (psS32 k = 0; k < nZsum; k++) {
             ySum = zSum;
-            for (int j = 0; j < nYsum; j++) {
+            for (psS32 j = 0; j < nYsum; j++) {
                 xSum = ySum;
-                for (int i = 0; i < nXsum; i++) {
+                for (psS32 i = 0; i < nXsum; i++) {
                     sums[i][j][k] = xSum;
                     xSum *= x;
@@ -221,9 +187,9 @@
     } else {
         xSum = 1.0;
-        for (int i = 0; i < nXsum; i++) {
+        for (psS32 i = 0; i < nXsum; i++) {
             ySum = xSum;
-            for (int j = 0; j < nYsum; j++) {
+            for (psS32 j = 0; j < nYsum; j++) {
                 zSum = ySum;
-                for (int k = 0; k < nZsum; k++) {
+                for (psS32 k = 0; k < nZsum; k++) {
                     sums[i][j][k] = zSum;
                     zSum *= z;
@@ -240,5 +206,5 @@
 /******************************************************************************
     BuildSums4D(sums, x, y, z, t, nXterm, nYterm, nZterm, nTterm). equiv to
-    BuildSums2D(). The result is returned as a double ****
+    BuildSums2D(). The result is returned as a psF64 ****
 *****************************************************************************/
 static psF64 ****BuildSums4D(
@@ -268,9 +234,9 @@
     if (sums == NULL) {
         sums = (psF64 ****) psAlloc (nXsum*sizeof(psF64));
-        for (int i = 0; i < nXsum; i++) {
+        for (psS32 i = 0; i < nXsum; i++) {
             sums[i] = (psF64 ***) psAlloc (nYsum*sizeof(psF64));
-            for (int j = 0; j < nYsum; j++) {
+            for (psS32 j = 0; j < nYsum; j++) {
                 sums[i][j] = (psF64 **) psAlloc (nZsum*sizeof(psF64));
-                for (int k = 0; k < nZsum; k++) {
+                for (psS32 k = 0; k < nZsum; k++) {
                     sums[i][j][k] = (psF64 *) psAlloc (nTsum*sizeof(psF64));
                 }
@@ -281,11 +247,11 @@
 
     tSum = 1.0;
-    for (int m = 0; m < nTsum; m++) {
+    for (psS32 m = 0; m < nTsum; m++) {
         zSum = tSum;
-        for (int k = 0; k < nZsum; k++) {
+        for (psS32 k = 0; k < nZsum; k++) {
             ySum = zSum;
-            for (int j = 0; j < nYsum; j++) {
+            for (psS32 j = 0; j < nYsum; j++) {
                 xSum = ySum;
-                for (int i = 0; i < nXsum; i++) {
+                for (psS32 i = 0; i < nXsum; i++) {
                     sums[i][j][k][m] = xSum;
                     xSum *= x;
@@ -299,4 +265,11 @@
     return (sums);
 }
+
+/******************************************************************************
+ ******************************************************************************
+ Analytical 1-D fitting routines.
+ ******************************************************************************
+ *****************************************************************************/
+
 
 /******************************************************************************
@@ -356,7 +329,7 @@
         NUM_DATA = x->n;
     }
-    // XX: psTrace    printf("vectorFitPolynomial1DChebySlow(): NUM_DATA is %d\n", NUM_DATA);
-
-    psPolynomial1D **chebPolys = createChebyshevPolys(NUM_POLY);
+    // psTrace    printf("vectorFitPolynomial1DChebySlow(): NUM_DATA is %d\n", NUM_DATA);
+
+    psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(NUM_POLY);
     if (0) {
         for (psS32 j = 0; j < NUM_POLY; j++) {
@@ -416,14 +389,16 @@
     }
 
-    // GaussJordan version
     if (0) {
-        // does the solution in place
-        // XXX: Check error codes!
-        psGaussJordan (A, B);
-
-        // the first nTerm entries in B correspond directly to the desired
-        // polynomial coefficients.  this is only true for the 1D case
-        for (psS32 k = 0; k < NUM_POLY; k++) {
-            myPoly->coeff[k] = B->data.F64[k];
+        // GaussJordan version
+        if (false == psGaussJordan(A, B)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
+            psFree(myPoly);
+            myPoly = NULL;
+        } else {
+            // the first nTerm entries in B correspond directly to the desired
+            // polynomial coefficients.  this is only true for the 1D case
+            for (psS32 k = 0; k < NUM_POLY; k++) {
+                myPoly->coeff[k] = B->data.F64[k];
+            }
         }
     } else {
@@ -433,4 +408,5 @@
         psVector* coeffs = NULL;
 
+        // XXX: Check return codes.
         ALUD = psImageAlloc(NUM_POLY, NUM_POLY, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
@@ -468,5 +444,5 @@
     the [-1:1] interval.
  
-XXX: mask, maskValue, yErr are currently ignored.
+    NOTE: mask, maskValue, yErr are ignored with this function.
 *****************************************************************************/
 static psPolynomial1D *vectorFitPolynomial1DChebyFast(
@@ -501,5 +477,5 @@
     tmpScalar.type.type = PS_TYPE_F64;
 
-    // XXX: These assignments appear too simple to warrant code and
+    // These assignments appear too simple to warrant code and
     // variable declarations.  I retain them here to maintain coherence
     // with the NR code.
@@ -538,10 +514,15 @@
             fScalar = p_psVectorInterpolate((psVector *) x, (psVector *) y,
                                             3, &tmpScalar);
+            if (fScalar == NULL) {
+                psError(PS_ERR_UNKNOWN, false, "Could not perform vector interpolation.  Returning NULL.\n");
+                psFree(myPoly)
+                return(NULL);
+            }
+
             f->data.F64[i] = fScalar->data.F64;
             psFree(fScalar);
         }
 
-        psTrace(".psLib.dataManip.vectorFitPolynomial1DCheby", 6,
-                "(x, X, y, f(X)) is (%f, %f, %f, %f)\n",
+        psTrace(__func__, 6, "(x, X, y, f(X)) is (%f, %f, %f, %f)\n",
                 x->data.F64[i], X, y->data.F64[i], f->data.F64[i]);
     }
@@ -551,5 +532,4 @@
 
     fac = 2.0/((psF32) n);
-    // XXX: is this loop bound correct?
     for (j=0;j<myPoly->nX+1;j++) {
         sum = 0.0;
@@ -573,5 +553,5 @@
 PS_TYPE_F64.
  *****************************************************************************/
-psPolynomial1D* VectorFitPolynomial1DOrd(
+static psPolynomial1D* VectorFitPolynomial1DOrd(
     psPolynomial1D* myPoly,
     const psVector *mask,
@@ -581,6 +561,5 @@
     const psVector *x)
 {
-    psTrace(__func__, 4, "---- VectorFitPolynomial1DOrd() begin ----\n");
-    // XXX: these ASSERTS are redundant.
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
@@ -600,6 +579,6 @@
     }
 
-    psImage*     A = NULL;
-    psVector*    B = NULL;
+    psImage* A = NULL;
+    psVector* B = NULL;
     psVector* xSums = NULL;
     psS32 nTerm;
@@ -607,5 +586,5 @@
     psF64 wt;
 
-    if (psTraceGetLevel (__func__) >= 6) {
+    if (psTraceGetLevel(__func__) >= 6) {
         psTrace(__func__, 6, "VectorFitPolynomial1D()\n");
         for (psS32 i = 0; i < f->n; i++) {
@@ -627,13 +606,15 @@
     nOrder = nTerm - 1;
 
-    A     = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
-    B     = psVectorAlloc(nTerm, PS_TYPE_F64);
-
-    //
+    A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
+    B = psVectorAlloc(nTerm, PS_TYPE_F64);
     // Initialize data structures.
-    // XXX: Use psLib function.
-    //
-    PS_VECTOR_SET_F64(B, 0.0);
-    PS_IMAGE_SET_F64(A, 0.0);
+    if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) {
+        psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B.  Returning NULL.\n");
+        psFree(myPoly);
+        psFree(A);
+        psFree(B);
+        psTrace(__func__, 4, "---- %s() End ----\n", __func__);
+        return(NULL);
+    }
 
     // xSums look like: 1, x, x^2, ... x^(2n+1)
@@ -641,5 +622,5 @@
     // XXX EAM : use temp pointers eg vB = B->data.F64 to save redirects
     // XXX EAM : this function is only valid for data vectors of F64
-    for (int k = 0; k < f->n; k++) {
+    for (psS32 k = 0; k < f->n; k++) {
         if ((mask != NULL) && (mask->data.U8[k] && maskValue)) {
             continue;
@@ -657,5 +638,5 @@
             wt = (fErr->data.F64[k] == 0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]);
         }
-        for (int i = 0; i < nTerm; i++) {
+        for (psS32 i = 0; i < nTerm; i++) {
             B->data.F64[i] += f->data.F64[k] * xSums->data.F64[i] * wt;
         }
@@ -663,6 +644,6 @@
         // we could skip half of the array and assign at the end
         // we must handle masked orders
-        for (int i = 0; i < nTerm; i++) {
-            for (int j = 0; j < nTerm; j++) {
+        for (psS32 i = 0; i < nTerm; i++) {
+            for (psS32 j = 0; j < nTerm; j++) {
                 A->data.F64[i][j] += xSums->data.F64[i + j] * wt;
             }
@@ -670,15 +651,18 @@
     }
 
-    // GaussJordan version
     if (0) {
-        // does the solution in place
-        // XXX: Check error codes!
-        psGaussJordan (A, B);
-
-        // the first nTerm entries in B correspond directly to the desired
-        // polynomial coefficients.  this is only true for the 1D case
-        for (int k = 0; k < nTerm; k++) {
-            myPoly->coeff[k] = B->data.F64[k];
-        }
+        // GaussJordan version
+        if (false == psGaussJordan(A, B)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
+            psFree(myPoly);
+            myPoly = NULL;
+        } else {
+            // the first nTerm entries in B correspond directly to the desired
+            // polynomial coefficients.  this is only true for the 1D case
+            for (psS32 k = 0; k < nTerm; k++) {
+                myPoly->coeff[k] = B->data.F64[k];
+            }
+        }
+
     } else {
         // LUD version of the fit
@@ -687,8 +671,9 @@
         psVector* coeffs = NULL;
 
+        // XXX: Check return codes.
         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
         coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-        for (int k = 0; k < nTerm; k++) {
+        for (psS32 k = 0; k < nTerm; k++) {
             myPoly->coeff[k] = coeffs->data.F64[k];
         }
@@ -703,5 +688,5 @@
     psFree(xSums);
 
-    psTrace(__func__, 4, "---- VectorFitPolynomial1DOrd() End ----\n");
+    psTrace(__func__, 4, "---- %s() End ----\n", __func__);
     return (myPoly);
 }
@@ -713,6 +698,4 @@
 Types F32 and F64 are supported, however, type F32 is done via vector
 conversion only.
- 
-XXX: Put the asserts at the beginning.  Why is this not failing tests?
  *****************************************************************************/
 psPolynomial1D *psVectorFitPolynomial1D(
@@ -748,5 +731,5 @@
 
     if (f->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
     } else {
         f64 = (psVector *) f;
@@ -755,5 +738,5 @@
     if (x != NULL) {
         if (x->type.type != PS_TYPE_F64) {
-            PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+            x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
         } else {
             x64 = (psVector *) x;
@@ -763,5 +746,5 @@
     if (fErr != NULL) {
         if (fErr->type.type != PS_TYPE_F64) {
-            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
         } else {
             fErr64 = (psVector *) fErr;
@@ -784,5 +767,5 @@
         if (x == NULL) {
             // If x==NULL, create an x64 vector with x values set to (-1:1).
-            PS_VECTOR_GEN_CHEBY_INDEX(x64, f64->n);
+            PS_VECTOR_GEN_CHEBY_INDEX(x64, f64->n, PS_TYPE_F64);
         }
 
@@ -850,10 +833,11 @@
         } else if (poly->type == PS_POLYNOMIAL_CHEB) {
             if (f->type.type == PS_TYPE_F32) {
-                PS_VECTOR_GEN_CHEBY_INDEX_F32(x, f->n);
+                PS_VECTOR_GEN_CHEBY_INDEX(x, f->n, PS_TYPE_F32);
             } else if (f->type.type == PS_TYPE_F64) {
-                PS_VECTOR_GEN_CHEBY_INDEX_F64(x, f->n);
+                PS_VECTOR_GEN_CHEBY_INDEX(x, f->n, PS_TYPE_F64);
             }
         } else {
-            printf("XXX: error, bad poly type.\n");
+            psError(PS_ERR_UNKNOWN, false, "Error, bad poly type.\n");
+            return(NULL);
         }
     }
@@ -885,7 +869,7 @@
 
     //
-    for (int N = 0; N < stats->clipIter; N++) {
+    for (psS32 N = 0; N < stats->clipIter; N++) {
         psTrace(__func__, 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n");
-        int Nkeep = 0;
+        psS32 Nkeep = 0;
         if (psTraceGetLevel(__func__) >= 6) {
             if (mask != NULL) {
@@ -897,10 +881,17 @@
         poly = psVectorFitPolynomial1D(poly, mask, maskValue, f, fErr, x);
         if (poly == NULL) {
-            printf("XXX: psVectorFitPolynomial1D() returned NULL: Generate error, free data.\n");
+            psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
+            if (xIn == NULL) {
+                psFree(x);
+            }
             return(NULL);
         }
 
-        // XXX: Check error codes
         fit = psPolynomial1DEvalVector(poly, x);
+        if (fit == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
+            psFree(resid)
+            return(NULL);
+        }
         for (psS32 i = 0 ; i < f->n ; i++) {
             if (f->type.type == PS_TYPE_F64) {
@@ -930,5 +921,5 @@
         // we are masking out any point which is out of range
         // recovery is not allowed with this scheme
-        for (int i = 0; i < resid->n; i++) {
+        for (psS32 i = 0; i < resid->n; i++) {
             if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
                 continue;
@@ -954,5 +945,5 @@
 
         //
-        // XXX: We should probably exit this loop if no new elements were masked
+        // We should probably exit this loop if no new elements were masked
         // since the polynomial fit won't change.
         //
@@ -984,7 +975,6 @@
 pairs.  All non-NULL vectors must be of type PS_TYPE_F64.
  
-// XXX: What about the maskValue?
  *****************************************************************************/
-psPolynomial2D* VectorFitPolynomial2DOrd(
+static psPolynomial2D* VectorFitPolynomial2DOrd(
     psPolynomial2D* myPoly,
     const psVector* mask,
@@ -995,6 +985,5 @@
     const psVector *y)
 {
-    psTrace(__func__, 4, "---- VectorFitPolynomial2DOrd() begin ----\n");
-    // These ASSERTS are redundant.
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
@@ -1025,5 +1014,4 @@
     psS32 nTerm;
 
-    // XXX:Watch for changes to the psPolys: nTerm != nOrder.
     psS32 nXterm = 1 + myPoly->nX;
     psS32 nYterm = 1 + myPoly->nY;
@@ -1032,16 +1020,18 @@
     A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
     B = psVectorAlloc(nTerm, PS_TYPE_F64);
-
-    //
     // Initialize data structures.
-    // XXX: Use psLib function.
-    //
-    PS_VECTOR_SET_F64(B, 0.0);
-    PS_IMAGE_SET_F64(A, 0.0);
+    if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) {
+        psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B.  Returning NULL.\n");
+        psFree(myPoly);
+        psFree(A);
+        psFree(B);
+        psTrace(__func__, 4, "---- %s() End ----\n", __func__);
+        return(NULL);
+    }
 
     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)
 
     // Build the B and A data structs.
-    for (int k  = 0; k < x->n; k++) {
+    for (psS32 k  = 0; k < x->n; k++) {
         if ((mask != NULL) && (mask->data.U8[k] & maskValue)) {
             continue;
@@ -1059,14 +1049,14 @@
         // we could skip half of the array and assign at the end
         // we must handle masked orders
-        for (int n = 0; n < nXterm; n++) {
-            for (int m = 0; m < nYterm; m++) {
+        for (psS32 n = 0; n < nXterm; n++) {
+            for (psS32 m = 0; m < nYterm; m++) {
                 B->data.F64[n+m*nXterm] += f->data.F64[k] * Sums->data.F64[n][m] * wt;
             }
         }
 
-        for (int i = 0; i < nXterm; i++) {
-            for (int j = 0; j < nYterm; j++) {
-                for (int n = 0; n < nXterm; n++) {
-                    for (int m = 0; m < nYterm; m++) {
+        for (psS32 i = 0; i < nXterm; i++) {
+            for (psS32 j = 0; j < nYterm; j++) {
+                for (psS32 n = 0; n < nXterm; n++) {
+                    for (psS32 m = 0; m < nYterm; m++) {
                         A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
                     }
@@ -1076,12 +1066,14 @@
     }
 
-    // does the solution in place
-    // XXX: Check error codes!
-    psGaussJordan (A, B);
-
-    // select the appropriate solution entries
-    for (int n = 0; n < nXterm; n++) {
-        for (int m = 0; m < nYterm; m++) {
-            myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
+    if (false == psGaussJordan(A, B)) {
+        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
+        psFree(myPoly);
+        myPoly = NULL;
+    } else {
+        // select the appropriate solution entries
+        for (psS32 n = 0; n < nXterm; n++) {
+            for (psS32 m = 0; m < nYterm; m++) {
+                myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
+            }
         }
     }
@@ -1091,5 +1083,5 @@
     psFree(Sums);
 
-    psTrace(__func__, 4, "---- VectorFitPolynomial2DOrd() end ----\n");
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return (myPoly);
 }
@@ -1139,5 +1131,5 @@
     //
     if (f->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
     } else {
         f64 = (psVector *) f;
@@ -1145,5 +1137,5 @@
 
     if (x->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+        x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
     } else {
         x64 = (psVector *) x;
@@ -1151,5 +1143,5 @@
 
     if (y->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(y64, y);
+        y64 = psVectorCopy(NULL, y, PS_TYPE_F64);
     } else {
         y64 = (psVector *) y;
@@ -1161,5 +1153,5 @@
     if (fErr != NULL) {
         if (fErr->type.type != PS_TYPE_F64) {
-            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
         } else {
             fErr64 = (psVector *) fErr;
@@ -1272,6 +1264,6 @@
 
     // clipping range defined by min and max and/or clipSigma
-    float minClipSigma;
-    float maxClipSigma;
+    psF32 minClipSigma;
+    psF32 maxClipSigma;
     if (isfinite(stats->max)) {
         maxClipSigma = fabs(stats->max);
@@ -1294,7 +1286,7 @@
     psTrace(__func__, 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
 
-    for (int N = 0; N < stats->clipIter; N++) {
+    for (psS32 N = 0; N < stats->clipIter; N++) {
         psTrace(__func__, 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n");
-        int Nkeep = 0;
+        psS32 Nkeep = 0;
         if (psTraceGetLevel(__func__) >= 6) {
             if (mask != NULL) {
@@ -1306,6 +1298,16 @@
 
         poly = psVectorFitPolynomial2D(poly, mask, maskValue, f, fErr, x, y);
-        // XXX: Check error codes
+        if (poly == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data.  Returning NULL.\n");
+            psFree(resid)
+            return(NULL);
+        }
+
         psVector *fit = psPolynomial2DEvalVector(poly, x, y);
+        if (fit == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
+            psFree(resid)
+            return(NULL);
+        }
 
         for (psS32 i = 0 ; i < f->n ; i++) {
@@ -1328,15 +1330,20 @@
         }
 
-        // XXX: Check error codes
         stats  = psVectorStats(stats, resid, NULL, mask, maskValue);
+        if (stats == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
+            psFree(resid)
+            psFree(fit)
+            return(NULL);
+        }
         psTrace(__func__, 6, "Median is %f\n", stats->sampleMedian);
         psTrace(__func__, 6, "Stdev is %f\n", stats->sampleStdev);
-        float minClipValue = -minClipSigma*stats->sampleStdev;
-        float maxClipValue = +maxClipSigma*stats->sampleStdev;
+        psF32 minClipValue = -minClipSigma*stats->sampleStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
 
         // set mask if pts are not valid
         // we are masking out any point which is out of range
         // recovery is not allowed with this scheme
-        for (int i = 0; i < resid->n; i++) {
+        for (psS32 i = 0; i < resid->n; i++) {
             if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
                 continue;
@@ -1388,9 +1395,6 @@
 y, z)-(f) pairs.  All non-NULL vectors must be of type PS_TYPE_F64.
  
-XXX: This routine needs to be written.  Currently, this is simply a shell.  We
-can assume that all vectors have been converted to F64, that (f, x, y, z) are
-non-null and F64.  fErr may be NULL, but will be F64 is not.
  *****************************************************************************/
-psPolynomial3D* VectorFitPolynomial3DOrd(
+static psPolynomial3D* VectorFitPolynomial3DOrd(
     psPolynomial3D* myPoly,
     const psVector* mask,
@@ -1402,6 +1406,5 @@
     const psVector *z)
 {
-    psTrace(__func__, 4, "---- VectorFitPolynomial3DOrd() begin ----\n");
-    // These ASSERTS are redundant.
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
@@ -1435,5 +1438,4 @@
     psS32 nTerm;
 
-    // XXX:Watch for changes to the psPolys: nTerm != nOrder.
     psS32 nXterm = 1 + myPoly->nX;
     psS32 nYterm = 1 + myPoly->nY;
@@ -1443,13 +1445,18 @@
     A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
     B = psVectorAlloc(nTerm, PS_TYPE_F64);
-
     // Initialize data structures.
-    psVectorInit (B, 0.0);
-    psImageInit (A, 0.0);
+    if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) {
+        psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B.  Returning NULL.\n");
+        psFree(myPoly);
+        psFree(A);
+        psFree(B);
+        psTrace(__func__, 4, "---- %s() End ----\n", __func__);
+        return(NULL);
+    }
 
     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ...
 
     // Build the B and A data structs.
-    for (int k  = 0; k < x->n; k++) {
+    for (psS32 k  = 0; k < x->n; k++) {
         if ((mask != NULL) && (mask->data.U8[k] & maskValue)) {
             continue;
@@ -1467,11 +1474,11 @@
         // we could skip half of the array and assign at the end
         // we must handle masked orders
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
                     if (myPoly->mask[ix][iy][iz]) {
                         continue;
                     } else {
-                        int nx = ix + iy*nXterm + iz*nXterm*nYterm;
+                        psS32 nx = ix + iy*nXterm + iz*nXterm*nYterm;
                         B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz] * wt;
                     }
@@ -1480,16 +1487,16 @@
         }
 
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
                     if (myPoly->mask[ix][iy][iz])
                         continue;
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
-                    for (int jx = 0; jx < nXterm; jx++) {
-                        for (int jy = 0; jy < nYterm; jy++) {
-                            for (int jz = 0; jz < nZterm; jz++) {
+                    psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
+                    for (psS32 jx = 0; jx < nXterm; jx++) {
+                        for (psS32 jy = 0; jy < nYterm; jy++) {
+                            for (psS32 jz = 0; jz < nZterm; jz++) {
                                 if (myPoly->mask[jx][jy][jz])
                                     continue;
-                                int ny = jx+jy*nXterm+jz*nXterm*nYterm;
+                                psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm;
                                 A->data.F64[nx][ny] += Sums[ix+jx][iy+jy][iz+jz] * wt;
                             }
@@ -1501,15 +1508,15 @@
     }
 
-    for (int ix = 0; ix < nXterm; ix++) {
-        for (int iy = 0; iy < nYterm; iy++) {
-            for (int iz = 0; iz < nZterm; iz++) {
+    for (psS32 ix = 0; ix < nXterm; ix++) {
+        for (psS32 iy = 0; iy < nYterm; iy++) {
+            for (psS32 iz = 0; iz < nZterm; iz++) {
                 if (!myPoly->mask[ix][iy][iz])
                     continue;
-                int nx = ix+iy*nXterm+iz*nXterm*nYterm;
+                psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
                 B->data.F64[nx] = 0;
-                for (int jx = 0; jx < nXterm; jx++) {
-                    for (int jy = 0; jy < nYterm; jy++) {
-                        for (int jz = 0; jz < nZterm; jz++) {
-                            int ny = jx+jy*nXterm+jz*nXterm*nYterm;
+                for (psS32 jx = 0; jx < nXterm; jx++) {
+                    for (psS32 jy = 0; jy < nYterm; jy++) {
+                        for (psS32 jz = 0; jz < nZterm; jz++) {
+                            psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm;
                             A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
                         }
@@ -1523,10 +1530,10 @@
         // does the solution in place
         // The matrices were overflowing, so I switched to LUD.
-        if (false == psGaussJordan (A, B)) {
+        if (false == psGaussJordan(A, B)) {
             psFree(A);
             psFree(B);
 
-            for (int ix = 0; ix < 2*nXterm; ix++) {
-                for (int iy = 0; iy < 2*nYterm; iy++) {
+            for (psS32 ix = 0; ix < 2*nXterm; ix++) {
+                for (psS32 iy = 0; iy < 2*nYterm; iy++) {
                     psFree(Sums[ix][iy]);
                 }
@@ -1538,8 +1545,8 @@
         }
         // select the appropriate solution entries
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
                     myPoly->coeff[ix][iy][iz] = B->data.F64[nx];
                     myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
@@ -1553,4 +1560,5 @@
         psVector* coeffs = NULL;
 
+        // XXX: Check return codes.
         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
@@ -1558,8 +1566,8 @@
 
         // select the appropriate solution entries
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
                     myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx];
                     myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
@@ -1575,6 +1583,6 @@
     psFree(B);
 
-    for (int ix = 0; ix < 2*nXterm; ix++) {
-        for (int iy = 0; iy < 2*nYterm; iy++) {
+    for (psS32 ix = 0; ix < 2*nXterm; ix++) {
+        for (psS32 iy = 0; iy < 2*nYterm; iy++) {
             psFree(Sums[ix][iy]);
         }
@@ -1583,5 +1591,5 @@
     psFree(Sums);
 
-    psTrace(__func__, 4, "---- VectorFitPolynomial3DOrd() end ----\n");
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return (myPoly);
 }
@@ -1621,11 +1629,8 @@
     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    //    PS_ASSERT_VECTOR_TYPE(x, f->type.type, NULL);
     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    //    PS_ASSERT_VECTOR_TYPE(y, f->type.type, NULL);
     PS_ASSERT_VECTOR_NON_NULL(z, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    //    PS_ASSERT_VECTOR_TYPE(z, f->type.type, NULL);
     if (fErr != NULL) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
@@ -1637,15 +1642,15 @@
     //
     if (f->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
     } else {
         f64 = (psVector *) f;
     }
     if (x->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+        x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
     } else {
         x64 = (psVector *) x;
     }
     if (y->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(y64, y);
+        y64 = psVectorCopy(NULL, y, PS_TYPE_F64);
     } else {
         y64 = (psVector *) y;
@@ -1653,5 +1658,5 @@
 
     if (z->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(z64, z);
+        z64 = psVectorCopy(NULL, z, PS_TYPE_F64);
     } else {
         z64 = (psVector *) z;
@@ -1660,5 +1665,5 @@
     if (fErr != NULL) {
         if (fErr->type.type != PS_TYPE_F64) {
-            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
         } else {
             fErr64 = (psVector *) fErr;
@@ -1788,6 +1793,6 @@
 
     // clipping range defined by min and max and/or clipSigma
-    float minClipSigma;
-    float maxClipSigma;
+    psF32 minClipSigma;
+    psF32 maxClipSigma;
     if (isfinite(stats->max)) {
         maxClipSigma = fabs(stats->max);
@@ -1811,7 +1816,7 @@
     psTrace(__func__, 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
 
-    for (int N = 0; N < stats->clipIter; N++) {
+    for (psS32 N = 0; N < stats->clipIter; N++) {
         psTrace(__func__, 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n");
-        int Nkeep = 0;
+        psS32 Nkeep = 0;
         if (psTraceGetLevel(__func__) >= 6) {
             if (mask != NULL) {
@@ -1823,6 +1828,16 @@
 
         poly = psVectorFitPolynomial3D(poly, mask, maskValue, f, fErr, x, y, z);
-        // XXX: Check error codes
+        if (poly == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data.  Returning NULL.\n");
+            psFree(resid)
+            psFree(fit)
+            return(NULL);
+        }
         fit = psPolynomial3DEvalVector(poly, x, y, z);
+        if (fit == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
+            psFree(resid)
+            return(NULL);
+        }
         for (psS32 i = 0 ; i < f->n ; i++) {
             if (f->type.type == PS_TYPE_F64) {
@@ -1844,15 +1859,21 @@
         }
 
-        // XXX: Check error codes
         stats  = psVectorStats(stats, resid, NULL, mask, maskValue);
+        if (stats == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
+            psFree(resid)
+            psFree(fit)
+            return(NULL);
+        }
+
         psTrace(__func__, 6, "Median is %f\n", stats->sampleMedian);
         psTrace(__func__, 6, "Stdev is %f\n", stats->sampleStdev);
-        float minClipValue = -minClipSigma*stats->sampleStdev;
-        float maxClipValue = +maxClipSigma*stats->sampleStdev;
+        psF32 minClipValue = -minClipSigma*stats->sampleStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
 
         // set mask if pts are not valid
         // we are masking out any point which is out of range
         // recovery is not allowed with this scheme
-        for (int i = 0; i < resid->n; i++) {
+        for (psS32 i = 0; i < resid->n; i++) {
             if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
                 continue;
@@ -1902,9 +1923,6 @@
 y, z, t)-(f) pairs.  All non-NULL vectors must be of type PS_TYPE_F64.
  
-XXX: This routine needs to be written.  Currently, this is simply a shell.  We
-can assume that all vectors have been converted to F64, that (f, x, y, z) are
-non-null and F64.  fErr may be NULL, but will be F64 is not.
  *****************************************************************************/
-psPolynomial4D* VectorFitPolynomial4DOrd(
+static psPolynomial4D* VectorFitPolynomial4DOrd(
     psPolynomial4D* myPoly,
     const psVector* mask,
@@ -1918,5 +1936,4 @@
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    // These ASSERTS are redundant.
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
@@ -1954,5 +1971,4 @@
     psS32 nTerm;
 
-    // XXX:Watch for changes to the psPolys: nTerm != nOrder.
     psS32 nXterm = 1 + myPoly->nX;
     psS32 nYterm = 1 + myPoly->nY;
@@ -1963,13 +1979,18 @@
     A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
     B = psVectorAlloc(nTerm, PS_TYPE_F64);
-
     // Initialize data structures.
-    psVectorInit (B, 0.0);
-    psImageInit (A, 0.0);
+    if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) {
+        psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B.  Returning NULL.\n");
+        psFree(myPoly);
+        psFree(A);
+        psFree(B);
+        psTrace(__func__, 4, "---- %s() End ----\n", __func__);
+        return(NULL);
+    }
 
     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ...
 
     // Build the B and A data structs.
-    for (int k  = 0; k < x->n; k++) {
+    for (psS32 k  = 0; k < x->n; k++) {
         if ((mask != NULL) && (mask->data.U8[k] & maskValue)) {
             continue;
@@ -1987,11 +2008,11 @@
         // we could skip half of the array and assign at the end
         // we must handle masked orders
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    for (int it = 0; it < nTterm; it++) {
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    for (psS32 it = 0; it < nTterm; it++) {
                         if (myPoly->mask[ix][iy][iz][it])
                             continue;
-                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+                        psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
                         B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz][it] * wt;
                     }
@@ -2000,18 +2021,18 @@
         }
 
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    for (int it = 0; it < nTterm; it++) {
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    for (psS32 it = 0; it < nTterm; it++) {
                         if (myPoly->mask[ix][iy][iz][it])
                             continue;
-                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
-                        for (int jx = 0; jx < nXterm; jx++) {
-                            for (int jy = 0; jy < nYterm; jy++) {
-                                for (int jz = 0; jz < nZterm; jz++) {
-                                    for (int jt = 0; jt < nTterm; jt++) {
+                        psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+                        for (psS32 jx = 0; jx < nXterm; jx++) {
+                            for (psS32 jy = 0; jy < nYterm; jy++) {
+                                for (psS32 jz = 0; jz < nZterm; jz++) {
+                                    for (psS32 jt = 0; jt < nTterm; jt++) {
                                         if (myPoly->mask[jx][jy][jz][jt])
                                             continue;
-                                        int ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
+                                        psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
                                         A->data.F64[nx][ny]+= Sums[ix+jx][iy+jy][iz+jz][it+jt] * wt;
                                     }
@@ -2025,17 +2046,17 @@
     }
 
-    for (int ix = 0; ix < nXterm; ix++) {
-        for (int iy = 0; iy < nYterm; iy++) {
-            for (int iz = 0; iz < nZterm; iz++) {
-                for (int it = 0; it < nTterm; it++) {
+    for (psS32 ix = 0; ix < nXterm; ix++) {
+        for (psS32 iy = 0; iy < nYterm; iy++) {
+            for (psS32 iz = 0; iz < nZterm; iz++) {
+                for (psS32 it = 0; it < nTterm; it++) {
                     if (!myPoly->mask[ix][iy][iz][it])
                         continue;
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+                    psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
                     B->data.F64[nx] = 0;
-                    for (int jx = 0; jx < nXterm; jx++) {
-                        for (int jy = 0; jy < nYterm; jy++) {
-                            for (int jz = 0; jz < nZterm; jz++) {
-                                for (int jt = 0; jt < nTterm; jt++) {
-                                    int ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
+                    for (psS32 jx = 0; jx < nXterm; jx++) {
+                        for (psS32 jy = 0; jy < nYterm; jy++) {
+                            for (psS32 jz = 0; jz < nZterm; jz++) {
+                                for (psS32 jt = 0; jt < nTterm; jt++) {
+                                    psS32 ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
                                     A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
                                 }
@@ -2051,11 +2072,11 @@
     if (0) {
         // does the solution in place
-        // XXX: The GaussJordan version was overflowing, so I'm using LUD.
+        // The GaussJordan version was overflowing, so I'm using LUD.
         if (false == psGaussJordan(A, B)) {
             psFree(A);
             psFree(B);
-            for (int ix = 0; ix < 2*nXterm; ix++) {
-                for (int iy = 0; iy < 2*nYterm; iy++) {
-                    for (int iz = 0; iz < 2*nZterm; iz++) {
+            for (psS32 ix = 0; ix < 2*nXterm; ix++) {
+                for (psS32 iy = 0; iy < 2*nYterm; iy++) {
+                    for (psS32 iz = 0; iz < 2*nZterm; iz++) {
                         psFree(Sums[ix][iy][iz]);
                     }
@@ -2070,9 +2091,9 @@
 
         // select the appropriate solution entries
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    for (int it = 0; it < nTterm; it++) {
-                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    for (psS32 it = 0; it < nTterm; it++) {
+                        psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
                         myPoly->coeff[ix][iy][iz][it] = B->data.F64[nx];
                         myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
@@ -2087,4 +2108,5 @@
         psVector* coeffs = NULL;
 
+        // XXX: Check return codes.
         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
@@ -2092,9 +2114,9 @@
 
         // select the appropriate solution entries
-        for (int ix = 0; ix < nXterm; ix++) {
-            for (int iy = 0; iy < nYterm; iy++) {
-                for (int iz = 0; iz < nZterm; iz++) {
-                    for (int it = 0; it < nTterm; it++) {
-                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+        for (psS32 ix = 0; ix < nXterm; ix++) {
+            for (psS32 iy = 0; iy < nYterm; iy++) {
+                for (psS32 iz = 0; iz < nZterm; iz++) {
+                    for (psS32 it = 0; it < nTterm; it++) {
+                        psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
                         myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx];
                         myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
@@ -2113,7 +2135,7 @@
     psFree(B);
 
-    for (int ix = 0; ix < 2*nXterm; ix++) {
-        for (int iy = 0; iy < 2*nYterm; iy++) {
-            for (int iz = 0; iz < 2*nZterm; iz++) {
+    for (psS32 ix = 0; ix < 2*nXterm; ix++) {
+        for (psS32 iy = 0; iy < 2*nYterm; iy++) {
+            for (psS32 iz = 0; iz < 2*nZterm; iz++) {
                 psFree(Sums[ix][iy][iz]);
             }
@@ -2180,25 +2202,25 @@
     //
     if (f->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+        f64 = psVectorCopy(NULL, f, PS_TYPE_F64);
     } else {
         f64 = (psVector *) f;
     }
     if (x->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+        x64 = psVectorCopy(NULL, x, PS_TYPE_F64);
     } else {
         x64 = (psVector *) x;
     }
     if (y->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(y64, y);
+        y64 = psVectorCopy(NULL, y, PS_TYPE_F64);
     } else {
         y64 = (psVector *) y;
     }
     if (z->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(z64, z);
+        z64 = psVectorCopy(NULL, z, PS_TYPE_F64);
     } else {
         z64 = (psVector *) z;
     }
     if (t->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(t64, t);
+        t64 = psVectorCopy(NULL, t, PS_TYPE_F64);
     } else {
         t64 = (psVector *) t;
@@ -2209,5 +2231,5 @@
     if (fErr != NULL) {
         if (fErr->type.type != PS_TYPE_F64) {
-            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+            fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64);
         } else {
             fErr64 = (psVector *) fErr;
@@ -2310,5 +2332,4 @@
 
 
-// XXX: Add F64 support.
 psPolynomial4D *psVectorClipFitPolynomial4D(
     psPolynomial4D *poly,
@@ -2356,6 +2377,6 @@
 
     // clipping range defined by min and max and/or clipSigma
-    float minClipSigma;
-    float maxClipSigma;
+    psF32 minClipSigma;
+    psF32 maxClipSigma;
     if (isfinite(stats->max)) {
         maxClipSigma = fabs(stats->max);
@@ -2379,7 +2400,7 @@
     psTrace(__func__, 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
 
-    for (int N = 0; N < stats->clipIter; N++) {
+    for (psS32 N = 0; N < stats->clipIter; N++) {
         psTrace(__func__, 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n");
-        int Nkeep = 0;
+        psS32 Nkeep = 0;
         if (psTraceGetLevel(__func__) >= 6) {
             if (mask != NULL) {
@@ -2391,7 +2412,17 @@
 
         poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t);
-        // XXX: Check error codes
+        if (poly == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data.  Returning NULL.\n");
+            psFree(resid)
+            psFree(fit)
+            return(NULL);
+        }
+
         fit = psPolynomial4DEvalVector (poly, x, y, z, t);
-        // XXX: This is coded differently for 1D, 2D.  Consolidate.
+        if (fit == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector().  Returning NULL.\n");
+            psFree(resid)
+            return(NULL);
+        }
         for (psS32 i = 0 ; i < f->n ; i++) {
             if (f->type.type == PS_TYPE_F64) {
@@ -2401,5 +2432,4 @@
             }
         }
-        // XXX: Check error codes
 
         if (psTraceGetLevel(__func__) >= 6) {
@@ -2414,15 +2444,20 @@
         }
 
-        // XXX: Check error codes
         stats  = psVectorStats (stats, resid, NULL, mask, maskValue);
+        if (stats == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
+            psFree(resid)
+            psFree(fit)
+            return(NULL);
+        }
         psTrace(__func__, 6, "Median is %f\n", stats->sampleMedian);
         psTrace(__func__, 6, "Stdev is %f\n", stats->sampleStdev);
-        float minClipValue = -minClipSigma*stats->sampleStdev;
-        float maxClipValue = +maxClipSigma*stats->sampleStdev;
+        psF32 minClipValue = -minClipSigma*stats->sampleStdev;
+        psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
 
         // set mask if pts are not valid
         // we are masking out any point which is out of range
         // recovery is not allowed with this scheme
-        for (int i = 0; i < resid->n; i++) {
+        for (psS32 i = 0; i < resid->n; i++) {
             if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
                 continue;