Changeset 4991

Timestamp:

Sep 11, 2005, 12:18:40 PM (21 years ago)

Author:

gusciora

Message:

....

Location:

trunk/psLib/src

Files:

• astro/psCoord.c (modified) (5 diffs)
• math/psConstants.h (modified) (3 diffs)
• math/psMinimize.c (modified) (40 diffs)
• math/psMinimize.h (modified) (11 diffs)
• math/psPolynomial.c (modified) (2 diffs)
• math/psSpline.c (modified) (3 diffs)
• math/psSpline.h (modified) (2 diffs)
• math/psStats.c (modified) (7 diffs)

trunk/psLib/src/astro/psCoord.c

@@ -10 +10 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.85 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-07 21:33:48 $
+*  @version $Revision: 1.86 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-11 22:18:40 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -545 +545 @@
         int nSamples)
 {
-
-    // XXX: This does not yet use region and nSamples:  need to modify -rdd
-
     PS_ASSERT_PTR_NON_NULL(trans1, NULL);
     PS_ASSERT_PTR_NON_NULL(trans2, NULL);
@@ -905 +902 @@
 )
 {
-    /*
-        PS_ASSERT_PTR_NON_NULL(transformation, NULL);
-        PS_ASSERT_POLY_NON_NULL(transformation->x, NULL);
-        PS_ASSERT_POLY_NON_NULL(transformation->y, NULL);
-        PS_ASSERT_PTR_NON_NULL(coord, NULL);
-     
-        if (out == NULL) {
-            out = psPlaneAlloc();
-        }
-     
-        out->x = 0.0;
-        out->y = 0.0;
-        out->xErr = 0.0;
-        out->yErr = 0.0;
-     
-        psPolynomial2D *xPoly = transformation->x;
-        psPolynomial2D *yPoly = transformation->y;
-     
-        //
-        // Calculate the derivative with respect to x.
-        //
-        psF32 xSum = 1.0;
-        psF32 ySum = 1.0;
-     
-        // This loop starts at loop_x=1 since the derivative of the loop_x=0 terms are all 0.0
-        for (psS32 loop_x = 1; loop_x < xPoly->nX; loop_x++) {
-            ySum = 1.0;
-            for (psS32 loop_y = 0; loop_y < xPoly->nY; loop_y++) {
-                //
-                // For each iteration of the loop, we multiple the (x, y) coefficient
-                // by (coord->x^(loop_x-1) * coord->y^loop_y)
-                //
-     
-                out->x+= xPoly->coeff[loop_x][loop_y] * xSum * ySum;
-                ySum*= coord->y;
-            }
-            xSum*= coord->x;
-        }
-     
-        //
-        // Calculate the derivative with respect to x.
-        //
-        xSum = 1.0;
-     
-        // This loop starts at loop_y=1 since the derivative of the loop_y=0 terms are all 0.0
-        for (psS32 loop_x = 0; loop_x < yPoly->nX; loop_x++) {
-            ySum = 1.0;
-            for (psS32 loop_y = 1; loop_y < yPoly->nY; loop_y++) {
-                //
-                // For each iteration of the loop, we multiple the (x, y) coefficient
-                // by (coord->x^(loop_x-1) * coord->y^loop_y)
-                //
-     
-                out->y+= yPoly->coeff[loop_x][loop_y] * xSum * ySum;
-                ySum*= coord->y;
-            }
-            xSum*= coord->x;
-        }
-    */
+    PS_ASSERT_PTR_NON_NULL(transformation, NULL);
+    PS_ASSERT_POLY_NON_NULL(transformation->x, NULL);
+    PS_ASSERT_POLY_NON_NULL(transformation->y, NULL);
+    PS_ASSERT_PTR_NON_NULL(coord, NULL);
+
+    if (out == NULL) {
+        out = psPlaneAlloc();
+    }
+
+    out->x = 0.0;
+    out->y = 0.0;
+    out->xErr = 0.0;
+    out->yErr = 0.0;
+
+    psPolynomial2D *xPoly = transformation->x;
+    psPolynomial2D *yPoly = transformation->y;
+
+    //
+    // Calculate the derivative with respect to x.
+    //
+    psF32 xSum = 1.0;
+    psF32 ySum = 1.0;
+
+    // This loop starts at loop_x=1 since the derivative of the loop_x=0 terms are all 0.0
+    for (psS32 loop_x = 1; loop_x < xPoly->nX; loop_x++) {
+        ySum = 1.0;
+        for (psS32 loop_y = 0; loop_y < xPoly->nY; loop_y++) {
+            //
+            // For each iteration of the loop, we multiple the (x, y) coefficient
+            // by (coord->x^(loop_x-1) * coord->y^loop_y)
+            //
+
+            out->x+= xPoly->coeff[loop_x][loop_y] * xSum * ySum;
+            ySum*= coord->y;
+        }
+        xSum*= coord->x;
+    }
+
+    //
+    // Calculate the derivative with respect to x.
+    //
+    xSum = 1.0;
+
+    // This loop starts at loop_y=1 since the derivative of the loop_y=0 terms are all 0.0
+    for (psS32 loop_x = 0; loop_x < yPoly->nX; loop_x++) {
+        ySum = 1.0;
+        for (psS32 loop_y = 1; loop_y < yPoly->nY; loop_y++) {
+            //
+            // For each iteration of the loop, we multiple the (x, y) coefficient
+            // by (coord->x^(loop_x-1) * coord->y^loop_y)
+            //
+
+            out->y+= yPoly->coeff[loop_x][loop_y] * xSum * ySum;
+            ySum*= coord->y;
+        }
+        xSum*= coord->x;
+    }
     return(out);
 }
