Index: trunk/psLib/src/math/psPolynomial.c
===================================================================
--- trunk/psLib/src/math/psPolynomial.c	(revision 42824)
+++ trunk/psLib/src/math/psPolynomial.c	(revision 42825)
@@ -806,71 +806,4 @@
 }
 
-// XXX add 1D, 3D, 4D versions
-bool psPolynomial2DRecycle(psPolynomial2D *poly,
-                           psPolynomialType type,
-                           unsigned int nX,
-                           unsigned int nY)
-{
-    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
-
-    bool match = true;
-    match &= (poly->type == type);
-    match &= (poly->nX == type);
-    match &= (poly->nY == type);
-
-    if (!match) {
-        for (int i = 0; i < poly->nX + 1; i++) {
-            psFree (poly->coeff[i]);
-            psFree (poly->coeffErr[i]);
-            psFree (poly->coeffMask[i]);
-        }
-        psFree (poly->coeff);
-        psFree (poly->coeffErr);
-        psFree (poly->coeffMask);
-
-        poly->type = type;
-        poly->nX = nX;
-        poly->nY = nY;
-
-        poly->coeff = psAlloc((1 + nX) * sizeof(psF64 *));
-        poly->coeffErr = psAlloc((1 + nX) * sizeof(psF64 *));
-        poly->coeffMask = (psMaskType **)psAlloc((1 + nX) * sizeof(psMaskType *));
-        for (int i = 0; i < (1 + nX); i++) {
-            poly->coeff[i] = psAlloc((1 + nY) * sizeof(psF64));
-            poly->coeffErr[i] = psAlloc((1 + nY) * sizeof(psF64));
-            poly->coeffMask[i] = (psMaskType *)psAlloc((1 + nY) * sizeof(psMaskType));
-        }
-    }
-    for (int i = 0; i < (1 + nX); i++) {
-        for (int j = 0; j < (1 + nY); j++) {
-            poly->coeff[i][j] = 0.0;
-            poly->coeffErr[i][j] = 0.0;
-            poly->coeffMask[i][j] = PS_POLY_MASK_NONE;
-        }
-    }
-    return(true);
-}
-
-// XXX add 1D, 3D, 4D versions
-psPolynomial2D *psPolynomial2DCopy(psPolynomial2D *out,
-                                   psPolynomial2D *poly)
-{
-    if (out == NULL) {
-        out = psPolynomial2DAlloc (poly->type, poly->nX, poly->nY);
-    } else {
-        psPolynomial2DRecycle (out, poly->type, poly->nX, poly->nY);
-    }
-
-    for (int i = 0; i < (1 + poly->nX); i++) {
-        for (int j = 0; j < (1 + poly->nY); j++) {
-            out->coeff[i][j] = poly->coeff[i][j];
-            out->coeffErr[i][j] = poly->coeffErr[i][j];
-            out->coeffMask[i][j] = poly->coeffMask[i][j];
-        }
-    }
-    return(out);
-}
-
 psPolynomial3D* psPolynomial3DAlloc(
     psPolynomialType type,
@@ -1002,4 +935,261 @@
 }
 
