Index: trunk/psLib/src/dataManip/psMinimize.c
===================================================================
--- trunk/psLib/src/dataManip/psMinimize.c	(revision 1189)
+++ trunk/psLib/src/dataManip/psMinimize.c	(revision 1190)
@@ -22,4 +22,5 @@
 #include "psSort.h"
 #include "psMinimize.h"
+#include "psMatrix.h"
 #include "float.h"
 #include <math.h>
@@ -739,12 +740,4 @@
 }
 
-
-
-
-
-
-
-
-
 /******************************************************************************
     This routine must fit a polynomial of degree myPoly to the data points
@@ -753,47 +746,76 @@
  *****************************************************************************/
 psPolynomial1D *
-psGetArrayPolynomial(psPolynomial1D myPoly,
+psGetArrayPolynomial(psPolynomial1D *myPoly,
                      const psVector *restrict x,
                      const psVector *restrict y,
                      const psVector *restrict yErr)
 {
-    /*
-        int polyOrder = myPoly->n;
-        float **A;
-     
-     
-        // Numerical Recipes routines are all index offset 1.
-        B = (float *) psAlloc((polyOrder+2) * sizeof(float));
-        ludIndex = (int *) psAlloc((polyOrder+2) * sizeof(int));
-     
-        A = (float **) psAlloc((polyOrder+2) * sizeof(float *));
-        for(i=0;i<(polyOrder+2);i++) {
-            A[i] = (float *) psAlloc((polyOrder+2) * sizeof(float));
-        }
-     
-        // Initialize data structures.
-        for(i=0;i<(polyOrder+2);i++) {
-            B[i] = 0.0;
-            ludIndex[i] = 0;
-            for(j=0;j<(polyOrder+2);j++) {
-                A[i][j] = 0.0;
-            }
-        }
-     
-        for(k=1;k<(polyOrder+2);k++) {
-            for (i=0;i<x->n;i++) {
-                B[k]+= y->data.F32[i] * (pow(x->data.F32[i], k));
-            }
-        }
-        for(k=1;k<(polyOrder+2);k++) {
-            for(i=1;i<(polyOrder+2);i++) {
-                for (i=0;i<x->n;i++) {
-                    A[k][j]+= (pow(y->data.F32[i], k) * pow(x->data.F32[i], j));
-                }
-            }
-        }
-     
-     
-    */
+    int 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;
+    int i = 0;
+    int j = 0;
+    int k = 0;
+
+    if (x->n != y->n) {
+        psAbort(__func__, "x and y arguments have different sizes.\n");
+    }
+    if (x->n != yErr->n) {
+        psAbort(__func__, "y and yErr arguments have different sizes.\n");
+    }
+
+    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);
+    outPerm = psVectorAlloc(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;
+        outPerm->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];
+    }
+
+    for(k=0;k<(polyOrder);k++) {
+        for (i=0;i<X->n;i++) {
+            B->data.F64[k]+= y->data.F64[i] * (pow(X->data.F64[i], (float) k));
+        }
+    }
+
+    for(k=0;k<(polyOrder);k++) {
+        for(j=0;j<(polyOrder);j++) {
+            for (i=0;i<X->n;i++) {
+                A->data.F64[k][j]+= pow(X->data.F64[i], (float) (k + j));
+            }
+        }
+    }
+
+    ALUD = psMatrixLUD(ALUD, outPerm, A);
+    coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
+
+    for(k=0;k<(polyOrder);k++) {
+        myPoly->coeff[k] = coeffs->data.F64[k];
+    }
+
+    psFree(A);
+    psFree(ALUD);
+    psFree(B);
+    psFree(coeffs);
+    psFree(X);
+    psFree(outPerm);
+
     return(NULL);
 }