@@ -970 +965 @@
 {
     if(sphere == NULL) {
-        //        psError();
+        psError( PS_ERR_UNKNOWN, true, "psSphere argument is NULL.  Returning NULL.\n");
         return NULL;
     }
@@ -990 +985 @@
 {
     if(cube == NULL) {
-        //        psError();
+        psError( PS_ERR_UNKNOWN, true, "psCube argument is NULL.  Returning NULL.\n");
         return NULL;
     }

trunk/psLib/src/math/psConstants.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.75 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-07-20 01:21:13 $
+ *  @version $Revision: 1.76 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-11 22:18:40 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -357 +357 @@
 
 #define PS_VECTOR_PRINT_F32(NAME) \
-for (int my_i=0;my_i<(NAME)->n;my_i++) { \
-    printf("%s->data.F32[%d] is %f\n", #NAME, my_i, (NAME)->data.F32[my_i]); \
-} \
-printf("\n"); \
+if (NAME != NULL) { \
+    for (int my_i=0;my_i<(NAME)->n;my_i++) { \
+        printf("%s->data.F32[%d] is %f\n", #NAME, my_i, (NAME)->data.F32[my_i]); \
+    } \
+    printf("\n"); \
+} else {\
+    printf("MACRO WARNING: vector %s is NULL.\n", #NAME); \
+}\
 
 #define PS_VECTOR_PRINT_F64(NAME) \
-for (int my_i=0;my_i<(NAME)->n;my_i++) { \
-    printf("%s->data.F64[%d] is %f\n", #NAME, my_i, (NAME)->data.F64[my_i]); \
-} \
-printf("\n"); \
-
+if (NAME != NULL) { \
+    for (int my_i=0;my_i<(NAME)->n;my_i++) { \
+        printf("%s->data.F64[%d] is %f\n", #NAME, my_i, (NAME)->data.F64[my_i]); \
+    } \
+    printf("\n"); \
+} else {\
+    printf("MACRO WARNING: vector %s is NULL.\n", #NAME); \
+}\
 
 #define PS_VECTOR_CONVERT_F64_TO_F32_STATIC(OLD, NEW_PTR32, NEW_STATIC32) \
@@ -742 +749 @@
 #define PS_SQR(A) \
 ((A) * (A))
+
+# define PS_SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;}

trunk/psLib/src/math/psMinimize.c

@@ -8 +8 @@
  *
  *  @author GLG, MHPCC
+ *  @author EAM, IfA
  *
- *  @version $Revision: 1.134 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-07 21:35:50 $
+ *  @version $Revision: 1.135 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-11 22:18:40 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -32 +33 @@
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
+/* What are the following macros for?                                        */
 /*****************************************************************************/
 #define PS_SEG psLib
@@ -57 +59 @@
 /*****************************************************************************/
 
+/******************************************************************************
+ ******************************************************************************
+ Levenberg-Marquadt routines.
+ ******************************************************************************
+ *****************************************************************************/
 // measure the distance to the minimum assuming a linear model
+// XXX: Move this to the pub area.
 bool psMinimizeGaussNewtonDelta (psVector *delta,
                                  const psVector *params,
@@ -112 +120 @@
 }
 
-/******************************************************************************
-buildSums1D(x, polyOrder, sums): this routine calculates the powers of input
-parameter "x" between 0 and input parameter polyOrder.  The result is returned
-as a psVector sums.
- 
-XXX: Use a static vector.
- 
-XXX: should the argument be polyOrder or numTerms?
- *****************************************************************************/
-static void buildSums1D(psF64 x,
-                        psS32 polyOrder,
-                        psVector* sums)
-{
-    psS32 i = 0;
-    psF64 xSum = 0.0;
-    psS32 numTerms = polyOrder + 1;
-
-    if (sums == NULL) {
-        sums = psVectorAlloc(numTerms, PS_TYPE_F64);
-    } else if (numTerms > sums->n) {
-        sums = psVectorRealloc(sums, numTerms);
-    }
-
-    xSum = 1.0;
-    for (i = 0; i < numTerms; i++) {
-        sums->data.F64[i] = xSum;
-        xSum *= x;
-    }
-}
-
-/*****************************************************************************
-CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a
-tabulated function at n points, this routine calculates the second
-derivatives of the interpolating cubic splines at those n points.
- 
-The first and second derivatives at the endpoints, undefined in the SDR, are
-here defined to be 0.0.  They can be modified via ypo and yp1.
- 
-This routine assumes that vectors x and y are of the appropriate types/sizes
-(F32).
- 
-XXX: This algorithm is derived from the Numerical Recipes.
-XXX: use recycled vectors for internal data.
-XXX: do an F64 version?
- *****************************************************************************/
-static psF32 *calculateSecondDerivs(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
-                                    const psVector* y)        ///< Coordinates
-{
-    psTrace(".psLib.dataManip.calculateSecondDerivs", 4,
-            "---- calculateSecondDerivs() begin ----\n");
-
-    psS32 i;
-    psS32 k;
-    psF32 sig;
-    psF32 p;
-    psS32 n = y->n;
-    psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32));
-    psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32));
-    psF32 *X = (psF32 *) & (x->data.F32[0]);
-    psF32 *Y = (psF32 *) & (y->data.F32[0]);
-    psF32 qn;
-
-    // XXX: The second derivatives at the endpoints, undefined in the SDR,
-    // are set in psConstants.h: PS_LEFT_SPLINE_DERIV, PS_RIGHT_SPLINE_DERIV.
-    derivs2[0] = -0.5;
-    u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - PS_LEFT_SPLINE_DERIV);
-
-    for (i=1;i<=(n-2);i++) {
-        sig = (X[i] - X[i-1]) / (X[i+1] - X[i-1]);
-        p = sig * derivs2[i-1] + 2.0;
-        derivs2[i] = (sig - 1.0) / p;
-        u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1]));
-        u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p;
-
-        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
-                "X[%d] is %f\n", i, X[i]);
-        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
-                "Y[%d] is %f\n", i, Y[i]);
-        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
-                "u[%d] is %f\n", i, u[i]);
-    }
-
-    qn = 0.5;
-    u[n-1] = (3.0/(X[n-1]-X[n-2])) * (PS_RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2]));
-    derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0);
-
-    for (k=(n-2);k>=0;k--) {
-        derivs2[k] = derivs2[k] * derivs2[k+1] + u[k];
-
-        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
-                "derivs2[%d] is %f\n", k, derivs2[k]);
-    }
-
-    psFree(u);
-    psTrace(".psLib.dataManip.calculateSecondDerivs", 4,
-            "---- calculateSecondDerivs() end ----\n");
-    return(derivs2);
-}
-
-/******************************************************************************
-p_psNRSpline1DEval(): This routine does NR-style evaluation of cubic splines.
-It takes advantage of the 2nd derivatives of the cubic splines, which are
-stored in the psSPline1D data structure, and computes the interpolated value
-directly, without computing (or using) the interpolating cubic spline
-polynomial.
- 
-This routine is here mostly for a sanity check on the psLib function
-evalSpline() which computes the interpolated value based on the cubic spline
-polynomials which are stored in psSpline1D.
- 
-XXX: This is F32 only
- 
-XXX: spline->knots must be psF32
- 
-XXXX: Remove this for next code shipment.
- *****************************************************************************/
-/*
-psF32 p_psNRSpline1DEval(psSpline1D *spline,
-                         const psVector* x,
-                         const psVector* y,
-                         psF32 X)
-{
-    PS_ASSERT_PTR_NON_NULL(spline, NAN);
-    PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(spline->domains, NAN);
-    PS_ASSERT_PTR_NON_NULL(spline->p_psDeriv2, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
-    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(y, NAN);
-    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NAN);
-    PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NULL);
- 
-    psS32 n;
-    psS32 klo;
-    psS32 khi;
-    psF32 H;
-    psF32 A;
-    psF32 B;
-    psF32 C;
-    psF32 D;
-    psF32 Y;
- 
-    n = spline->n;
-    klo = p_psVectorBinDisect32(spline->knots->data.F32, (spline->n)+1, X);
-    if (klo < 0) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: psMinimize.c: p_psNRSpline1DEval(): p_psVectorBinDisect32 returned an error (%d).\n", klo);
-        return(NAN);
-    }
-    khi = klo + 1;
-    H = spline->knots->data.F32[khi] - spline->knots->data.F32[klo];
-    A = (spline->knots->data.F32[khi] - X) / H;
-    B = (X - spline->knots->data.F32[klo]) / H;
-    C = ((A*A*A)-A) * (H*H/6.0);
-    D = ((B*B*B)-B) * (H*H/6.0);
- 
-    Y = (A * y->data.F32[klo]) +
-        (B * y->data.F32[khi]) +
-        (C * (spline->p_psDeriv2)[klo]) +
-        (D * (spline->p_psDeriv2)[khi]);
- 
-    return(Y);
-}
-*/
-/*****************************************************************************/
-/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
-/*****************************************************************************/
-
-/*****************************************************************************
-psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y
-vectors, this routine generates the linear or cublic splines which satisfy
-those data points.
- 
-The formula for calculating the spline polynomials is derived from Numerical
-Recipes in C.  The basic idea is that the polynomial is
- (1)     y = (A * y[0]) +
- (2)         (B * y[1]) +
- (3)         ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 +
- (4)         ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0
-Where:
- H = x[1]-x[0]
- A = (x[1]-x)/H
- B = (x-x[0])/H
-The bulk of the code in this routine is the expansion of the above equation
-into a polynomial in terms of x, and then saving the coefficients of the
-powers of x in the spline polynomials.  This gets pretty complicated.
- 
-XXX: usage of yErr is not specified in IfA documentation.
- 
-XXX: Is the x argument redundant?  What do we do if the x argument is
-supplied, but does not equal the knots specified in mySpline?
- 
-XXX: can psSpline be NULL?
- 
-XXX: reimplement this assuming that mySpline is NULL?
- 
-XXX: What happens if X is NULL, then an index vector is generated for X, but
-that index vector lies outside the range vectors in mySpline?
- 
-XXX: Assumes mySpline->knots is psF32.  Must add psU32 and psF64.
- *****************************************************************************/
-psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,     ///< The spline which will be generated.
-                                const psVector* x,        ///< Ordinates (or NULL to just use the indices)
-                                const psVector* y,        ///< Coordinates
-                                const psVector* yErr)     ///< Errors in coordinates, or NULL
-{
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
-    if (mySpline != NULL) {
-        PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL);
-    }
-    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
-            "---- psVectorFitSpline1D() begin ----\n");
-    psS32 numSplines = (y->n)-1;
-    psF32 tmp;
-    psF32 H;
-    psS32 i;
-    psF32 slope;
-    psVector *x32 = NULL;
-    psVector *y32 = NULL;
-    psVector *yErr32 = NULL;
-    static psVector *x32Static = NULL;
-    static psVector *y32Static = NULL;
-    static psVector *yErr32Static = NULL;
-
-    PS_VECTOR_CONVERT_F64_TO_F32_STATIC(y, y32, y32Static);
-
-    // If yErr==NULL, set all errors equal.
-    if (yErr == NULL) {
-        PS_VECTOR_GEN_YERR_STATIC_F32(yErr32Static, y->n);
-        yErr32 = yErr32Static;
-    } else {
-        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(yErr, NULL);
-        PS_VECTOR_CONVERT_F64_TO_F32_STATIC(yErr, yErr32, yErr32Static);
-    }
-
-    // If x==NULL, create an x32 vector with x values set to (0:n).
-    if (x == NULL) {
-        PS_VECTOR_GEN_X_INDEX_STATIC_F32(x32Static, y->n);
-        x32 = x32Static;
-    } else {
-        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
-        PS_VECTOR_CONVERT_F64_TO_F32_STATIC(x, x32, x32Static);
-    }
-    PS_ASSERT_VECTORS_SIZE_EQUAL(x32, y32, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(yErr32, y32, NULL);
-
-    /*
-        XXX:
-        This can not be implemented until SDR states what order spline should be
-        created.
-        Should we error if mySpline is not NULL?
-        Should we error if mySPline is not NULL?
-    */
-    if (mySpline == NULL) {
-        mySpline = psSpline1DAllocGeneric(x32, 3);
-    }
-    PS_ASSERT_PTR_NON_NULL(mySpline, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(mySpline->n, NULL);
-
-    if (y32->n != (1 + mySpline->n)) {
-        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
-                "data size / spline size mismatch (%d %d)\n",
-                y32->n, mySpline->n);
-        return(NULL);
-    }
-
-    // If these are linear splines, which means their polynomials will have
-    // two coefficients, then we do the simple calculation.
-    if (2 == (mySpline->spline[0])->n) {
-        for (i=0;i<mySpline->n;i++) {
-            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
-                    (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]);
-            (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
-                                              (slope * mySpline->knots->data.F32[i]);
-
-            (mySpline->spline[i])->coeff[1] = slope;
-            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
-                    "---- mySpline %d coeffs are (%f, %f)\n", i,
-                    (mySpline->spline[i])->coeff[0],
-                    (mySpline->spline[i])->coeff[1]);
-        }
-        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
-                "---- Exiting psVectorFitSpline1D()()\n");
-        return((psSpline1D *) mySpline);
-    }
-
-    // Check if these are cubic splines (n==4).  If not, psError.
-    if (4 != (mySpline->spline[0])->n) {
-        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
-                "Don't know how to generate %d-order splines.",
-                (mySpline->spline[0])->n-1);
-        return(NULL);
-    }
-
-    // If we get here, then we know these are cubic splines.  We first
-    // generate the second derivatives at each data point.
-    mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
-    for (i=0;i<y32->n;i++)
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
-
-    // We generate the coefficients of the spline polynomials.  I can't
-    // concisely explain how this code works.  See above function comments
-    // and Numerical Recipes in C.
-    for (i=0;i<numSplines;i++) {
-        H = x32->data.F32[i+1] - x32->data.F32[i];
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
-                "x data (%f - %f) (%f)\n",
-                x32->data.F32[i],
-                x32->data.F32[i+1], H);
-        //
-        // ******** Calculate 0-order term ********
-        //
-        // From (1)
-        (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
-        // From (2)
-        ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
-        // From (3)
-        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
-        tmp-= (x32->data.F32[i+1] / H);
-        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[0])+= tmp;
-        // From (4)
-        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
-        tmp+= (x32->data.F32[i] / H);
-        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[0])+= tmp;
-
-        //
-        // ******** Calculate 1-order term ********
-        //
-        // From (1)
-        (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
-        // From (2)
-        ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
-        // From (3)
-        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
-        tmp+= (1.0 / H);
-        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[1])+= tmp;
-        // From (4)
-        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
-        tmp-= (1.0 / H);
-        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[1])+= tmp;
-
-        //
-        // ******** Calculate 2-order term ********
-        //
-        // From (3)
-        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
-        // From (4)
-        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
-
-        //
-        // ******** Calculate 3-order term ********
-        //
-        // From (3)
-        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
-        // From (4)
-        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
-
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
-
-    }
-
-    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
-            "---- psVectorFitSpline1D() end ----\n");
-    return(mySpline);
-}
-
-
-
-
-
-
-
-
-
-/*****************************************************************************
-psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y
- 
-xF32 and yF32 are internal psVectors which are used to hold the psF32 versions
-of the input data, if necessary.  xPtr and yPtr are pointers to either xF32 or
-the x argument.  All computation is done on xPtr and yPtr.  xF32 and yF32 will
-simply be psFree() at the end.
- 
-XXX: nKnots makes no sense.  This number is always equal to the size of the x
-an y vectors.
- 
-XXX: How do we specify the spline order?  For now, order=3.
- *****************************************************************************/
-#define PS_XXX_SPLINE_ORDER 3
-psSpline1D *psVectorFitSpline1DNEW(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
-                                   const psVector* y,        ///< Coordinates
-                                   int nKnots)
-{
-    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n");
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
-    //    PS_ASSERT_INT_EQUAL(y->n, nKnots);
-
-    //
-    // The following code ensures that xPtr points to a psF32 version of the
-    // ordinate data.
-    //
-    psVector *xF32 = NULL;
-    psVector *xPtr = NULL;
-    if (x != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
-        if (PS_TYPE_F64 == x->type.type) {
-            xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
-            for (psS32 i = 0 ; i < x->n ; i++) {
-                xF32->data.F32[i] = (psF32) x->data.F64[i];
-            }
-            xPtr = xF32;
-        } else if (PS_TYPE_F32 == x->type.type) {
-            xPtr = (psVector *) x;
-        } else {
-            printf("XXX: Gen Error message: x is wrong type.\n");
-            return(NULL);
-        }
-    } else {
-        // If x==NULL, create an x32 vector with x values set to (0:n).
-        xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
-        for (psS32 i = 0 ; i < x->n ; i++) {
-            xF32->data.F32[i] = (psF32) i;
-        }
-        xPtr = xF32;
-    }
-
-    //
-    // If y is of type psF64, then create a new vector yF32 and convert the
-    // y elements.  Regardless of y's type, we create a yPtr which will be
-    // used in the remainder of this function.
-    //
-    psVector *yF32 = NULL;
-    psVector *yPtr = NULL;
-    if (PS_TYPE_F64 == y->type.type) {
-        yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
-        for (psS32 i = 0 ; i < y->n ; i++) {
-            yF32->data.F32[i] = (psF32) y->data.F64[i];
-        }
-        yPtr = yF32;
-    } else {
-        yPtr = (psVector *) y;
-    }
-
-    psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER);
-
-    psS32 numSplines = nKnots - 1;
-    psF32 tmp;
-    psF32 H;
-    psS32 i;
-    psF32 slope;
-    // XXX: get rid of x32 and y32 (this is from old code)
-    psVector *x32 = xPtr;
-    psVector *y32 = yPtr;
-
-    // If these are linear splines, which means their polynomials will have
-    // two coefficients, then we do the simple calculation.
-    if (1 == PS_XXX_SPLINE_ORDER) {
-        for (i=0;i<mySpline->n;i++) {
-            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
-                    (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]);
-            (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
-                                              (slope * mySpline->knots->data.F32[i]);
-
-            (mySpline->spline[i])->coeff[1] = slope;
-            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
-                    "---- mySpline %d coeffs are (%f, %f)\n", i,
-                    (mySpline->spline[i])->coeff[0],
-                    (mySpline->spline[i])->coeff[1]);
-        }
-        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
-                "---- Exiting psVectorFitSpline1D()()\n");
-        return((psSpline1D *) mySpline);
-    }
-
-    //
-    // Check if these are cubic splines (n==4).  If not, psError.
-    //
-    if (3 != PS_XXX_SPLINE_ORDER) {
-        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
-                "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER);
-        return(NULL);
-    }
-
-    //
-    // If we get here, then we know these are cubic splines.  We first
-    // generate the second derivatives at each data point.
-    //
-    mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
-    for (i=0;i<y32->n;i++)
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
-
-    //
-    // We generate the coefficients of the spline polynomials.  I can't
-    // concisely explain how this code works.  See above function comments
-    // and Numerical Recipes in C.
-    //
-    for (i=0;i<numSplines;i++) {
-        H = x32->data.F32[i+1] - x32->data.F32[i];
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
-                "x data (%f - %f) (%f)\n",
-                x32->data.F32[i],
-                x32->data.F32[i+1], H);
-        //
-        // ******** Calculate 0-order term ********
-        //
-        // From (1)
-        (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
-        // From (2)
-        ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
-        // From (3)
-        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
-        tmp-= (x32->data.F32[i+1] / H);
-        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[0])+= tmp;
-        // From (4)
-        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
-        tmp+= (x32->data.F32[i] / H);
-        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[0])+= tmp;
-
-        //
-        // ******** Calculate 1-order term ********
-        //
-        // From (1)
-        (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
-        // From (2)
-        ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
-        // From (3)
-        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
-        tmp+= (1.0 / H);
-        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[1])+= tmp;
-        // From (4)
-        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
-        tmp-= (1.0 / H);
-        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[1])+= tmp;
-
-        //
-        // ******** Calculate 2-order term ********
-        //
-        // From (3)
-        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
-        // From (4)
-        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
-
-        //
-        // ******** Calculate 3-order term ********
-        //
-        // From (3)
-        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
-        // From (4)
-        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
-
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
-
-    }
-
-    psFree(xF32);
-    psFree(yF32);
-
-    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
-            "---- psVectorFitSpline1D() end ----\n");
-
-    return(mySpline);
-}
-
-
-
-
-//XXX: What's this for?
-psF64 p_psImageGetElementF64(psImage *a, int i, int j);
 /******************************************************************************
 psMinimizeLMChi2():  This routine will take an procedure which calculates
@@ -717 +132 @@
      will have to convert all F32 input vectors to F64 regardless.  So,
      the F64 port might be.
+ 
+XXX: Change the whole thing to F64, if input data is F32, convert it.
  *****************************************************************************/
 psBool psMinimizeLMChi2(psMinimization *min,
@@ -816 +233 @@
         if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) {
             // XXX: Where is this?
-            //            p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
+            // p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
         }
         # endif /* PS_NO_TRACE */
