Index: trunk/psLib/src/math/psMinimizePowell.c
===================================================================
--- trunk/psLib/src/math/psMinimizePowell.c	(revision 6484)
+++ trunk/psLib/src/math/psMinimizePowell.c	(revision 7764)
@@ -11,6 +11,6 @@
  *  NOTE: XXX: The SDR is silent about data types.  F32 is implemented here.
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-24 23:43:15 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-06-30 01:56:22 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -64,14 +64,4 @@
 /*****************************************************************************/
 /* TYPE DEFINITIONS                                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* GLOBAL VARIABLES                                                          */
-/*****************************************************************************/
-static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
-static psVector *PowellValue;
-static psVector *PowellError;
-/*****************************************************************************/
-/* FILE STATIC VARIABLES                                                     */
 /*****************************************************************************/
 
@@ -239,4 +229,5 @@
 psTrace(__func__, 4, "Final bracket (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc); \
 psTrace(__func__, 4, "---- p_psDetermineBracket() end ----\n"); \
+psFree(tmp); \
 return(bracket); \
 
@@ -256,5 +247,5 @@
     psF32 fc = 0.0;
     psS32 iter = 0;
-    PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32);
+    psVector *tmp = psVectorAlloc(params->n, PS_TYPE_F32);
     psBool boolLineIsNull = true;
     psF32 prevMin = 0.0;
@@ -268,4 +259,5 @@
         psTrace(__func__, 2, "p_psDetermineBracket() called with zero line vector.\n");
         psTrace(__func__, 4, "---- p_psDetermineBracket() end (NULL) ----\n");
+        psFree(tmp);
         return(NULL);
     }
@@ -374,8 +366,4 @@
     psF32 fn = 0.0;
     psF32 mul = 0.0;
-    PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32);
-    PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32);
-    PS_VECTOR_GEN_STATIC_RECYCLED(tmpc, params->n, PS_TYPE_F32);
-    PS_VECTOR_GEN_STATIC_RECYCLED(tmpn, params->n, PS_TYPE_F32);
     psS32 i = 0;
     psS32 boolLineIsNull = true;
@@ -405,4 +393,10 @@
     }
     numIterations = 0;
+
+    psVector *tmpa = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *tmpb = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *tmpc = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *tmpn = psVectorAlloc(params->n, PS_TYPE_F32);
+
     while (numIterations < PS_LINEMIN_MAX_ITERATIONS) {
         numIterations++;
@@ -462,4 +456,8 @@
             psFree(bracket);
             psTrace(__func__, 4, "---- LineMin() end.a (%f) (%f) ----\n", mul, min->value);
+            psFree(tmpa);
+            psFree(tmpb);
+            psFree(tmpc);
+            psFree(tmpn);
             return(mul);
         }
@@ -472,4 +470,8 @@
 
     psFree(bracket);
+    psFree(tmpa);
+    psFree(tmpb);
+    psFree(tmpc);
+    psFree(tmpn);
     return(mul);
 }
@@ -503,7 +505,4 @@
     PS_ASSERT_PTR_NON_NULL(func, NULL);
     psS32 numDims = params->n;
-    PS_VECTOR_GEN_STATIC_RECYCLED(pQP, numDims, PS_TYPE_F32);
-    PS_VECTOR_GEN_STATIC_RECYCLED(u, numDims, PS_TYPE_F32);
-    PS_VECTOR_GEN_STATIC_RECYCLED(Q, numDims, PS_TYPE_F32);
     psS32 i = 0;
     psS32 j = 0;
@@ -533,4 +532,9 @@
     }
     PS_ASSERT_VECTORS_SIZE_EQUAL(params, myParamMask, NULL);
+
+
+    psVector *pQP = psVectorAlloc(numDims, PS_TYPE_F32);
+    psVector *u   = psVectorAlloc(numDims, PS_TYPE_F32);
+    psVector *Q   = psVectorAlloc(numDims, PS_TYPE_F32);
 
     // 1: Set v[i] to be the unit vectors for each dimension in params
@@ -578,4 +582,7 @@
                             "Could not perform line minimization.  Returning FALSE.\n");
                     psFree(v);
+                    psFree(pQP);
+                    psFree(u);
+                    psFree(Q);
                     return(false);
                 }
