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);
 }
Index: /trunk/psLib/src/dataManip/psMinimize.d
===================================================================
--- /trunk/psLib/src/dataManip/psMinimize.d	(revision 1189)
+++ /trunk/psLib/src/dataManip/psMinimize.d	(revision 1190)
@@ -3,3 +3,3 @@
   ../collections/psImage.h ../sysUtils/psTrace.h ../sysUtils/psError.h \
   ../sysUtils/psAbort.h psFunctions.h ../collections/psSort.h \
-  psMinimize.h
+  psMinimize.h psMatrix.h
Index: /trunk/psLib/src/dataManip/psMinimize.h
===================================================================
--- /trunk/psLib/src/dataManip/psMinimize.h	(revision 1189)
+++ /trunk/psLib/src/dataManip/psMinimize.h	(revision 1190)
@@ -38,5 +38,5 @@
 /** Derive a polynomial fit by chi^2 minimisation (analytically) */
 psPolynomial1D *
-psGetArrayPolynomial(psPolynomial1D myPoly, ///< Polynomial to fit
+psGetArrayPolynomial(psPolynomial1D *myPoly, ///< Polynomial to fit
                      const psVector *restrict x, ///< Ordinates (or NULL to just use the indices)
                      const psVector *restrict y, ///< Coordinates
Index: /trunk/psLib/src/math/psMinimize.c
===================================================================
--- /trunk/psLib/src/math/psMinimize.c	(revision 1189)
+++ /trunk/psLib/src/math/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);
 }
Index: /trunk/psLib/test/dataManip/tst_psMinimize04.c
===================================================================
--- /trunk/psLib/test/dataManip/tst_psMinimize04.c	(revision 1190)
+++ /trunk/psLib/test/dataManip/tst_psMinimize04.c	(revision 1190)
@@ -0,0 +1,79 @@
+/*****************************************************************************
+    This routine must ensure that the psGetArrayPolynomial works correctly.
+ *****************************************************************************/
+#include <stdio.h>
+#include "pslib.h"
+#include "psTest.h"
+#include "psMemory.h"
+#include "psVector.h"
+#include "psImage.h"
+#include "psFunctions.h"
+#include "psMinimize.h"
+#include <math.h>
+#define NUM_DATA 10
+#define POLY_ORDER 2
+
+double setData(double A,
+               double B,
+               double C,
+               double x)
+{
+    return(A + (B * x) + (C * x * x));
+}
+
+int main()
+{
+    psPolynomial1D *myPoly = NULL;
+    psVector *x = NULL;
+    psVector *y = NULL;
+    psVector *yErr = NULL;
+    int i = 0;
+    int currentId = psMemGetId();
+    int testStatus = true;
+    int memLeaks = 0;
+
+    myPoly = psPolynomial1DAlloc(POLY_ORDER+1);
+    x = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+    y = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+    yErr = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+
+    for (i=0;i<NUM_DATA;i++) {
+        x->data.F64[i] = (double) i;
+        y->data.F64[i] = setData(2.0, 3.0, 2.0, x->data.F64[i]);
+        yErr->data.F64[i] = 0.1;
+        printf("Original data %d: (%.1f %.1f)\n", i, x->data.F64[i], y->data.F64[i]);
+    }
+
+    psGetArrayPolynomial(myPoly, x, y, yErr);
+
+    for (i=0;i<POLY_ORDER+1;i++) {
+        printf("Polynomial coefficient %d is %0.1f\n", i, myPoly->coeff[i]);
+    }
+
+    for (i=0;i<NUM_DATA;i++) {
+        printf("Fitted data %d: (%.1f %.1f)\n", i, x->data.F64[i],
+               setData(myPoly->coeff[0], myPoly->coeff[1],
+                       myPoly->coeff[2], x->data.F64[i]));
+    }
+
+
+    psMemCheckCorruption(1);
+    printFooter(stdout,
+                "psMinimize functions",
+                "psMinimize(): no masks",
+                testStatus);
+
+    psMemCheckCorruption(1);
+    psFree(myPoly);
+    psFree(x);
+    psFree(y);
+    psFree(yErr);
+
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,NULL);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    return (!testStatus);
+}