@@ -876 +293 @@
 
 // XXX EAM: this needs to respect the mask on params
-// XXX EAM: check not NULL on alpha, beta, params
 // alpha, beta, params are already allocated
 psF64 p_psMinLM_SetABX (psImage  *alpha,
@@ -887 +303 @@
                         psMinimizeLMChi2Func func)
 {
+    PS_ASSERT_IMAGE_NON_NULL(alpha, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(beta, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(x, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(y, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(dy, NAN);
 
     psF64 chisq;
@@ -1011 +433 @@
 }
 
-// XXX: put this in constants.h
-# define PS_SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;}
-
 // XXX EAM : temporary gauss-jordan solver based on gene's
 // version based on the Numerical Recipes version
@@ -1044 +463 @@
         for (j = 0; j < Nx; j++) {
             if (!isfinite(matrix[i][j])) {
-                // XXX EAM: this should use the psError stack
-                fprintf (stderr, "GAUSSJ: NaN\n");
+                psError(PS_ERR_UNKNOWN, false, "Input matrix contains NaNs.\n");
                 goto fescape;
             }
@@ -1058 +476 @@
                     } else {
                         if (ipiv[k] > 1) {
-                            // XXX EAM: this should use the psError stack
-                            fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n");
+                            psError(PS_ERR_UNKNOWN, false, "Singular Matrix (1).\n");
                             goto fescape;
                         }
@@ -1076 +493 @@
         indxc[i] = icol;
         if (matrix[icol][icol] == 0.0) {
-            // XXX EAM: this should use the psError stack
-            fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n");
+            psError(PS_ERR_UNKNOWN, false, "Singular Matrix (2).\n");
             goto fescape;
         }
@@ -1118 +534 @@
 }
 
-/******************************************************************************
-vectorFitPolynomial1DCheb():  This routine will fit a Chebyshev polynomial of
-degree myPoly to the data points (x, y) and return the coefficients of that
-polynomial.
- 
-XXX: yErr is currently ignored.
-*****************************************************************************/
-static psPolynomial1D *vectorFitPolynomial1DCheby(psPolynomial1D* myPoly,
-        const psVector* x,
-        const psVector* y,
-        const psVector* yErr)
-{
-    psS32 j;
-    psS32 k;
-    psS32 n = x->n;
-    psF64 fac;
-    psF64 sum;
-    PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64);
-    psScalar *fScalar;
-    psScalar tmpScalar;
-    tmpScalar.type.type = PS_TYPE_F64;
-
-    // XXX: These assignments appear too simple to warrant code and
-    // variable declarations.  I retain them here to maintain coherence
-    // with the NR code.
-    psF64 min = -1.0;
-    psF64 max = 1.0;
-    psF64 bma = 0.5 * (max-min);  // 1
-    psF64 bpa = 0.5 * (max+min);  // 0
-
-    // In this loop, we first calculate the values of X for which the
-    // Chebyshev polynomials are zero (see NR, section 5.4).  Then we
-    // calculate the value of the function we are fitting the Chebyshev
-    // polynomials to at those values of X.  This is a bit tricky since
-    // we don't know that function.  So, we instead do 3-order LaGrange
-    // interpolation at the point X for the psVectors x,y for which we
-    // are fitting this ChebyShev polynomial to.
-
-    for (psS32 i=0;i<n;i++) {
-        // NR 5.8.4
-        psF64 Y = cos(M_PI * (0.5 + ((psF32) i)) / ((psF32) n));
-        psF64 X = (Y + bma + bpa) - 1.0;
-        tmpScalar.data.F64 = X;
-
-        // We interpolate against are tabluated x,y vectors to determine the
-        // function value at X.
-        fScalar = p_psVectorInterpolate((psVector *) x,
-                                        (psVector *) y,
-                                        3,
-                                        &tmpScalar);
-
-        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",
-                x->data.F64[i], X, y->data.F64[i], f->data.F64[i]);
-    }
-
-    // We have the values for f() at the zero points, we now calculate the
-    // coefficients of the Chebyshev polynomial: NR 5.8.7.
-
-    fac = 2.0/((psF32) n);
-    // XXX: is this loop bound correct?
-    for (j=0;j<myPoly->n;j++) {
-        sum = 0.0;
-        for (k=0;k<n;k++) {
-            sum+= f->data.F64[k] *
-                  cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n));
-        }
-
-        myPoly->coeff[j] = fac * sum;
-    }
-
-    return(myPoly);
-}
-
-/******************************************************************************
-VectorFitPolynomial1DOrd():  This routine will fit an ordinary polynomial of
-degree myPoly to the data points (x, y) and return the coefficients of that
-polynomial.
- 
-XXX: Use private name?
-XXX: Use recycled vectors.
- *****************************************************************************/
-static psPolynomial1D* vectorFitPolynomial1DOrd(psPolynomial1D* myPoly,
-        const psVector* x,
-        const psVector* y,
-        const psVector* yErr)
-{
-    psS32 polyOrder = myPoly->n;
-    psImage* A = NULL;
-    psImage* ALUD = NULL;
-    psVector* B = NULL;
-    psVector* outPerm = NULL;
-    psVector* X = NULL;         // NOTE: do we need this?
-    psVector* coeffs = NULL;
-    psS32 i = 0;
-    psS32 j = 0;
-    psS32 k = 0;
-    psVector* xSums = NULL;
-
-    psTrace(".psLib.dataManip.vectorFitPolynomial1DOrd", 4,
-            "---- vectorFitPolynomial1DOrd() begin ----\n");
-    // printf("VectorFitPolynomial1D()\n");
-    // for (i=0;i<x->n;i++) {
-    // printf("(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
-    // }
-
-    A = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64);
-    ALUD = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64);
-
-    B = psVectorAlloc(polyOrder, PS_TYPE_F64);
-    coeffs = psVectorAlloc(polyOrder, PS_TYPE_F64);
-    X = psVectorAlloc(x->n, PS_TYPE_F64);
-    xSums = psVectorAlloc(1 + 2 * polyOrder, PS_TYPE_F64);
-
-    // Initialize data structures.
-    for (i = 0; i < polyOrder; i++) {
-        B->data.F64[i] = 0.0;
-        coeffs->data.F64[i] = 0.0;
-        for (j = 0; j < polyOrder; j++) {
-            A->data.F64[i][j] = 0.0;
-            ALUD->data.F64[i][j] = 0.0;
-        }
-    }
-    for (i = 0; i < X->n; i++) {
-        X->data.F64[i] = x->data.F64[i];
-    }
-
-    // Build the B and A data structs.
-    if (yErr == NULL) {
-        for (i = 0; i < X->n; i++) {
-            buildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
-
-            for (k = 0; k < polyOrder; k++) {
-                B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k];
-            }
-
-            for (k = 0; k < polyOrder; k++) {
-                for (j = 0; j < polyOrder; j++) {
-                    A->data.F64[k][j] += xSums->data.F64[k + j];
-                }
-            }
-        }
-    } else {
-        for (i = 0; i < X->n; i++) {
-            buildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
-
-            for (k = 0; k < polyOrder; k++) {
-                B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k] /
-                                  yErr->data.F64[i];
-            }
-
-            for (k = 0; k < polyOrder; k++) {
-                for (j = 0; j < polyOrder; j++) {
-                    A->data.F64[k][j] += xSums->data.F64[k + j] /
-                                         yErr->data.F64[i];
-                }
-            }
-        }
-    }
-
-    // XXX: How do we know if these routines were successful?
-    ALUD = psMatrixLUD(ALUD, &outPerm, A);
-    coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-
-    for (k = 0; k < polyOrder; k++) {
-        myPoly->coeff[k] = coeffs->data.F64[k];
-        // printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]);
-    }
-
-    psFree(A);
-    psFree(ALUD);
-    psFree(B);
-    psFree(coeffs);
-    psFree(X);
-    psFree(outPerm);
-    psFree(xSums);
-
-    psTrace(".psLib.dataManip.vectorFitPolynomial1DOrd", 4,
-            "---- vectorFitPolynomial1DOrd() begin ----\n");
-    return (myPoly);
-}
-
-/******************************************************************************
-psVectorFitPolynomial1D():  This routine must fit a polynomial of degree
-myPoly to the data points (x, y) and return the coefficients of that
-polynomial.
- 
-XXX: type F32 is done via vector conversion only.
- 
-XXX: Get rid of this.  Use new argument list.
- *****************************************************************************/
-psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* poly,
-                                        const psVector* x,
-                                        const psVector* y,
-                                        const psVector* yErr)
-{
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(poly->n, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTOR_NON_EMPTY(y, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
-
-    psS32 i;
-    psVector *x64 = NULL;
-    psVector *y64 = NULL;
-    psVector *yErr64 = NULL;
-    static psVector *x64Static = NULL;
-    static psVector *y64Static = NULL;
-    static psVector *yErr64Static = NULL;
-
-    PS_VECTOR_CONVERT_F32_TO_F64_STATIC(y, y64, y64Static);
-    // If yErr==NULL, set all errors equal.
-    if (yErr == NULL) {
-        PS_VECTOR_GEN_YERR_STATIC_F64(yErr64Static, y->n);
-        yErr64 = yErr64Static;
-    } else {
-        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(yErr, NULL);
-        PS_VECTOR_CONVERT_F32_TO_F64_STATIC(yErr, yErr64, yErr64Static);
-    }
-
-    // If x==NULL, create an x64 vector with x values set to (0:n).
-    if (x == NULL) {
-        PS_VECTOR_GEN_X_INDEX_STATIC_F64(x64Static, y->n);
-        if (poly->type == PS_POLYNOMIAL_CHEB) {
-            p_psNormalizeVectorRangeF64(x64Static, -1.0, 1.0);
-        }
-        x64 = x64Static;
-    } else {
-        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
-        PS_VECTOR_CONVERT_F32_TO_F64_STATIC(x, x64, x64Static);
-        if (poly->type == PS_POLYNOMIAL_CHEB) {
-            p_psNormalizeVectorRangeF64(x64, -1.0, 1.0);
-        }
-    }
-    PS_ASSERT_VECTORS_SIZE_EQUAL(x64, y64, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(yErr64, y64, NULL);
-
-    // Call the appropriate vector fitting routine.
-    psPolynomial1D *rc = NULL;
-    if (poly->type == PS_POLYNOMIAL_CHEB) {
-        rc = vectorFitPolynomial1DCheby(poly, x64, y64, yErr64);
-    } else if (poly->type == PS_POLYNOMIAL_ORD) {
-        rc = vectorFitPolynomial1DOrd(poly, x64, y64, yErr64);
-    } else {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                "unknown polynomial type.\n");
-        return(NULL);
-    }
-    if (rc == NULL) {
-        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
-        return(NULL);
-    }
-
-    return(poly);
-}
-
 static void minimizationFree(psMinimization *min)
 {
-    // There are non dynamic allocated items
+    // There are no dynamically allocated items
 }
 
@@ -1437 +594 @@
 
 /******************************************************************************
+ ******************************************************************************
+Powell routines.
+******************************************************************************
+*****************************************************************************/
+
+/******************************************************************************
 p_psDetermineBracket():  This routine takes as input an arbitrary function,
 and the parameter to vary, and the line along which it must vary.  This
@@ -1446 +609 @@
 Algorithm:
  
-XXX completely ad hoc:
-start with the user-supplied starting parameter and
+XXX completely ad hoc: start with the user-supplied starting parameter and
 call that b.  Calculate a/c as a fractional amount smaller/larger than b.
 Repeat this process until a local minimum is found.
  
-XXX:
-new algorithm:
-start at x=0, expand in one direction until the function
+XXX: new algorithm:  start at x=0, expand in one direction until the function
 decreases.  Then you have two points in the bracket.  Keep going until it
 increases, or x is too large.  If thst does not work, expand in the other
 direction.
  
-XXX:
-This is F32 only.
+XXX: This is F32 only.  Must add F64 support (actually, make the defaults F64,
+and convert F32 vectors to F64).
  
-XXX:
-output bracket vector should be an input as well.
+XXX: output bracket vector should be an input as well.
 *****************************************************************************/
 psVector *p_psDetermineBracket(psVector *params,
@@ -1716 +875 @@
 XXX: This routine is not very efficient in terms of total evaluations of the
 function.
-XXX: This is F32 only
+XXX: This is F32 only (make it F64).
 XXX: Since this is an internal function, many of the parameter checks are
      redundant.
@@ -1871 +1030 @@
  
 XXX: The SDR is silent about data types.  F32 is implemented here.
- 
-XXX: Check for F32 types?
+Reimplement in F64, convert F32 vectors to F64.
  *****************************************************************************/
 #define PS_MINIMIZE_POWELL_LINEMIN_MAX_ITERATIONS 20
@@ -2083 +1241 @@
 This functions uses global variables to receive the function pointer, the
 data values, and the data errors.
-XXX: This is F32 only
+XXX: This is F32 only.  Must implement F64.
  *****************************************************************************/
 static psF32 myPowellChi2Func(const psVector *params,
@@ -2149 +1307 @@
 /******************************************************************************
  ******************************************************************************
- ******************************************************************************
- ******************************************************************************
-EAM Code:
- ******************************************************************************
- ******************************************************************************
+ Analytical 1-D fitting routines.
  ******************************************************************************
  *****************************************************************************/
-
-static psVector *EAMBuildSums1D(
+// 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]; \
+} \
+
+#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; \
+}\
+/******************************************************************************
+BuildSums1D(sums, x, polyOrder, sums): this routine calculates the powers of
+input parameter "x" between 0 and input parameter nTerms*2.  The result is
+returned as a psVector sums.
+*****************************************************************************/
+static psVector *BuildSums1D(
     psVector* sums,
     psF64 x,
@@ -2166 +1336 @@
 
     //
-    // XXX: Why do we multiply by 2 here?  Better to do it outside and have the
-    // definition of this function remain sensible.
+    // XXX: Why do we multiply by 2 here?  It's better to do it outside and
+    // have the definition of this function remain sensible.
     //
     nSum = 2*nTerm;
@@ -2184 +1354 @@
 }
 
-// XXX EAM : test version of 1d fitting
-psPolynomial1D* VectorFitPolynomial1DOrd_EAM(
+/******************************************************************************
+BuildSums2D(sums, x, y, nXterm, nYterm): this routine calculates the powers of
+input parameter "x" and "y" between 0 and input parameter nXterms*2 and
+nYterm*2.  The result is returned as a psImage sums.
+ *****************************************************************************/
+static psImage *BuildSums2D(
+    psImage *sums,
+    psF64 x,
+    psF64 y,
+    psS32 nXterm,
+    psS32 nYterm)
+{
+    psS32 nXsum = 0;
+    psS32 nYsum = 0;
+    psF64 xSum = 1.0;
+    psF64 ySum = 1.0;
+
+    nXsum = 2*nXterm;
+    nYsum = 2*nYterm;
+    if (sums == NULL) {
+        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
+    }
+    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
+        psFree (sums);
+        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
+    }
+
+    ySum = 1.0;
+    for (int j = 0; j < nYsum; j++) {
+        xSum = ySum;
+        for (int i = 0; i < nXsum; i++) {
+            sums->data.F64[i][j] = xSum;
+            xSum *= x;
+        }
+        ySum *= y;
+    }
+    return (sums);
+}
+
+
+/******************************************************************************
+Polynomial2DEvalVectorD(myPoly, x, y): This routine takes as input two
+psVectors x and y, and evaluates myPoly for each pair of (x, y), and stores it
+in the output vector.  This routine works on single-precision polynomials with
+double precision data.
+ *****************************************************************************/
+psVector *Polynomial2DEvalVectorD(
+    const psPolynomial2D *myPoly,
+    const psVector *x,
+    const psVector *y)
+{
+    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
+
+    //
+    // Set the length of the output vector to the minimum of the x,y vectors.
+    //
+    psS32 vecLen=x->n;
+    if (y->n < vecLen) {
+        vecLen = y->n;
+    }
+
+    //
+    // Create output vector to return
+    //
+    psVector *tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
+
+    //
+    // Evaluate the polynomial at the specified points
+    //
+    for (psS32 i=0; i<vecLen; i++) {
+        tmp->data.F64[i] = psPolynomial2DEval(myPoly, x->data.F64[i], y->data.F64[i]);
+    }
+
+    return(tmp);
+}
+
+/******************************************************************************
+ ******************************************************************************
+ 1-D Vector Fitting Code.
+ ******************************************************************************
+ *****************************************************************************/
+
+/******************************************************************************
+vectorFitPolynomial1DCheb():  This routine will fit a Chebyshev polynomial of
+degree myPoly to the data points (x, y) and return the coefficients of that
+polynomial.
+ 
+XXX: mask, maskValue, yErr are currently ignored.
+ 
+XXX: Change arg order to that of the psLib function.
+*****************************************************************************/
+static psPolynomial1D *vectorFitPolynomial1DCheby(
     psPolynomial1D* myPoly,
-    psVector *mask,
-    const psVector *x,
-    const psVector *y,
-    const psVector *yErr)
-{
-    // I think this is 1 dimension down
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector* y,
+    const psVector* yErr,
+    const psVector* x)
+{
+    // XXX: these ASSERTS are redundant.
+    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->n, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
+    if (yErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, yErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(yErr, PS_TYPE_F64, NULL);
+    }
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, NULL);
+        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
+    }
+
+    psS32 j;
+    psS32 k;
+    psS32 n = y->n;
+    psF64 fac;
+    psF64 sum;
+    PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64);
+    psScalar *fScalar;
+    psScalar tmpScalar;
+    tmpScalar.type.type = PS_TYPE_F64;
+
+    // XXX: These assignments appear too simple to warrant code and
+    // variable declarations.  I retain them here to maintain coherence
+    // with the NR code.
+    psF64 min = -1.0;
+    psF64 max = 1.0;
+    psF64 bma = 0.5 * (max-min);  // 1
+    psF64 bpa = 0.5 * (max+min);  // 0
+
+    // In this loop, we first calculate the values of X for which the
+    // Chebyshev polynomials are zero (see NR, section 5.4).  Then we
+    // calculate the value of the function we are fitting the Chebyshev
+    // polynomials to at those values of X.  This is a bit tricky since
+    // we don't know that function.  So, we instead do 3-order LaGrange
+    // interpolation at the point X for the psVectors x,y for which we
+    // are fitting this ChebyShev polynomial to.
+
+    for (psS32 i=0;i<n;i++) {
+        // NR 5.8.4
+        psF64 Y = cos(M_PI * (0.5 + ((psF32) i)) / ((psF32) n));
+        psF64 X = (Y + bma + bpa) - 1.0;
+        tmpScalar.data.F64 = X;
+
+        // We interpolate against the tabluated x,y vectors to determine the
+        // function value at X.
+        fScalar = p_psVectorInterpolate((psVector *) x,
+                                        (psVector *) y,
+                                        3,
+                                        &tmpScalar);
+
+        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",
+                x->data.F64[i], X, y->data.F64[i], f->data.F64[i]);
+    }
+
+    // We have the values for f() at the zero points, we now calculate the
+    // coefficients of the Chebyshev polynomial: NR 5.8.7.
+
+    fac = 2.0/((psF32) n);
+    // XXX: is this loop bound correct?
+    for (j=0;j<myPoly->n;j++) {
+        sum = 0.0;
+        for (k=0;k<n;k++) {
+            sum+= f->data.F64[k] *
+                  cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n));
+        }
+
+        myPoly->coeff[j] = fac * sum;
+    }
+
+    return(myPoly);
+}
+
+/******************************************************************************
+VectorFitPolynomial1DOrd(myPoly, *mask, maskValue, *y, *yErr, *x): This is a
+private routine which will fit a 1-D polynomial to a set of (x, f) pairs.  The
+x and fErr vectors may be NULL.  All non-NULL vectors must be of type
+PS_TYPE_F64.
+ *****************************************************************************/
+psPolynomial1D* VectorFitPolynomial1DOrd(
+    psPolynomial1D* myPoly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x)
+{
+    // XXX: these ASSERTS are redundant.
+    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->n, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL);
+    }
+
     psImage*     A = NULL;
     psVector*    B = NULL;