@@ -620,4 +627,7 @@
                     "Could not perform line minimization.  Returning FALSE.\n");
             psFree(v);
+            psFree(pQP);
+            psFree(u);
+            psFree(Q);
             return(false);
         }
@@ -626,4 +636,7 @@
         if (dummyMin.value > currFuncVal) {
             psFree(v);
+            psFree(pQP);
+            psFree(u);
+            psFree(Q);
             min->iter = iterationNumber;
             min->value = currFuncVal;
@@ -664,4 +677,7 @@
         if (fabs(baseFuncVal - currFuncVal) <= min->tol) {
             psFree(v);
+            psFree(pQP);
+            psFree(u);
+            psFree(Q);
             // XXX: Ensure that currFuncVal is the correct value to use here.
             min->value = currFuncVal;
@@ -674,4 +690,7 @@
 
     psFree(v);
+    psFree(pQP);
+    psFree(u);
+    psFree(Q);
     min->iter = iterationNumber;
     psTrace(__func__, 4, "---- psMinimizePowell() end (0) (false) ----\n");
@@ -695,6 +714,4 @@
     PS_ASSERT_VECTOR_NON_NULL(params, NAN);
     PS_ASSERT_VECTOR_NON_EMPTY(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(PowellValue, NAN);
-    PS_ASSERT_VECTOR_NON_EMPTY(PowellValue, NAN);
     PS_ASSERT_PTR_NON_NULL(coords, NAN);
 
@@ -704,13 +721,28 @@
     psVector *tmp;
 
-    tmp = Chi2PowellFunc(params, coords);
-    if (PowellError == NULL) {
+    psVector *values = coords->data[coords->n];
+    psVector *errors = coords->data[coords->n + 1];
+    psMinimizeChi2PowellFunc *func = coords->data[coords->n + 2];
+
+    PS_ASSERT_VECTOR_NON_NULL(values, NAN);
+    PS_ASSERT_VECTOR_NON_EMPTY(values, NAN);
+    PS_ASSERT_VECTOR_TYPE(values, PS_TYPE_F32, NAN);
+    if (errors) {
+        PS_ASSERT_VECTOR_NON_NULL(errors, NAN);
+        PS_ASSERT_VECTOR_NON_EMPTY(errors, NAN);
+        PS_ASSERT_VECTOR_TYPE(errors, PS_TYPE_F32, NAN);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(values, errors, NAN);
+    }
+
+    tmp = (*func)(params, coords);
+
+    if (errors == NULL) {
         for (i=0;i<coords->n;i++) {
-            d = (tmp->data.F32[i] - PowellValue->data.F32[i]);
+            d = (tmp->data.F32[i] - values->data.F32[i]);
             chi2+= d * d;
         }
     } else {
         for (i=0;i<coords->n;i++) {
-            d = (tmp->data.F32[i] - PowellValue->data.F32[i]) / PowellError->data.F32[i];
+            d = (tmp->data.F32[i] - values->data.F32[i]) / errors->data.F32[i];
             chi2+= d * d;
         }
@@ -741,10 +773,20 @@
     psMinimizeChi2PowellFunc model)
 {
-    PowellValue = (psVector *) value;
-    PowellError = (psVector *) error;
-
-    Chi2PowellFunc = model;
-
-    return(psMinimizePowell(min, params, constrain->paramMask, coords, myPowellChi2Func));
+    // Generate extended version of coords array, so we can pass in extra data to the chi^2 function
+    psArray *newCoords = psArrayAlloc(coords->n + 3);
+    for (long i = 0; i < coords->n; i++) {
+        newCoords->data[i] = psMemIncrRefCounter(coords->data[i]);
+    }
+    newCoords->n = coords->n;           // We deceive everyone else as to the length
+    // Casting away const: I'm not going to hurt you, just want to increment your reference counter is all
+    newCoords->data[coords->n] = psMemIncrRefCounter((psVector*)value);
+    newCoords->data[coords->n + 1] = psMemIncrRefCounter((psVector*)error);
+    newCoords->data[coords->n + 2] = &model;
+
+    bool success = psMinimizePowell(min, params, constrain->paramMask, newCoords, myPowellChi2Func);
+
+    newCoords->data[coords->n - 1] = NULL; // We can't free the array with a function pointer on it
+    psFree(newCoords);
+    return success;
 }