+// XXX add 3D, 4D versions
+bool psPolynomial1DRecycle(psPolynomial1D *poly,
+                           psPolynomialType type,
+                           unsigned int nX)
+{
+    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
+
+    bool match = true;
+    match &= (poly->type == type);
+    match &= (poly->nX   == nX);
+
+    if (!match) {
+	polynomial1DFree (poly); // frees the coeffs
+
+        poly->type = type;
+        poly->nX   = nX;
+
+        poly->coeff = psAlloc((1 + nX) * sizeof(psF64));
+        poly->coeffErr = psAlloc((1 + nX) * sizeof(psF64));
+        poly->coeffMask = (psMaskType *)psAlloc((1 + nX) * sizeof(psMaskType));
+    }
+    for (int i = 0; i < (1 + nX); i++) {
+	poly->coeff[i] = 0.0;
+	poly->coeffErr[i] = 0.0;
+	poly->coeffMask[i] = PS_POLY_MASK_NONE;
+    }
+    return(true);
+}
+
+bool psPolynomial2DRecycle(psPolynomial2D *poly,
+                           psPolynomialType type,
+                           unsigned int nX,
+                           unsigned int nY)
+{
+    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
+
+    bool match = true;
+    match &= (poly->type == type);
+    match &= (poly->nX ==   nX);
+    match &= (poly->nY ==   nY);
+
+    if (!match) {
+	polynomial2DFree (poly); // frees the coeffs
+
+        poly->type = type;
+        poly->nX   = nX;
+        poly->nY   = nY;
+
+        poly->coeff = psAlloc((1 + nX) * sizeof(psF64 *));
+        poly->coeffErr = psAlloc((1 + nX) * sizeof(psF64 *));
+        poly->coeffMask = (psMaskType **)psAlloc((1 + nX) * sizeof(psMaskType *));
+        for (int i = 0; i < (1 + nX); i++) {
+            poly->coeff[i] = psAlloc((1 + nY) * sizeof(psF64));
+            poly->coeffErr[i] = psAlloc((1 + nY) * sizeof(psF64));
+            poly->coeffMask[i] = (psMaskType *)psAlloc((1 + nY) * sizeof(psMaskType));
+        }
+    }
+    for (int i = 0; i < (1 + nX); i++) {
+        for (int j = 0; j < (1 + nY); j++) {
+            poly->coeff[i][j] = 0.0;
+            poly->coeffErr[i][j] = 0.0;
+            poly->coeffMask[i][j] = PS_POLY_MASK_NONE;
+        }
+    }
+    return(true);
+}
+
+bool psPolynomial3DRecycle(psPolynomial3D *poly,
+                           psPolynomialType type,
+                           unsigned int nX,
+                           unsigned int nY,
+                           unsigned int nZ)
+{
+    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nZ, NULL);
+
+    bool match = true;
+    match &= (poly->type == type);
+    match &= (poly->nX   == nX);
+    match &= (poly->nY   == nY);
+    match &= (poly->nZ   == nZ);
+
+    if (!match) {
+	polynomial3DFree (poly); // frees the coeffs
+
+        poly->type = type;
+        poly->nX = nX;
+        poly->nY = nY;
+        poly->nZ = nZ;
+
+	poly->coeff = psAlloc((nX + 1) * sizeof(psF64 **));
+	poly->coeffErr = psAlloc((nX + 1) * sizeof(psF64 **));
+	poly->coeffMask = (psMaskType ***)psAlloc((nX + 1) * sizeof(psMaskType **));
+	for (int ix = 0; ix < (1 + nX); ix++) {
+	    poly->coeff[ix] = psAlloc((nY + 1) * sizeof(psF64 *));
+	    poly->coeffErr[ix] = psAlloc((nY + 1) * sizeof(psF64 *));
+	    poly->coeffMask[ix] = (psMaskType **)psAlloc((nY + 1) * sizeof(psMaskType *));
+	    for (int iy = 0; iy < (nY + 1); iy++) {
+		poly->coeff[ix][iy] = psAlloc((nZ + 1) * sizeof(psF64));
+		poly->coeffErr[ix][iy] = psAlloc((nZ + 1) * sizeof(psF64));
+		poly->coeffMask[ix][iy] = (psMaskType *)psAlloc((nZ + 1) * sizeof(psMaskType));
+	    }
+	}
+    }
+    for (int ix = 0; ix < (1 + nX); ix++) {
+        for (int iy = 0; iy < (1 + nY); iy++) {
+	    for (int iz = 0; iz < (1 + nZ); iz++) {
+		poly->coeff[ix][iy][iz]     = 0.0;
+		poly->coeffErr[ix][iy][iz]  = 0.0;
+		poly->coeffMask[ix][iy][iz] = PS_POLY_MASK_NONE;
+	    }
+        }
+    }
+    return(true);
+}
+
+bool psPolynomial4DRecycle(psPolynomial4D *poly,
+                           psPolynomialType type,
+                           unsigned int nX,
+                           unsigned int nY,
+                           unsigned int nZ,
+                           unsigned int nT)
+{
+    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nZ, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(nT, NULL);
+
+    bool match = true;
+    match &= (poly->type == type);
+    match &= (poly->nX   == nX);
+    match &= (poly->nY   == nY);
+    match &= (poly->nZ   == nZ);
+    match &= (poly->nT   == nT);
+
+    if (!match) {
+	polynomial4DFree (poly); // frees the coeffs
+
+        poly->type = type;
+        poly->nX = nX;
+        poly->nY = nY;
+        poly->nZ = nZ;
+        poly->nT = nT;
+
+	poly->coeff = psAlloc((nX + 1) * sizeof(psF64 ***));
+	poly->coeffErr = psAlloc((nX + 1) * sizeof(psF64 ***));
+	poly->coeffMask = (psMaskType ****)psAlloc((nX + 1) * sizeof(psMaskType ***));
+	for (int ix = 0; ix < (nX + 1); ix++) {
+	    poly->coeff[ix] = psAlloc((nY + 1) * sizeof(psF64 **));
+	    poly->coeffErr[ix] = psAlloc((nY + 1) * sizeof(psF64 **));
+	    poly->coeffMask[ix] = (psMaskType ***)psAlloc((nY + 1) * sizeof(psMaskType **));
+	    for (int iy = 0; iy < (nY + 1); iy++) {
+		poly->coeff[ix][iy] = psAlloc((nZ + 1) * sizeof(psF64 *));
+		poly->coeffErr[ix][iy] = psAlloc((nZ + 1) * sizeof(psF64 *));
+		poly->coeffMask[ix][iy] = (psMaskType **)psAlloc((nZ + 1) * sizeof(psMaskType *));
+		for (int iz = 0; iz < (nZ + 1); iz++) {
+		    poly->coeff[ix][iy][iz] = psAlloc((nT + 1) * sizeof(psF64));
+		    poly->coeffErr[ix][iy][iz] = psAlloc((nT + 1) * sizeof(psF64));
+		    poly->coeffMask[ix][iy][iz] = (psMaskType *)psAlloc((nT + 1) * sizeof(psMaskType));
+		}
+	    }
+	}
+    }
+    for (int ix = 0; ix < (1 + nX); ix++) {
+        for (int iy = 0; iy < (1 + nY); iy++) {
+	    for (int iz = 0; iz < (1 + nZ); iz++) {
+		for (int it = 0; it < (1 + nT); it++) {
+		    poly->coeff[ix][iy][iz][it]     = 0.0;
+		    poly->coeffErr[ix][iy][iz][it]  = 0.0;
+		    poly->coeffMask[ix][iy][iz][it] = PS_POLY_MASK_NONE;
+		}
+	    }
+	}
+    }
+    return(true);
+}
+
+// ######## Copy polynomials ########
+psPolynomial1D *psPolynomial1DCopy(psPolynomial1D *out,
+				   psPolynomial1D *poly)
+{
+    if (out == NULL) {
+	out = psPolynomial1DAlloc (poly->type, poly->nX);
+    } else {
+	psPolynomial1DRecycle (out, poly->type, poly->nX);
+    }
+
+    for (int i = 0; i < (1 + poly->nX); i++) {
+	out->coeff[i] = poly->coeff[i];
+	out->coeffErr[i] = poly->coeffErr[i];
+	out->coeffMask[i] = poly->coeffMask[i];
+    }
+    return(out);
+}
+psPolynomial2D *psPolynomial2DCopy(psPolynomial2D *out,
+				   psPolynomial2D *poly)
+{
+    if (out == NULL) {
+	out = psPolynomial2DAlloc (poly->type, poly->nX, poly->nY);
+    } else {
+	psPolynomial2DRecycle (out, poly->type, poly->nX, poly->nY);
+    }
+
+    for (int i = 0; i < (1 + poly->nX); i++) {
+	for (int j = 0; j < (1 + poly->nY); j++) {
+	    out->coeff[i][j] = poly->coeff[i][j];
+	    out->coeffErr[i][j] = poly->coeffErr[i][j];
+	    out->coeffMask[i][j] = poly->coeffMask[i][j];
+	}
+    }
+    return(out);
+}
+psPolynomial3D *psPolynomial3DCopy(psPolynomial3D *out,
+				   psPolynomial3D *poly)
+{
+    if (out == NULL) {
+	out = psPolynomial3DAlloc (poly->type, poly->nX, poly->nY, poly->nZ);
+    } else {
+	psPolynomial3DRecycle (out, poly->type, poly->nX, poly->nY, poly->nZ);
+    }
+
+    for (int ix = 0; ix < (1 + poly->nX); ix++) {
+	for (int iy = 0; iy < (1 + poly->nY); iy++) {
+	    for (int iz = 0; iz < (1 + poly->nZ); iz++) {
+		out->coeff[ix][iy][iz] = poly->coeff[ix][iy][iz];
+		out->coeffErr[ix][iy][iz] = poly->coeffErr[ix][iy][iz];
+		out->coeffMask[ix][iy][iz] = poly->coeffMask[ix][iy][iz];
+	    }
+	}
+    }
+    return(out);
+}
+psPolynomial4D *psPolynomial4DCopy(psPolynomial4D *out,
+				   psPolynomial4D *poly)
+{
+    if (out == NULL) {
+	out = psPolynomial4DAlloc (poly->type, poly->nX, poly->nY, poly->nZ, poly->nT);
+    } else {
+	psPolynomial4DRecycle (out, poly->type, poly->nX, poly->nY, poly->nZ, poly->nT);
+    }
+
+    for (int ix = 0; ix < (1 + poly->nX); ix++) {
+	for (int iy = 0; iy < (1 + poly->nY); iy++) {
+	    for (int iz = 0; iz < (1 + poly->nZ); iz++) {
+		for (int it = 0; it < (1 + poly->nT); it++) {
+		    out->coeff[ix][iy][iz][it] = poly->coeff[ix][iy][iz][it];
+		    out->coeffErr[ix][iy][iz][it] = poly->coeffErr[ix][iy][iz][it];
+		    out->coeffMask[ix][iy][iz][it] = poly->coeffMask[ix][iy][iz][it];
+		}
+	    }
+	}
+    }
+    return(out);
+}
+
 /* note these functions accept unscaled values and apply the scaling saved on poly */
 psF64 psPolynomial1DEval(const psPolynomial1D* poly,
@@ -1034,5 +1224,5 @@
 
     switch (x->type.type) {
-    case PS_TYPE_F64:
+      case PS_TYPE_F64:
         tmp = psVectorAlloc(x->n, PS_TYPE_F64);
         for (unsigned int i=0;i<x->n;i++) {
@@ -1040,5 +1230,5 @@
         }
         break;
-    case PS_TYPE_F32:
+      case PS_TYPE_F32:
         tmp = psVectorAlloc(x->n, PS_TYPE_F32);
         for (unsigned int i=0;i<x->n;i++) {
@@ -1046,5 +1236,5 @@
         }
         break;
-    default:
+      default:
         psError(PS_ERR_UNKNOWN, false, "invalid input data type.\n");
         return (NULL);
Index: trunk/psLib/src/math/psPolynomial.h
===================================================================
--- trunk/psLib/src/math/psPolynomial.h	(revision 42824)
+++ trunk/psLib/src/math/psPolynomial.h	(revision 42825)
@@ -161,11 +161,31 @@
 ) PS_ATTR_MALLOC;
 
+bool psPolynomial1DRecycle(psPolynomial1D *poly,
+                           psPolynomialType type,
+                           unsigned int nX);
 bool psPolynomial2DRecycle(psPolynomial2D *poly,
                            psPolynomialType type,
                            unsigned int nX,
                            unsigned int nY);
-
+bool psPolynomial3DRecycle(psPolynomial3D *poly,
+                           psPolynomialType type,
+                           unsigned int nX,
+                           unsigned int nY,
+                           unsigned int nZ);
+bool psPolynomial4DRecycle(psPolynomial4D *poly,
+                           psPolynomialType type,
+                           unsigned int nX,
+                           unsigned int nY,
+                           unsigned int nZ,
+                           unsigned int nT);
+
+psPolynomial1D *psPolynomial1DCopy(psPolynomial1D *out,
+                                   psPolynomial1D *poly);
 psPolynomial2D *psPolynomial2DCopy(psPolynomial2D *out,
                                    psPolynomial2D *poly);
+psPolynomial3D *psPolynomial3DCopy(psPolynomial3D *out,
+                                   psPolynomial3D *poly);
+psPolynomial4D *psPolynomial4DCopy(psPolynomial4D *out,
+                                   psPolynomial4D *poly);
 
 /** Evaluates a 1-D polynomial at specific coordinates.