@@ -2204 +1579 @@
 
     if (psTraceGetLevel (".psLib.dataManip.VectorFitPolynomial1DOrd") >= 5) {
-        FILE *f = fdopen((psTraceGetDestination ()), "a+");
-        fprintf (f, "VectorFitPolynomial1D()\n");
-        for (int i = 0; i < x->n; i++) {
-            fprintf (f, "(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
-        }
-        fclose(f);
+        FILE *fp = fdopen((psTraceGetDestination ()), "a+");
+        fprintf(fp, "VectorFitPolynomial1D()\n");
+        for (int i = 0; i < f->n; i++) {
+            fprintf(fp, "(x, f, fErr) is (");
+            if (x != NULL) {
+                fprintf(fp, "%f, %f, ", x->data.F64[i], f->data.F64[i]);
+            } else {
+                fprintf(fp, "%f, %f, ", (psF64) i, f->data.F64[i]);
+            }
+            if (fErr != NULL) {
+                fprintf(fp, "%f)\n", fErr->data.F64[i]);
+            } else {
+                fprintf(fp, "NULL)\n");
+            }
+        }
+        fclose(fp);
     }
 
@@ -2227 +1612 @@
     // xSums look like: 1, x, x^2, ... x^(2n+1)
     // Build the B and A data structs.
-    for (int k = 0; k < x->n; k++) {
+    for (int k = 0; k < f->n; k++) {
         if ((mask != NULL) && mask->data.U8[k])
             continue;
-        xSums = EAMBuildSums1D(xSums, x->data.F64[k], nTerm);
-
-        if (yErr == NULL) {
+        if (x != NULL) {
+            xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
+        } else {
+            xSums = BuildSums1D(xSums, (psF64) k, nTerm);
+        }
+
+        if (fErr == NULL) {
             wt = 1.0;
         } else {
-            // this should filter yErr == 0 values
-            wt = 1.0 / PS_SQR(yErr->data.F64[k]);
+            // this should filter fErr == 0 values
+            wt = 1.0 / PS_SQR(fErr->data.F64[k]);
         }
         for (int i = 0; i < nTerm; i++) {
-            B->data.F64[i] += y->data.F64[k] * xSums->data.F64[i] * wt;
+            B->data.F64[i] += f->data.F64[k] * xSums->data.F64[i] * wt;
         }
 
@@ -2261 +1650 @@
             myPoly->coeff[k] = B->data.F64[k];
         }
-    } else
+    } else {
         // LUD version of the fit
-    {
         psImage *ALUD = NULL;
         psVector* outPerm = NULL;
@@ -2274 +1662 @@
             myPoly->coeff[k] = coeffs->data.F64[k];
         }
+        psFree(ALUD);
+        psFree(coeffs);
+        psFree(outPerm);
     }
 
@@ -2281 +1672 @@
 
     psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4,
-            "---- VectorFitPolynomial1DOrd() begin ----\n");
+            "---- VectorFitPolynomial1DOrd() End ----\n");
     return (myPoly);
 }
 
-// XXX EAM : this version uses the F64 vectors
-psVector *Polynomial2DEvalVectorD(
-    const psPolynomial2D *myPoly,
+
+/******************************************************************************
+psVectorFitPolynomial1D():  This routine fits a polynomial of arbitrary degree
+(specified in poly) to the data points (x, y) and return that polynomial.
+Types F32 and F64 are supported, however, type F32 is done via vector
+conversion only.
+ *****************************************************************************/
+psPolynomial1D *psVectorFitPolynomial1D(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x)
+{
+    // Internal pointers for possibly NULL or mis-typed vectors.
+    psVector *x64 = NULL;
+    psVector *f64 = NULL;
+    psVector *fErr64 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(poly->n, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_NON_EMPTY(f, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
+    if (f->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+    } else {
+        f64 = (psVector *) f;
+    }
+
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
+        if (x->type.type != PS_TYPE_F64) {
+            PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+        } else {
+            x64 = (psVector *) x;
+        }
+    }
+
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
+        if (fErr->type.type != PS_TYPE_F64) {
+            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+        } else {
+            fErr64 = (psVector *) fErr;
+        }
+    }
+
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        poly = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
+        if (poly == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
+            return(NULL);
+        }
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        if (mask != NULL) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
+        }
+        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);
+        }
+        // XXX: Change arg order.
+        // XXX: Must modify this routine so that x64 or fErr64 can be NULL.
+        poly = vectorFitPolynomial1DCheby(poly, NULL, 0, f64, fErr64, x64);
+        if (x == NULL) {
+            psFree(x64);
+        }
+    } else {
+        psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type.  Returning NULL.\n");
+    }
+
+    // Free psVectors that were created for NULL arguments.
+    if (f->type.type != PS_TYPE_F64) {
+        psFree(f64);
+    }
+
+    if ((x != NULL) && (x->type.type != PS_TYPE_F64)) {
+        psFree(x64);
+    }
+
+    if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
+        psFree(fErr64);
+    }
+
+    return(NULL);
+}
+
+
+psPolynomial1D *psVectorClipFitPolynomial1D(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *x32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_PTR_NON_NULL(stats, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    }
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    // Free psVectors that were created for NULL arguments.
+    if (x == NULL) {
+        psFree(x32);
+    }
+    return(NULL);
+}
+
+
+/******************************************************************************
+ ******************************************************************************
+ 2-D Vector Fitting Code.
+ ******************************************************************************
+ *****************************************************************************/
+
+/******************************************************************************
+VectorFitPolynomial2DOrd(myPoly, *mask, maskValue, *f, *fErr, *x, *y): This is
+a private routine which will fit a 2-D polynomial to a set of (x, y)-(f)
+pairs.  All non-NULL vectors must be of type PS_TYPE_F64.
+ 
+// XXX: What about the maskValue?
+ *****************************************************************************/
+psPolynomial2D* VectorFitPolynomial2DOrd(
+    psPolynomial2D* myPoly,
+    const psVector* mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
     const psVector *x,
     const psVector *y)
 {
+    // These ASSERTS are redundant.
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL);
+
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
     PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
-
-    psVector *tmp;
-    psS32 vecLen=x->n;
-
-    // Determine the length of the output vector to by the minimum of the x,y vectors
-    if (y->n < vecLen) {
-        vecLen = y->n;
-    }
-
-    // Create output vector to return
-    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
-
-    // Evaluate the polynomial at the specified points
-    for (psS32 i=0; i<vecLen; i++) {
-        tmp->data.F64[i] = psPolynomial2DEval(myPoly, x->data.F64[i], y->data.F64[i]);
-    }
-
-    // Return output vector
-    return(tmp);
-}
-
-// XXX EAM : EAMBuildSums2D in analogy with EAMBuildSums1D
-static psImage *EAMBuildSums2D(
-    psImage *sums,
-    psF64 x,
-    psF64 y,
-    psS32 nXterm,
-    psS32 nYterm)
-{
-    psS32 nXsum = 0;
-    psS32 nYsum = 0;
-    psF64 xSum = 1.0;
-    psF64 ySum = 1.0;
-
-    nXsum = 2*nXterm;
-    nYsum = 2*nYterm;
-    if (sums == NULL) {
-        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
-    }
-    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
-        psFree (sums);
-        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
-    }
-
-    ySum = 1.0;
-    for (int j = 0; j < nYsum; j++) {
-        xSum = ySum;
-        for (int i = 0; i < nXsum; i++) {
-            sums->data.F64[i][j] = xSum;
-            xSum *= x;
-        }
-        ySum *= y;
-    }
-    return (sums);
-}
-
-// XXX EAM : test version of 2d fitting
-psPolynomial2D* VectorFitPolynomial2DOrd_EAM(
-    psPolynomial2D* myPoly,
-    psVector* mask,
-    const psVector* x,
-    const psVector* y,
-    const psVector* z,
-    const psVector* zErr)
-{
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+
     // I think this is 1 dimension down
     psImage*     A = NULL;
@@ -2389 +1884 @@
         if ((mask != NULL) && mask->data.U8[k])
             continue;
-        Sums = EAMBuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
-
-        if (zErr == NULL) {
+        Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
+
+        if (fErr == NULL) {
             wt = 1.0;
         } else {
-            // XXX: this should probably by zErr^2 !!
-            // this should filter zErr == 0 values
-            // XXX: Why isn't this zErr^2?
-            wt = 1.0 / zErr->data.F64[k];
+            // XXX: this should probably by fErr^2 !!
+            // this should filter fErr == 0 values
+            // XXX: Why isn't this fErr^2?
+            wt = 1.0 / fErr->data.F64[k];
         }
 
@@ -2404 +1899 @@
         for (int n = 0; n < nXterm; n++) {
             for (int m = 0; m < nYterm; m++) {
-                B->data.F64[n+m*nXterm] += z->data.F64[k] * Sums->data.F64[n][m] * wt;
+                B->data.F64[n+m*nXterm] += f->data.F64[k] * Sums->data.F64[n][m] * wt;
             }
         }
@@ -2447 +1942 @@
     const psVector* x,
     const psVector* y,
-    const psVector* z,
-    const psVector* dz)
+    const psVector* f,
+    const psVector* df)
 {
     psVector *X;
@@ -2455 +1950 @@
     psVector *dZ;
 
-    psVector *zFit   = NULL;
-    psVector *zResid = NULL;
+    psVector *fFit   = NULL;
+    psVector *fResid = NULL;
     psStats  *stats  = NULL;
 
     X  = psVectorCopy (NULL, x, PS_TYPE_F64);
     Y  = psVectorCopy (NULL, y, PS_TYPE_F64);
-    Z  = psVectorCopy (NULL, z, PS_TYPE_F64);
-    dZ = psVectorCopy (NULL, dz, PS_TYPE_F64);
+    Z  = psVectorCopy (NULL, f, PS_TYPE_F64);
+    dZ = psVectorCopy (NULL, df, PS_TYPE_F64);
 
     for (int N = 0; N < 3; N++) {
         // XXX EAM : this would be better defined with an element mask
-        poly   = VectorFitPolynomial2DOrd_EAM(poly, NULL, X, Y, Z, dZ);
-        zFit   = Polynomial2DEvalVectorD(poly, x, y);
-        zResid = (psVector *) psBinaryOp(NULL, (void *) z, "-", (void *) zFit);
+        poly   = VectorFitPolynomial2DOrd(poly, NULL, 0, X, Y, Z, dZ);
+        fFit   = Polynomial2DEvalVectorD(poly, x, y);
+        fResid = (psVector *) psBinaryOp(NULL, (void *) f, "-", (void *) fFit);
 
         stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-        stats  = psVectorStats (stats, zResid, NULL, NULL, 0);
+        stats  = psVectorStats (stats, fResid, NULL, NULL, 0);
         psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g (%d pts)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues);
 
         // re-create X, Y, Z, dZ if pts are valid
         int n = 0;
-        for (int i = 0; i < zResid->n; i++) {
-            if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
+        for (int i = 0; i < fResid->n; i++) {
+            if (fabs(fResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
                 continue;
             }
             X->data.F64[n]  =  x->data.F64[i];
             Y->data.F64[n]  =  y->data.F64[i];
-            Z->data.F64[n]  =  z->data.F64[i];
-            dZ->data.F64[n] = dz->data.F64[i];
+            Z->data.F64[n]  =  f->data.F64[i];
+            dZ->data.F64[n] = df->data.F64[i];
             n++;
         }
@@ -2507 +2002 @@
 
 psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
-                            psVector* mask,
+                            const psVector* mask,
                             const psVector* x,
                             const psVector* y,
-                            const psVector* z,
-                            const psVector* dz)
+                            const psVector* f,
+                            const psVector* df)
 {
     PS_ASSERT_VECTOR_NON_NULL(mask, NULL);
     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(dz, NULL);
-
-    psVector *zFit   = NULL;
-    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(df, NULL);
+
+    psVector *fFit   = NULL;
+    psVector *fResid = psVectorAlloc (x->n, PS_TYPE_F64);
     psStats  *stats  = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
 
     for (int N = 0; N < 3; N++) {
-        poly   = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz);
-        zFit   = Polynomial2DEvalVectorD (poly, x, y);
-        zResid = (psVector *) psBinaryOp(zResid, (void *) z, "-", (void *) zFit);
-
-        stats  = psVectorStats (stats, zResid, NULL, mask, 1);
+        poly   = VectorFitPolynomial2DOrd (poly, mask, 0, x, y, f, df);
+        fFit   = Polynomial2DEvalVectorD (poly, x, y);
+        fResid = (psVector *) psBinaryOp(fResid, (void *) f, "-", (void *) fFit);
+
+        stats  = psVectorStats (stats, fResid, NULL, mask, 1);
         psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n",
                  stats->sampleMean, stats->sampleStdev);
@@ -2535 +2030 @@
         // we are masking out any point which is out of range
         // recovery is not allowed with this scheme
-        for (int i = 0; i < zResid->n; i++) {
+        for (int i = 0; i < fResid->n; i++) {
             if (mask->data.U8[i])
                 continue;
-            if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
+            if (fabs(fResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
                 mask->data.U8[i] = 1;
                 continue;
             }
         }
-        psFree (zFit);
-    }
-    psFree (zResid);
+        psFree (fFit);
+    }
+    psFree (fResid);
     psFree (stats);
     return (poly);
 }
 
+
 /******************************************************************************
- ******************************************************************************
- ******************************************************************************
- ******************************************************************************
-NEW Code:
- ******************************************************************************
- ******************************************************************************
- ******************************************************************************
+psVectorFitPolynomial2D():  This routine fits a 2D polynomial of arbitrary
+degree (specified in poly) to the data points (x, y)-(f) and returns that
+polynomial.  Types F32 and F64 are supported, however, type F32 is done via
+vector conversion only.
  *****************************************************************************/
-
-#define PS_VECTOR_GEN_INDEX_F32(VEC, SIZE) \
-VEC = psVectorAlloc(SIZE, PS_TYPE_F32); \
-for (psS32 i = 0 ; i < SIZE ; i++) { \
-    VEC->data.F32[i] = (psF32) i; \
-} \
-
-psPolynomial1D *psVectorFitPolynomial1D_NEW(
-    psPolynomial1D *poly,
-    const psVector *mask,
-    psMaskType maskValue,
-    const psVector *f,
-    const psVector *fErr,
-    const psVector *x)
-{
-    // Internal pointers for possibly NULL vectors.
-    psVector *x32 = NULL;
-    psVector *fErr32 = NULL;
-
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    if (x != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
-    } else {
-        PS_VECTOR_GEN_INDEX_F32(x32, f->n);
-    }
-
-    if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
-    } else {
-        fErr = psVectorAlloc(f->n, PS_TYPE_F32);
-        PS_VECTOR_SET_F32(fErr, 1.0);
-    }
-
-    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
-    if (poly->type == PS_POLYNOMIAL_ORD) {
-        // XXX: Call EAM code
-    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
-        // XXX: Call my code
-    } else {
-        printf("XXX: ERROR: incorrect polynomial type.\n");
-    }
-
-    // Free psVectors that were created for NULL arguments.
-    if (x == NULL) {
-        psFree(x32);
-    }
-    if (fErr == NULL) {
-        psFree(fErr32);
-    }
-
-    return(NULL);
-}
-
-
-psPolynomial2D *psVectorFitPolynomial2D_NEW(
-    psPolynomial1D *poly,
+psPolynomial2D *psVectorFitPolynomial2D(
+    psPolynomial2D *poly,
     const psVector *mask,
     psMaskType maskValue,
@@ -2630 +2061 @@
     const psVector *y)
 {
-    // Internal pointers for possibly NULL vectors.
-    psVector *fErr32 = NULL;
+    // Internal pointers for possibly NULL or mis-typed vectors.
+    psVector *x64 = NULL;
+    psVector *y64 = NULL;
+    psVector *f64 = NULL;
+    psVector *fErr64 = NULL;
 
     PS_ASSERT_POLY_NON_NULL(poly, NULL);
     PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+
+    //
+    // f
+    //
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
+    if (f->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+    } else {
+        f64 = (psVector *) f;
+    }
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+
+    //
+    // x
+    //
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+    if (x->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+    } else {
+        x64 = (psVector *) x;
+    }
+
+    //
+    // y
+    //
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    if (y->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(y64, y);
+    } else {
+        y64 = (psVector *) y;
+    }
+
+    //
+    // fErr
+    //
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
+        if (fErr->type.type != PS_TYPE_F64) {
+            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+        } else {
+            fErr64 = (psVector *) fErr;
+        }
+    }
+
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        poly = VectorFitPolynomial2DOrd(poly, mask, maskValue, f64, fErr64, x64, y64);
+        if (poly == NULL) {
+            psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
+            // Free psVectors that were created for NULL arguments.
+            if (f->type.type != PS_TYPE_F64) {
+                psFree(f64);
+            }
+
+            if (x->type.type != PS_TYPE_F64) {
+                psFree(x64);
+            }
+
+            if (y->type.type != PS_TYPE_F64) {
+                psFree(y64);
+            }
+
+            if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
+                psFree(fErr64);
+            }
+            return(NULL);
+        }
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        if (mask != NULL) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
+        }
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+        psFree(poly);
+        poly = NULL;
+    } else {
+        // Free psVectors that were created for NULL arguments.
+        if (f->type.type != PS_TYPE_F64) {
+            psFree(f64);
+        }
+
+        if (x->type.type != PS_TYPE_F64) {
+            psFree(x64);
+        }
+
+        if (y->type.type != PS_TYPE_F64) {
+            psFree(y64);
+        }
+
+        if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
+            psFree(fErr64);
+        }
+        psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type.  Returning NULL.\n");
+    }
+
+
+    // Free psVectors that were created for NULL arguments.
+    if (f->type.type != PS_TYPE_F64) {
+        psFree(f64);
+    }
+
+    if (x->type.type != PS_TYPE_F64) {
+        psFree(x64);
+    }
+
+    if (y->type.type != PS_TYPE_F64) {
+        psFree(y64);
+    }
+
+    if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
+        psFree(fErr64);
+    }
+
+    return(poly);
+}
+
+psPolynomial2D *psVectorClipFitPolynomial2D(
+    psPolynomial2D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y)
+{
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_PTR_NON_NULL(stats, NULL);
     PS_ASSERT_VECTOR_NON_NULL(f, NULL);
     PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
@@ -2648 +2216 @@
     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
     if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
         PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    }
+
+    if (mask == NULL) {
+        // XXX: Change argument order.
+        poly = RobustFit2D_nomask(poly, x, y, f, fErr);
     } else {
-        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
-        PS_VECTOR_SET_F32(fErr32, 1.0);
-    }
-
-    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
-    if (poly->type == PS_POLYNOMIAL_ORD) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D polynomial vector fitting has not been implemented.  Returning NULL.\n");
-    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
-    } else {
-        printf("XXX: ERROR: incorrect polynomial type.\n");
-    }
-    if (fErr == NULL) {
-        psFree(fErr32);
-    }
-    return(NULL);
-}
-
-psPolynomial3D *psVectorFitPolynomial3D_NEW(
-    psPolynomial1D *poly,
+        // XXX: Use maskValue.
+        // XXX: Change argument order.
+        poly = RobustFit2D(poly, mask, x, y, f, fErr);
+    }
+
+    if (poly == NULL) {
+        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
+        return(NULL);
+    }
+    return(poly);
+}
+
+/******************************************************************************
+ ******************************************************************************
+ 3-D Vector Fitting Code.
+ ******************************************************************************
+ *****************************************************************************/
+
+/******************************************************************************
+VectorFitPolynomial3DOrd(myPoly, *mask, maskValue, *f, *fErr, *x, *y, *z):
+This is a private routine which will fit a 3-D polynomial to a set of (x,
+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(
+    psPolynomial3D* myPoly,
+    const psVector* mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z)
+{
+    // These ASSERTS are redundant.
+    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL);
+
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+
+    psError(PS_ERR_UNKNOWN, true, "3-D Polynomial Fitting is not Implemented.\n");
+    return (NULL);
+}
+
+/******************************************************************************
+psVectorFitPolynomial3D():  This routine fits a 3D polynomial of arbitrary
+degree (specified in poly) to the data points (x, y, z)-(f) and returns that
+polynomial.  Types F32 and F64 are supported, however, type F32 is done via
+vector conversion only.
+ *****************************************************************************/
+psPolynomial3D *psVectorFitPolynomial3D(
+    psPolynomial3D *poly,
     const psVector *mask,
     psMaskType maskValue,
@@ -2679 +2307 @@
     const psVector *z)
 {
-    // Internal pointers for possibly NULL vectors.
-    psVector *fErr32 = NULL;
+    // Internal pointers for possibly NULL or mis-typed vectors.
+    psVector *x64 = NULL;
+    psVector *y64 = NULL;
+    psVector *z64 = NULL;
+    psVector *f64 = NULL;
+    psVector *fErr64 = NULL;
 
     PS_ASSERT_POLY_NON_NULL(poly, NULL);
     PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+
+    //
+    // f
+    //
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
+    if (f->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+    } else {
+        f64 = (psVector *) f;
+    }
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+
+    //
+    // x
+    //
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+    if (x->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+    } else {
+        x64 = (psVector *) x;
+    }
+
+    //
+    // y
+    //
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    if (y->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(y64, y);
+    } else {
+        y64 = (psVector *) y;
+    }
+
+    //
+    // z
+    //
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    if (z->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(y64, z);
+    } else {
+        z64 = (psVector *) z;
+    }
+
+    //
+    // fErr
+    //
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
+        if (fErr->type.type != PS_TYPE_F64) {
+            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+        } else {
+            fErr64 = (psVector *) fErr;
+        }
+    }
+
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        poly = VectorFitPolynomial3DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64);
+        if (poly == NULL) {
+            psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
+            // Free psVectors that were created for NULL arguments.
+            if (f->type.type != PS_TYPE_F64) {
+                psFree(f64);
+            }
+
+            if (x->type.type != PS_TYPE_F64) {
+                psFree(x64);
+            }
+
+            if (y->type.type != PS_TYPE_F64) {
+                psFree(y64);
+            }
+
+            if (z->type.type != PS_TYPE_F64) {
+                psFree(z64);
+            }
+
+            if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
+                psFree(fErr64);
+            }
+            return(NULL);
+        }
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        if (mask != NULL) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
+        }
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+        psFree(poly);
+        poly = NULL;
+    } else {
+        // Free psVectors that were created for NULL arguments.
+        if (f->type.type != PS_TYPE_F64) {
+            psFree(f64);
+        }
+
+        if (x->type.type != PS_TYPE_F64) {
+            psFree(x64);
+        }
+
+        if (y->type.type != PS_TYPE_F64) {
+            psFree(y64);
+        }
+
+        if (z->type.type != PS_TYPE_F64) {
+            psFree(z64);
+        }
+
+        if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
+            psFree(fErr64);
+        }
+        psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type.  Returning NULL.\n");
+    }
+
+
+    // Free psVectors that were created for NULL arguments.
+    if (f->type.type != PS_TYPE_F64) {
+        psFree(f64);
+    }
+
+    if (x->type.type != PS_TYPE_F64) {
+        psFree(x64);
+    }
+
+    if (y->type.type != PS_TYPE_F64) {
+        psFree(y64);
+    }
+
+    if (z->type.type != PS_TYPE_F64) {
+        psFree(z64);
+    }
+
+    if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
+        psFree(fErr64);
+    }
+
+    return(poly);
+}
+
+psPolynomial3D *psVectorClipFitPolynomial3D(
+    psPolynomial3D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z)
+{
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_PTR_NON_NULL(stats, NULL);
     PS_ASSERT_VECTOR_NON_NULL(f, NULL);
     PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
@@ -2696 +2486 @@
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
     if (fErr != NULL) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
         PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
-    } else {
-        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
-        PS_VECTOR_SET_F32(fErr32, 1.0);
     }
 
     psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
-    if (poly->type == PS_POLYNOMIAL_ORD) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D polynomial vector fitting has not been implemented.  Returning NULL.\n");
-    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
-    } else {
-        printf("XXX: ERROR: incorrect polynomial type.\n");
-    }
-    if (fErr == NULL) {
-        psFree(fErr32);
-    }
     return(NULL);
 }
 
-psPolynomial4D *psVectorFitPolynomial4D_NEW(
-    psPolynomial1D *poly,
+
+/******************************************************************************
+ ******************************************************************************
+ 4-D Vector Fitting Code.
+ ******************************************************************************
+ *****************************************************************************/
+/******************************************************************************
+VectorFitPolynomial4DOrd(myPoly, *mask, maskValue, *f, *fErr, *x, *y, *z, *t):
+This is a private routine which will fit a 4-D polynomial to a set of (x,
+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(
+    psPolynomial4D* myPoly,
+    const psVector* mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z,
+    const psVector *t)
+{
+    // These ASSERTS are redundant.
+    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(myPoly->nT, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
+    PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+
+    psError(PS_ERR_UNKNOWN, true, "4-D Polynomial Fitting is not Implemented.\n");
+    return (NULL);
+}
+
+/******************************************************************************
+psVectorFitPolynomial4D():  This routine fits a 4D polynomial of arbitrary
+degree (specified in poly) to the data points (x, y, z, t)-(f) and returns
+that polynomial.  Types F32 and F64 are supported, however, type F32 is done
+via vector conversion only.
+ *****************************************************************************/
+psPolynomial4D *psVectorFitPolynomial4D(
+    psPolynomial4D *poly,
     const psVector *mask,
     psMaskType maskValue,
@@ -2732 +2574 @@
     const psVector *t)
 {
-    // Internal pointers for possibly NULL vectors.
-    psVector *fErr32 = NULL;
+    // Internal pointers for possibly NULL or mis-typed vectors.
+    psVector *x64 = NULL;
+    psVector *y64 = NULL;
+    psVector *z64 = NULL;
+    psVector *t64 = NULL;
+    psVector *f64 = NULL;
+    psVector *fErr64 = NULL;
 
     PS_ASSERT_POLY_NON_NULL(poly, NULL);
     PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+
+    //
+    // f
+    //
     PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
+    if (f->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(f64, f);
+    } else {
+        f64 = (psVector *) f;
+    }
     if (mask != NULL) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
         PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
     }
+
+    //
+    // x
+    //
     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+    if (x->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(x64, x);
+    } else {
+        x64 = (psVector *) x;
+    }
+
+    //
+    // y
+    //
     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    if (y->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(y64, y);
+    } else {
+        y64 = (psVector *) y;
+    }
+
+    //
+    // z
+    //
     PS_ASSERT_VECTOR_NON_NULL(z, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    if (z->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(y64, z);
+    } else {
+        z64 = (psVector *) z;
+    }
+
+    //
+    // t
+    //
     PS_ASSERT_VECTOR_NON_NULL(t, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
-    PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, NULL);
-    if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    if (t->type.type != PS_TYPE_F64) {
+        PS_VECTOR_GEN_F64_FROM_F32(y64, t);
     } else {
-        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
-        PS_VECTOR_SET_F32(fErr32, 1.0);
-    }
-
-    if (poly->type == PS_POLYNOMIAL_ORD) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D polynomial vector fitting has not been implemented.  Returning NULL.\n");
-    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
-    } else {
-        printf("XXX: ERROR: incorrect polynomial type.\n");
-    }
-    if (fErr == NULL) {
-        psFree(fErr32);
-    }
-    return(NULL);
-}
-
-
-psPolynomial1D *psVectorClipFitPolynomial1D_NEW(
-    psPolynomial1D *poly,
-    psStats *stats,
-    const psVector *mask,
-    psMaskType maskValue,
-    const psVector *f,
-    const psVector *fErr,
-    const psVector *x)
-{
-    // Internal pointers for possibly NULL vectors.
-    psVector *x32 = NULL;
-    psVector *fErr32 = NULL;
-
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    if (x != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
-    }
+        t64 = (psVector *) t;
+    }
+
+    //
+    // fErr
+    //
     if (fErr != NULL) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
+        if (fErr->type.type != PS_TYPE_F64) {
+            PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr);
+        } else {
+            fErr64 = (psVector *) fErr;
+        }
+    }
+
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        poly = VectorFitPolynomial4DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64, t64);
+        if (poly == NULL) {
+            psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial.  Returning NULL.\n");
+            // Free psVectors that were created for NULL arguments.
+            if (f->type.type != PS_TYPE_F64) {
+                psFree(f64);
+            }
+
+            if (x->type.type != PS_TYPE_F64) {
+                psFree(x64);
+            }
+
+            if (y->type.type != PS_TYPE_F64) {
+                psFree(y64);
+            }
+
+            if (z->type.type != PS_TYPE_F64) {
+                psFree(z64);
+            }
+
+            if (t->type.type != PS_TYPE_F64) {
+                psFree(t64);
+            }
+
+            if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
+                psFree(fErr64);
+            }
+            return(NULL);
+        }
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        if (mask != NULL) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n");
+        }
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+        psFree(poly);
+        poly = NULL;
     } else {
-        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
-        PS_VECTOR_SET_F32(fErr32, 1.0);
-    }
-
-    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+        // Free psVectors that were created for NULL arguments.
+        if (f->type.type != PS_TYPE_F64) {
+            psFree(f64);
+        }
+
+        if (x->type.type != PS_TYPE_F64) {
+            psFree(x64);
+        }
+
+        if (y->type.type != PS_TYPE_F64) {
+            psFree(y64);
+        }
+
+        if (z->type.type != PS_TYPE_F64) {
+            psFree(z64);
+        }
+
+        if (t->type.type != PS_TYPE_F64) {
+            psFree(t64);
+        }
+
+        if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
+            psFree(fErr64);
+        }
+        psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type.  Returning NULL.\n");
+    }
+
+
     // Free psVectors that were created for NULL arguments.
-    if (x == NULL) {
-        psFree(x32);
-    }
-    if (fErr == NULL) {
-        psFree(fErr32);
-    }
-    return(NULL);
-}
-
-psPolynomial2D *psVectorClipFitPolynomial2D_NEW(
-    psPolynomial1D *poly,
-    psStats *stats,
-    const psVector *mask,
-    psMaskType maskValue,
-    const psVector *f,
-    const psVector *fErr,
-    const psVector *x,
-    const psVector *y)
-{
-    // Internal pointers for possibly NULL vectors.
-    psVector *fErr32 = NULL;
-
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
-    if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
-    } else {
-        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
-        PS_VECTOR_SET_F32(fErr32, 1.0);
-    }
-
-    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
-    if (fErr == NULL) {
-        psFree(fErr32);
-    }
-    return(NULL);
-}
-
-psPolynomial3D *psVectorClipFitPolynomial3D_NEW(
-    psPolynomial1D *poly,
-    psStats *stats,
-    const psVector *mask,
-    psMaskType maskValue,
-    const psVector *f,
-    const psVector *fErr,
-    const psVector *x,
-    const psVector *y,
-    const psVector *z)
-{
-    // Internal pointers for possibly NULL vectors.
-    psVector *fErr32 = NULL;
-
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
-    if (fErr != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
-    } else {
-        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
-        PS_VECTOR_SET_F32(fErr32, 1.0);
-    }
-
-    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
-    if (fErr == NULL) {
-        psFree(fErr32);
-    }
-    return(NULL);
-}
-
-psPolynomial4D *psVectorClipFitPolynomial4D_NEW(
-    psPolynomial1D *poly,
+    if (f->type.type != PS_TYPE_F64) {
+        psFree(f64);
+    }
+
+    if (x->type.type != PS_TYPE_F64) {
+        psFree(x64);
+    }
+
+    if (y->type.type != PS_TYPE_F64) {
+        psFree(y64);
+    }
+
+    if (z->type.type != PS_TYPE_F64) {
+        psFree(z64);
+    }
+
+    if (t->type.type != PS_TYPE_F64) {
+        psFree(t64);
+    }
+
+    if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) {
+        psFree(fErr64);
+    }
+
+    return(poly);
+}
+
+
+psPolynomial4D *psVectorClipFitPolynomial4D(
+    psPolynomial4D *poly,
     psStats *stats,
     const psVector *mask,
@@ -2921 +2764 @@
     const psVector *t)
 {
-    // Internal pointers for possibly NULL vectors.
-    psVector *fErr32 = NULL;
-
     PS_ASSERT_POLY_NON_NULL(poly, NULL);
     PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_PTR_NON_NULL(stats, NULL);
     PS_ASSERT_VECTOR_NON_NULL(f, NULL);
     PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
@@ -2938 +2779 @@
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
     PS_ASSERT_VECTOR_NON_NULL(t, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
@@ -2947 +2788 @@
         PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
         PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
-    } else {
-        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
-        PS_VECTOR_SET_F32(fErr32, 1.0);
     }
 
     psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
-    if (fErr == NULL) {
-        psFree(fErr32);
-    }
+
     return(NULL);
 }

trunk/psLib/src/math/psMinimize.h

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-07 23:46:04 $
+ *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-11 22:18:40 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -85 +85 @@
  *  @return psPolynomial1D*    polynomial fit
  */
-psPolynomial1D* psVectorFitPolynomial1D(
-    psPolynomial1D* poly,            ///< Polynomial to fit
-    const psVector* x,                 ///< Ordinates (or NULL to just use the indices)
-    const psVector* y,                 ///< Coordinates
-    const psVector* yErr               ///< Errors in coordinates, or NULL
-);
-
-
-
-
-
-psPolynomial1D *psVectorFitPolynomial1D_NEW(
+
+psPolynomial1D *psVectorFitPolynomial1D(
     psPolynomial1D *poly,
     const psVector *mask,
@@ -105 +95 @@
 );
 
-psPolynomial2D *psVectorFitPolynomial2D_NEW(
-    psPolynomial1D *poly,
+psPolynomial2D *psVectorFitPolynomial2D(
+    psPolynomial2D *poly,
     const psVector *mask,
     psMaskType maskValue,
@@ -115 +105 @@
 );
 
-psPolynomial3D *psVectorFitPolynomial3D_NEW(
-    psPolynomial1D *poly,
+psPolynomial3D *psVectorFitPolynomial3D(
+    psPolynomial3D *poly,
     const psVector *mask,
     psMaskType maskValue,
@@ -126 +116 @@
 );
 
-psPolynomial4D *psVectorFitPolynomial4D_NEW(
-    psPolynomial1D *poly,
+psPolynomial4D *psVectorFitPolynomial4D(
+    psPolynomial4D *poly,
     const psVector *mask,
     psMaskType maskValue,
@@ -139 +129 @@
 
 
-psPolynomial1D *psVectorClipFitPolynomial1D_NEW(
+psPolynomial1D *psVectorClipFitPolynomial1D(
     psPolynomial1D *poly,
     psStats *stats,
@@ -149 +139 @@
 );
 
-psPolynomial2D *psVectorClipFitPolynomial2D_NEW(
-    psPolynomial1D *poly,
+psPolynomial2D *psVectorClipFitPolynomial2D(
+    psPolynomial2D *poly,
     psStats *stats,
     const psVector *mask,
@@ -160 +150 @@
 );
 
-psPolynomial3D *psVectorClipFitPolynomial3D_NEW(
-    psPolynomial1D *poly,
+psPolynomial3D *psVectorClipFitPolynomial3D(
+    psPolynomial3D *poly,
     psStats *stats,
     const psVector *mask,
@@ -172 +162 @@
 );
 
-psPolynomial4D *psVectorClipFitPolynomial4D_NEW(
-    psPolynomial1D *poly,
+psPolynomial4D *psVectorClipFitPolynomial4D(
+    psPolynomial4D *poly,
     psStats *stats,
     const psVector *mask,
@@ -183 +173 @@
     const psVector *z,
     const psVector *t
-);
-
-
-
-
-/** Derive a one-dimensional spline fit.
- *
- *  Given a psSpline1D data structure and a set of x,y vectors, this routine
- *  generates the linear splines which satisfy those data points.
- *
- *  @return psSpline1D*:  the calculated one-dimensional splines
- */
-psSpline1D *psVectorFitSpline1D(
-    psSpline1D *mySpline,              ///< The spline which will be generated.
-    const psVector* x,                 ///< Ordinates (or NULL to just use the indices)
-    const psVector* y,                 ///< Coordinates
-    const psVector* yErr               ///< Errors in coordinates, or NULL
-);
-psSpline1D *psVectorFitSpline1DNEW(
-    const psVector* x,                 ///< Ordinates (or NULL to just use the indices)
-    const psVector* y,                 ///< Coordinates
-    int nKnots
 );
 
@@ -233 +201 @@
     const psVector *yErr,              ///< Errors in the measurement coordinates
     psMinimizeLMChi2Func func          ///< Specified function
+);
+
+bool psMinimizeGaussNewtonDelta (
+    psVector *delta,
+    const psVector *params,
+    const psVector *paramMask,
+    const psArray  *x,
+    const psVector *y,
+    const psVector *yErr,
+    psMinimizeLMChi2Func func
 );
 

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.119 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-08 00:14:31 $
+*  @version $Revision: 1.120 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-11 22:18:40 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -564 +564 @@
  of numbers with the specified mean and sigma.
  
-XXX: It's possible to have a different seed everutime.  However, for now,
+XXX: It's possible to have a different seed everytime.  However, for now,
 for testability, we use a common seed.
  *****************************************************************************/

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
 *  splines.
 *
-*  @version $Revision: 1.122 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-08 00:14:32 $
+*  @version $Revision: 1.123 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-11 22:18:40 $
 *
 *
@@ -207 +207 @@
 }
 
+
+/*****************************************************************************
+CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a
+tabulated function at n points, this routine calculates the second
+derivatives of the interpolating cubic splines at those n points.
+ 
+The first and second derivatives at the endpoints, undefined in the SDR, are
+here defined to be 0.0.  They can be modified via ypo and yp1.
+ 
+This routine assumes that vectors x and y are of the appropriate types/sizes
+(F32).
+ 
+XXX: This algorithm is derived from the Numerical Recipes.
+XXX: use recycled vectors for internal data.
+XXX: do an F64 version?
+ *****************************************************************************/
+static psF32 *calculateSecondDerivs(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
+                                    const psVector* y)        ///< Coordinates
+{
+    psTrace(".psLib.dataManip.calculateSecondDerivs", 4,
+            "---- calculateSecondDerivs() begin ----\n");
+
+    psS32 i;
+    psS32 k;
+    psF32 sig;
+    psF32 p;
+    psS32 n = y->n;
+    psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32));
+    psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32));
+    psF32 *X = (psF32 *) & (x->data.F32[0]);
+    psF32 *Y = (psF32 *) & (y->data.F32[0]);
+    psF32 qn;
+
+    // XXX: The second derivatives at the endpoints, undefined in the SDR,
+    // are set in psConstants.h: PS_LEFT_SPLINE_DERIV, PS_RIGHT_SPLINE_DERIV.
+    derivs2[0] = -0.5;
+    u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - PS_LEFT_SPLINE_DERIV);
+
+    for (i=1;i<=(n-2);i++) {
+        sig = (X[i] - X[i-1]) / (X[i+1] - X[i-1]);
+        p = sig * derivs2[i-1] + 2.0;
+        derivs2[i] = (sig - 1.0) / p;
+        u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1]));
+        u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p;
+
+        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
+                "X[%d] is %f\n", i, X[i]);
+        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
+                "Y[%d] is %f\n", i, Y[i]);
+        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
+                "u[%d] is %f\n", i, u[i]);
+    }
+
+    qn = 0.5;
+    u[n-1] = (3.0/(X[n-1]-X[n-2])) * (PS_RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2]));
+    derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0);
+
+    for (k=(n-2);k>=0;k--) {
+        derivs2[k] = derivs2[k] * derivs2[k+1] + u[k];
+
+        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
+                "derivs2[%d] is %f\n", k, derivs2[k]);
+    }
+
+    psFree(u);
+    psTrace(".psLib.dataManip.calculateSecondDerivs", 4,
+            "---- calculateSecondDerivs() end ----\n");
+    return(derivs2);
+}
+
+/*****************************************************************************/
+/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
+/*****************************************************************************/
+
+/*****************************************************************************
+psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y
+vectors, this routine generates the linear or cublic splines which satisfy
+those data points.
+ 
+The formula for calculating the spline polynomials is derived from Numerical
+Recipes in C.  The basic idea is that the polynomial is
+ (1)     y = (A * y[0]) +
+ (2)         (B * y[1]) +
+ (3)         ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 +
+ (4)         ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0
+Where:
+ H = x[1]-x[0]
+ A = (x[1]-x)/H
+ B = (x-x[0])/H
+The bulk of the code in this routine is the expansion of the above equation
+into a polynomial in terms of x, and then saving the coefficients of the
+powers of x in the spline polynomials.  This gets pretty complicated.
+ 
+XXX: usage of yErr is not specified in IfA documentation.
+ 
+XXX: Is the x argument redundant?  What do we do if the x argument is
+supplied, but does not equal the knots specified in mySpline?
+ 
+XXX: can psSpline be NULL?
+ 
+XXX: reimplement this assuming that mySpline is NULL?
+ 
+XXX: What happens if X is NULL, then an index vector is generated for X, but
+that index vector lies outside the range vectors in mySpline?
+ 
+XXX: Assumes mySpline->knots is psF32.  Must add psU32 and psF64.
+ *****************************************************************************/
+psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,     ///< The spline which will be generated.
+                                const psVector* x,        ///< Ordinates (or NULL to just use the indices)
+                                const psVector* y,        ///< Coordinates
+                                const psVector* yErr)     ///< Errors in coordinates, or NULL
+{
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
+    if (mySpline != NULL) {
+        PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL);
+    }
+    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
+            "---- psVectorFitSpline1D() begin ----\n");
+    psS32 numSplines = (y->n)-1;
+    psF32 tmp;
+    psF32 H;
+    psS32 i;
+    psF32 slope;
+    psVector *x32 = NULL;
+    psVector *y32 = NULL;
+    psVector *yErr32 = NULL;
+    static psVector *x32Static = NULL;
+    static psVector *y32Static = NULL;
+    static psVector *yErr32Static = NULL;
+
+    PS_VECTOR_CONVERT_F64_TO_F32_STATIC(y, y32, y32Static);
+
+    // If yErr==NULL, set all errors equal.
+    if (yErr == NULL) {
+        PS_VECTOR_GEN_YERR_STATIC_F32(yErr32Static, y->n);
+        yErr32 = yErr32Static;
+    } else {
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(yErr, NULL);
+        PS_VECTOR_CONVERT_F64_TO_F32_STATIC(yErr, yErr32, yErr32Static);
+    }
+
+    // If x==NULL, create an x32 vector with x values set to (0:n).
+    if (x == NULL) {
+        PS_VECTOR_GEN_X_INDEX_STATIC_F32(x32Static, y->n);
+        x32 = x32Static;
+    } else {
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
+        PS_VECTOR_CONVERT_F64_TO_F32_STATIC(x, x32, x32Static);
+    }
+    PS_ASSERT_VECTORS_SIZE_EQUAL(x32, y32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(yErr32, y32, NULL);
+
+    /*
+        XXX:
+        This can not be implemented until SDR states what order spline should be
+        created.
+        Should we error if mySpline is not NULL?
+        Should we error if mySPline is not NULL?
+    */
+    if (mySpline == NULL) {
+        mySpline = psSpline1DAllocGeneric(x32, 3);
+    }
+    PS_ASSERT_PTR_NON_NULL(mySpline, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(mySpline->n, NULL);
+
+    if (y32->n != (1 + mySpline->n)) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                "data size / spline size mismatch (%d %d)\n",
+                y32->n, mySpline->n);
+        return(NULL);
+    }
+
+    // If these are linear splines, which means their polynomials will have
+    // two coefficients, then we do the simple calculation.
+    if (2 == (mySpline->spline[0])->n) {
+        for (i=0;i<mySpline->n;i++) {
+            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
+                    (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]);
+            (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
+                                              (slope * mySpline->knots->data.F32[i]);
+
+            (mySpline->spline[i])->coeff[1] = slope;
+            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
+                    "---- mySpline %d coeffs are (%f, %f)\n", i,
+                    (mySpline->spline[i])->coeff[0],
+                    (mySpline->spline[i])->coeff[1]);
+        }
+        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
+                "---- Exiting psVectorFitSpline1D()()\n");
+        return((psSpline1D *) mySpline);
+    }
+
+    // Check if these are cubic splines (n==4).  If not, psError.
+    if (4 != (mySpline->spline[0])->n) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                "Don't know how to generate %d-order splines.",
+                (mySpline->spline[0])->n-1);
+        return(NULL);
+    }
+
+    // If we get here, then we know these are cubic splines.  We first
+    // generate the second derivatives at each data point.
+    mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
+    for (i=0;i<y32->n;i++)
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
+
+    // We generate the coefficients of the spline polynomials.  I can't
+    // concisely explain how this code works.  See above function comments
+    // and Numerical Recipes in C.
+    for (i=0;i<numSplines;i++) {
+        H = x32->data.F32[i+1] - x32->data.F32[i];
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
+                "x data (%f - %f) (%f)\n",
+                x32->data.F32[i],
+                x32->data.F32[i+1], H);
+        //
+        // ******** Calculate 0-order term ********
+        //
+        // From (1)
+        (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
+        // From (2)
+        ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
+        // From (3)
+        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
+        tmp-= (x32->data.F32[i+1] / H);
+        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[0])+= tmp;
+        // From (4)
+        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
+        tmp+= (x32->data.F32[i] / H);
+        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[0])+= tmp;
+
+        //
+        // ******** Calculate 1-order term ********
+        //
+        // From (1)
+        (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
+        // From (2)
+        ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
+        // From (3)
+        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
+        tmp+= (1.0 / H);
+        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[1])+= tmp;
+        // From (4)
+        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
+        tmp-= (1.0 / H);
+        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[1])+= tmp;
+
+        //
+        // ******** Calculate 2-order term ********
+        //
+        // From (3)
+        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
+        // From (4)
+        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
+
+        //
+        // ******** Calculate 3-order term ********
+        //
+        // From (3)
+        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
+        // From (4)
+        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
+
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
+
+    }
+
+    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
+            "---- psVectorFitSpline1D() end ----\n");
+    return(mySpline);
+}
+
+
+
+
+
+
+
+
+
+/*****************************************************************************
+psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y
+ 
+xF32 and yF32 are internal psVectors which are used to hold the psF32 versions
+of the input data, if necessary.  xPtr and yPtr are pointers to either xF32 or
+the x argument.  All computation is done on xPtr and yPtr.  xF32 and yF32 will
+simply be psFree() at the end.
+ 
+XXX: nKnots makes no sense.  This number is always equal to the size of the x
+an y vectors.
+ 
+XXX: How do we specify the spline order?  For now, order=3.
+ *****************************************************************************/
+#define PS_XXX_SPLINE_ORDER 3
+psSpline1D *psVectorFitSpline1DNEW(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
+                                   const psVector* y,        ///< Coordinates
+                                   int nKnots)
+{
+    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n");
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
+    //    PS_ASSERT_INT_EQUAL(y->n, nKnots);
+
+    //
+    // The following code ensures that xPtr points to a psF32 version of the
+    // ordinate data.
+    //
+    psVector *xF32 = NULL;
+    psVector *xPtr = NULL;
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
+        if (PS_TYPE_F64 == x->type.type) {
+            xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
+            for (psS32 i = 0 ; i < x->n ; i++) {
+                xF32->data.F32[i] = (psF32) x->data.F64[i];
+            }
+            xPtr = xF32;
+        } else if (PS_TYPE_F32 == x->type.type) {
+            xPtr = (psVector *) x;
+        } else {
+            psError(PS_ERR_UNKNOWN, true, "psVector x is wrong type.\n");
+            return(NULL);
+        }
+    } else {
+        // If x==NULL, create an x32 vector with x values set to (0:n).
+        xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
+        for (psS32 i = 0 ; i < x->n ; i++) {
+            xF32->data.F32[i] = (psF32) i;
+        }
+        xPtr = xF32;
+    }
+
+    //
+    // If y is of type psF64, then create a new vector yF32 and convert the
+    // y elements.  Regardless of y's type, we create a yPtr which will be
+    // used in the remainder of this function.
+    //
+    psVector *yF32 = NULL;
+    psVector *yPtr = NULL;
+    if (PS_TYPE_F64 == y->type.type) {
+        yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
+        for (psS32 i = 0 ; i < y->n ; i++) {
+            yF32->data.F32[i] = (psF32) y->data.F64[i];
+        }
+        yPtr = yF32;
+    } else {
+        yPtr = (psVector *) y;
+    }
+
+    psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER);
+
+    psS32 numSplines = nKnots - 1;
+    psF32 tmp;
+    psF32 H;
+    psS32 i;
+    psF32 slope;
+    // XXX: get rid of x32 and y32 (this is from old code)
+    psVector *x32 = xPtr;
+    psVector *y32 = yPtr;
+
+    // If these are linear splines, which means their polynomials will have
+    // two coefficients, then we do the simple calculation.
+    if (1 == PS_XXX_SPLINE_ORDER) {
+        for (i=0;i<mySpline->n;i++) {
+            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
+                    (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]);
+            (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
+                                              (slope * mySpline->knots->data.F32[i]);
+
+            (mySpline->spline[i])->coeff[1] = slope;
+            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
+                    "---- mySpline %d coeffs are (%f, %f)\n", i,
+                    (mySpline->spline[i])->coeff[0],
+                    (mySpline->spline[i])->coeff[1]);
+        }
+        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
+                "---- Exiting psVectorFitSpline1D()()\n");
+        return((psSpline1D *) mySpline);
+    }
+
+    //
+    // Check if these are cubic splines (n==4).  If not, psError.
+    //
+    if (3 != PS_XXX_SPLINE_ORDER) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER);
+        return(NULL);
+    }
+
+    //
+    // If we get here, then we know these are cubic splines.  We first
+    // generate the second derivatives at each data point.
+    //
+    mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
+    for (i=0;i<y32->n;i++)
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
+
+    //
+    // We generate the coefficients of the spline polynomials.  I can't
+    // concisely explain how this code works.  See above function comments
+    // and Numerical Recipes in C.
+    //
+    for (i=0;i<numSplines;i++) {
+        H = x32->data.F32[i+1] - x32->data.F32[i];
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
+                "x data (%f - %f) (%f)\n",
+                x32->data.F32[i],
+                x32->data.F32[i+1], H);
+        //
+        // ******** Calculate 0-order term ********
+        //
+        // From (1)
+        (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
+        // From (2)
+        ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
+        // From (3)
+        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
+        tmp-= (x32->data.F32[i+1] / H);
+        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[0])+= tmp;
+        // From (4)
+        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
+        tmp+= (x32->data.F32[i] / H);
+        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[0])+= tmp;
+
+        //
+        // ******** Calculate 1-order term ********
+        //
+        // From (1)
+        (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
+        // From (2)
+        ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
+        // From (3)
+        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
+        tmp+= (1.0 / H);
+        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[1])+= tmp;
+        // From (4)
+        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
+        tmp-= (1.0 / H);
+        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[1])+= tmp;
+
+        //
+        // ******** Calculate 2-order term ********
+        //
+        // From (3)
+        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
+        // From (4)
+        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
+
+        //
+        // ******** Calculate 3-order term ********
+        //
+        // From (3)
+        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
+        // From (4)
+        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
+
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
+
+    }
+
+    psFree(xF32);
+    psFree(yF32);
+
+    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
+            "---- psVectorFitSpline1D() end ----\n");
+
+    return(mySpline);
+}
+
+
+
+
+
 /*****************************************************************************/
 /*  FUNCTION IMPLEMENTATION - PUBLIC                                         */
@@ -283 +781 @@
         }
     } else {
-        printf("XXX: Generate an error here.\n");
+        psError(PS_ERR_UNKNOWN, true, "psVector in has an incorrect type.\n");
         return(NULL);
     }

trunk/psLib/src/math/psSpline.h

@@ -10 +10 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-08 00:02:48 $
+ *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-11 22:18:40 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -118 +118 @@
 );
 
+/** Derive a one-dimensional spline fit.
+ *
+ *  Given a psSpline1D data structure and a set of x,y vectors, this routine
+ *  generates the linear splines which satisfy those data points.
+ *
+ *  @return psSpline1D*:  the calculated one-dimensional splines
+ */
+psSpline1D *psVectorFitSpline1D(
+    psSpline1D *mySpline,              ///< The spline which will be generated.
+    const psVector* x,                 ///< Ordinates (or NULL to just use the indices)
+    const psVector* y,                 ///< Coordinates
+    const psVector* yErr               ///< Errors in coordinates, or NULL
+);
+
+psSpline1D *psVectorFitSpline1DNEW(
+    const psVector* x,                 ///< Ordinates (or NULL to just use the indices)
+    const psVector* y,                 ///< Coordinates
+    int nKnots
+);
+
 /** \} */ // End of MathGroup Functions
 

trunk/psLib/src/math/psStats.c

@@ -17 +17 @@
  *
  *
- *  @version $Revision: 1.144 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-07 21:40:09 $
+ *  @version $Revision: 1.145 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-11 22:18:40 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -746 +746 @@
 XXX: Write a general routine which smoothes a psVector.  This routine should
 call that.  Is that possible?
+ 
+XXX: This is, or will be, prettier than the previous
+psVectorSmoothHistGaussian().  However, it is not being used, yet.
  *****************************************************************************/
 psVector *p_psVectorSmoothHistGaussianNEW(psHistogram *histogram,
@@ -1549 +1552 @@
 
         // Determine the coefficients of the polynomial.
-        myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
+        //        myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
+        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x);
         if (myPoly == NULL) {
             psError(PS_ERR_UNEXPECTED_NULL,
@@ -1825 +1829 @@
         psPolynomial1D *tmpPoly = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
         // XXX: What about the NULL x argument?
-        tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, y, NULL);
+        tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, 0, y, NULL, NULL);
         if (tmpPoly == NULL) {
             psLogMsg(__func__, PS_LOG_WARN,
@@ -1936 +1940 @@
         }
     } else {
-        printf("XXX: Generate an error here.\n");
+        psError(PS_ERR_UNKNOWN, true, "Unallowable vector type.\n");
         return(NULL);
     }
@@ -1950 +1954 @@
 psVector *Fit1DGaussian(psVector *x, psVector*y)
 {
-    printf("XXX: Generate an error here.\n");
-    printf("XXX: Error: This function was previously part of psStats.c, was removed, was purged from CVS, and now needs to be retrieved from tape.\n");
+    psError(PS_ERR_UNKNOWN, true, "This code has not been implemented yet.\n");
+    // XXX: This function was previously part of psStats.c, was removed, was
+    // purged from CVS, and now needs to be retrieved from tape.
+
     return(NULL);
 }
@@ -2168 +2174 @@
             stats->robustLQ = binLo25F32;
             stats->robustUQ = binHi25F32;
-            // XXX: No idea how to calculate stats->stdev
 
             // PS_BIN_MIDPOINT(robustHistogram, modeBinNum);
-
             // XXX: I think sumNfit == sumN50 here.
             stats->robustNfit = -1;