Index: /trunk/psLib/src/dataManip/psFunctions.c
===================================================================
--- /trunk/psLib/src/dataManip/psFunctions.c	(revision 1830)
+++ /trunk/psLib/src/dataManip/psFunctions.c	(revision 1831)
@@ -7,6 +7,6 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-17 02:14:52 $
+ *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-18 01:50:45 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1809,10 +1809,10 @@
 p_psInterpolate1D(): This routine will take as input n-element floating
 point arrays domain and range, and the x value, assumed to lie with the
-domain vector.  It produces as output the n-order LaGrange interpolated
+domain vector.  It produces as output the (n-1)-order LaGrange interpolated
 value of x.
  
 XXX: do we error check for non-distinct domain values?
  *****************************************************************************/
-float p_psFullInterpolate1DF32(float *domain,
+float p_ps1DFullInterpolateF32(float *domain,
                                float *range,
                                int n,
@@ -1822,21 +1822,39 @@
     int m;
     static psVector *p = NULL;
-    psVectorRecycle(p, n, PS_TYPE_F32);
+    p = psVectorRecycle(p, n, PS_TYPE_F32);
     p_psMemSetPersistent(p, true);
+    p_psMemSetPersistent(p->data.F32, true);
+
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
+            "---- p_ps1DFullInterpolateF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n);
+
+    for (i=0;i<n;i++) {
+        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                "domain/range is (%f %f)\n", domain[i], range[i]);
+    }
+
+    for (i=0;i<n;i++) {
+        p->data.F32[i] = range[i];
+        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
+
+    }
 
     // From NR, during each iteration of the m loop, we are computing the
     // p_{i ... i+m} terms.
-    for (m=0;m<n;m++) {
+    for (m=1;m<n;m++) {
         for (i=0;i<n-m;i++) {
-            if (m == 0) {
-                p->data.F32[i] = range[i];
-            } else {
-                // From NR: we are computing P_{i ... i+m}
-                p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) +
-                                  ((range[i]-x) * p->data.F32[i+1])) /
-                                 (domain[i] - domain[i+m]);
-            }
-        }
-    }
+            // From NR: we are computing P_{i ... i+m}
+            p->data.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) +
+                              ((domain[i]-x) * p->data.F32[i+1])) /
+                             (domain[i] - domain[i+m]);
+            //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]);
+            psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                    "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
+        }
+    }
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
+            "---- p_ps1DFullInterpolateF32() end ----\n");
+
     return(p->data.F32[0]);
 }
@@ -1844,5 +1862,5 @@
 
 // This is a
-float p_psInterpolate1DF32(float *domain,
+float p_ps1DInterpolateF32(float *domain,
                            float *range,
                            int n,
@@ -1854,5 +1872,9 @@
     int origin;
 
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
+            "---- p_ps1DInterpolateF32() begin ----\n");
+
     binNum = VectorBinDisectF32(domain, n, x);
+
     if (0 == numIntPoints%2) {
         origin = binNum - ((numIntPoints/2) - 1);
@@ -1871,5 +1893,7 @@
     }
 
-    return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x));
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
+            "---- p_ps1DInterpolateF32() end ----\n");
+    return(p_ps1DFullInterpolateF32(&domain[origin], &range[origin], order+1, x));
 }
 
@@ -1882,9 +1906,12 @@
 XXX: This stuff does not work with a mask.
  *****************************************************************************/
-float p_psInterpolate1D(psVector *domain,
-                        psVector *range,
-                        int order,
-                        psScalar *x)
-{
+psScalar *p_psVectorInterpolate(psVector *domain,
+                                psVector *range,
+                                int order,
+                                psScalar *x)
+{
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+            "---- p_psVectorInterpolate() begin ----\n");
+
     if (domain->type.type != range->type.type != x->type.type) {
         // XXX psError
@@ -1898,11 +1925,18 @@
 
     if (x->type.type == PS_TYPE_F32) {
-        return(p_psInterpolate1DF32(domain->data.F32, range->data.F32,
-                                    domain->n, order, x->data.F32));
+        psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+                "---- p_psVectorInterpolate() end ----\n");
+        return(psScalarAlloc(p_ps1DInterpolateF32(domain->data.F32,
+                             range->data.F32,
+                             domain->n,
+                             order,
+                             x->data.F32), PS_TYPE_F32));
     } else {
         // XXX psError: type not supported
     }
 
-    return(-1.0);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+            "---- p_psVectorInterpolate() end ----\n");
+    return(NULL);
 }
 
Index: /trunk/psLib/src/dataManip/psFunctions.h
===================================================================
--- /trunk/psLib/src/dataManip/psFunctions.h	(revision 1830)
+++ /trunk/psLib/src/dataManip/psFunctions.h	(revision 1831)
@@ -12,6 +12,6 @@
 *  @author George Gusciora, MHPCC
 *
-*  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-09-17 02:14:52 $
+*  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-09-18 01:50:45 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -409,4 +409,9 @@
                         psScalar *x);
 
+psScalar *p_psVectorInterpolate(psVector *domain,
+                                psVector *range,
+                                int order,
+                                psScalar *x);
+
 /* \} */// End of MathGroup Functions
 
Index: /trunk/psLib/src/dataManip/psMinimize.c
===================================================================
--- /trunk/psLib/src/dataManip/psMinimize.c	(revision 1830)
+++ /trunk/psLib/src/dataManip/psMinimize.c	(revision 1831)
@@ -9,6 +9,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-16 19:30:24 $
+ *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-18 01:50:45 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -844,10 +844,9 @@
     return (initialGuess);
 }
+
 /******************************************************************************
 This routine will take an procedure which calculates an arbitrary function
 and it's derivative and minimize it.
  
-GUS
- 
 XXX: Do this:
  After checking that all entries in the paramMask are 1 or 0, when
@@ -855,6 +854,4 @@
  
      A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
- 
- 
  *****************************************************************************/
 bool psMinimizeLM(psMinimization *min,
@@ -982,5 +979,5 @@
 
         oldValue = min->value;
-        newValue = func(newDeriv, NRparams, coords);
+        newValue = func(deriv, NRparams, coords);
         psTrace(".psLib.dataManip.psMinimizeLM", 4,
                 "old/new values are (%f, %f)\n", oldValue, newValue);
@@ -1029,4 +1026,238 @@
     psTrace(".psLib.dataManip.psMinimizeLM", 4,
             "---- psMinimizeLM() end (false) ----\n");
+    return(false);
+}
+
+/******************************************************************************
+This routine will take an procedure which calculates an arbitrary function
+and it's derivative and minimize the chi-squared match between that function
+at the specified coords and the specified value at those coords.
+ 
+XXX: Do this:
+ After checking that all entries in the paramMask are 1 or 0, when
+ forming the A matrix from alpha, try this:
+ 
+     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
+ *****************************************************************************/
+bool psMinimizeLMChi2(psMinimization *min,
+                      psVector *params,
+                      const psVector *paramMask,
+                      psImage *covar,
+                      const psArray *coords,
+                      const psVector *value,
+                      psMinimizeLMChi2Func func)
+{
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+            "---- psMinimizeLMChi2() begin ----\n");
+    int numData = value->n;
+    int numParams = params->n;
+    int i;
+    int j;
+    int k;
+    int l;
+    int n;
+    int p;
+    psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64);
+    psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64);
+
+    psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64);
+    psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32);
+    psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32);
+
+    psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32);
+    psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+    psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+
+    psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));
+    for (i=0;i<numData;i++) {
+        deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);
+    }
+
+    psVector *currValueVec = psVectorAlloc(value->n, PS_TYPE_F32);
+    psVector *newValueVec = psVectorAlloc(value->n, PS_TYPE_F32);
+
+    float currChi2 = 0.0;
+    float newChi2 = 0.0;
+    float lamda = 0.00005;
+
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            "min->maxIter is %d\n", min->maxIter);
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            "min->tol is %f\n", min->tol);
+
+    for (p=0;p<numParams;p++) {
+        origParams->data.F32[p] = params->data.F32[p];
+    }
+
+    min->lastDelta = HUGE;
+    min->lastDelta = 12345.0;
+    min->iter = 0;
+
+    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "------------------------------------------------------\n");
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
+
+        //
+        // Calculate the current values and chi-squared of the function.
+        //
+        currChi2 = 0.0;
+        for (n=0;n<numData;n++) {
+            currValueVec->data.F32[n] = func(deriv[n],
+                                             params,
+                                             (psVector *) coords->data[n]);
+            currChi2+= (currValueVec->data.F32[n] * currValueVec->data.F32[n]);
+        }
+
+        for (p=0;p<numParams;p++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);
+        }
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                "Current chi-squared is (%f)\n", currChi2);
+
+        //
+        // Mask elements of the derivative for each data point.
+        //
+        for (p=0;p<numParams;p++) {
+            if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
+                for (n=0;n<numData;n++) {
+                    (deriv[n])->data.F32[p] = 0.0;
+                }
+            }
+        }
+
+        //
+        // Calculate the BETA vector.
+        //
+        for (p=0;p<numParams;p++) {
+            beta->data.F64[p] = 0.0;
+            for (n=0;n<numData;n++) {
+                (beta->data.F64[p])+=
+                    (value->data.F32[n] - currValueVec->data.F32[n]) *
+                    (deriv[n])->data.F32[p];
+            }
+            // XXX: multiple by -1 here?
+            (beta->data.F64[p])*= -1.0;
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);
+        }
+
+        //
+        // Calculate the ALPHA matrix.
+        //
+        for (k=0;k<numParams;k++) {
+            for (l=0;l<numParams;l++) {
+                alpha->data.F32[k][l] = 0.0;
+                for (n=0;n<numData;n++) {
+                    alpha->data.F32[k][l]+= (deriv[n])->data.F32[k] *
+                                            (deriv[n])->data.F32[l];
+                }
+            }
+        }
+
+        //
+        // Calculate the matrix A.
+        //
+        for (j=0;j<numParams;j++) {
+            for (k=0;k<numParams;k++) {
+                if (j == k) {
+                    A->data.F64[j][k] =
+                        (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
+                } else {
+                    A->data.F64[j][k] = (double) alpha->data.F32[j][k];
+                }
+            }
+        }
+        for (j=0;j<numParams;j++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
+            for (k=0;k<numParams;k++) {
+                psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "%f ", A->data.F64[j][k]);
+            }
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
+        }
+
+        //
+        // Solve A * alpha = Beta
+        //
+        aOut = psMatrixLUD(aOut, perm, A);
+        paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);
+
+        //
+        // Mask any masked parameters.
+        //
+        for (i=0;i<numParams;i++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);
+            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
+                newParams->data.F32[i] = origParams->data.F32[i];
+            } else {
+                newParams->data.F32[i] = params->data.F32[i] -
+                                         (float) paramDeltasF64->data.F64[i];
+            }
+        }
+
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                "Calling func() with new parameters:\n");
+        for (i=0;i<numParams;i++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);
+        }
+
+
+        //
+        // Calculate new function values.
+        //
+        newChi2 = 0.0;
+        for (n=0;n<numData;n++) {
+            newValueVec->data.F32[n] = func(deriv[n], newParams,
+                                            (psVector *) coords->data[n]);
+            newChi2+= (newValueVec->data.F32[n] * newValueVec->data.F32[n]);
+        }
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);
+
+        //
+        // If the new chi-squared is lower, then keep it.
+        //
+        if (currChi2 > newChi2) {
+            min->lastDelta = currChi2 - newChi2;
+            min->value = newChi2;
+
+            // We already masked params.
+            for (i=0;i<numParams;i++) {
+                params->data.F32[i] = (float) newParams->data.F32[i];
+            }
+            lamda*= 0.1;
+        } else {
+            lamda*= 10.0;
+        }
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "lamda is %f\n", lamda);
+        min->iter++;
+    }
+    psFree(beta);
+    psFree(perm);
+    psFree(paramDeltasF64);
+    psFree(origParams);
+    psFree(newParams);
+    psFree(alpha);
+    psFree(A);
+    psFree(aOut);
+    for (i=0;i<numData;i++) {
+        psFree(deriv[i]);
+    }
+    psFree(deriv);
+    psFree(currValueVec);
+    psFree(newValueVec);
+
+    if ((min->iter < min->maxIter) ||
+            (min->lastDelta <= min->tol)) {
+        return(true);
+    }
+
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+            "---- psMinimizeLMChi2() end (false) ----\n");
     return(false);
 }
Index: /trunk/psLib/src/dataManip/psMinimize.h
===================================================================
--- /trunk/psLib/src/dataManip/psMinimize.h	(revision 1830)
+++ /trunk/psLib/src/dataManip/psMinimize.h	(revision 1831)
@@ -7,6 +7,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-08 23:32:03 $
+ *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-18 01:50:45 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -137,4 +137,8 @@
                                    const psArray *coords);
 
+typedef float (*psMinimizeLMChi2Func) (psVector *deriv,
+                                       const psVector *params,
+                                       const psVector *coords);
+
 bool psMinimizeLM(psMinimization *min,
                   psImage *covar,
@@ -143,4 +147,12 @@
                   const psArray *coords,
                   psMinimizeLMFunc func);
+
+bool psMinimizeLMChi2(psMinimization *min,
+                      psVector *params,
+                      const psVector *paramMask,
+                      psImage *covar,
+                      const psArray *coords,
+                      const psVector *value,
+                      psMinimizeLMChi2Func func);
 
 
Index: /trunk/psLib/src/math/psMinimize.c
===================================================================
--- /trunk/psLib/src/math/psMinimize.c	(revision 1830)
+++ /trunk/psLib/src/math/psMinimize.c	(revision 1831)
@@ -9,6 +9,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-16 19:30:24 $
+ *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-18 01:50:45 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -844,10 +844,9 @@
     return (initialGuess);
 }
+
 /******************************************************************************
 This routine will take an procedure which calculates an arbitrary function
 and it's derivative and minimize it.
  
-GUS
- 
 XXX: Do this:
  After checking that all entries in the paramMask are 1 or 0, when
@@ -855,6 +854,4 @@
  
      A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
- 
- 
  *****************************************************************************/
 bool psMinimizeLM(psMinimization *min,
@@ -982,5 +979,5 @@
 
         oldValue = min->value;
-        newValue = func(newDeriv, NRparams, coords);
+        newValue = func(deriv, NRparams, coords);
         psTrace(".psLib.dataManip.psMinimizeLM", 4,
                 "old/new values are (%f, %f)\n", oldValue, newValue);
@@ -1029,4 +1026,238 @@
     psTrace(".psLib.dataManip.psMinimizeLM", 4,
             "---- psMinimizeLM() end (false) ----\n");
+    return(false);
+}
+
+/******************************************************************************
+This routine will take an procedure which calculates an arbitrary function
+and it's derivative and minimize the chi-squared match between that function
+at the specified coords and the specified value at those coords.
+ 
+XXX: Do this:
+ After checking that all entries in the paramMask are 1 or 0, when
+ forming the A matrix from alpha, try this:
+ 
+     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
+ *****************************************************************************/
+bool psMinimizeLMChi2(psMinimization *min,
+                      psVector *params,
+                      const psVector *paramMask,
+                      psImage *covar,
+                      const psArray *coords,
+                      const psVector *value,
+                      psMinimizeLMChi2Func func)
+{
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+            "---- psMinimizeLMChi2() begin ----\n");
+    int numData = value->n;
+    int numParams = params->n;
+    int i;
+    int j;
+    int k;
+    int l;
+    int n;
+    int p;
+    psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64);
+    psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64);
+
+    psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64);
+    psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32);
+    psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32);
+
+    psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32);
+    psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+    psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+
+    psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));
+    for (i=0;i<numData;i++) {
+        deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);
+    }
+
+    psVector *currValueVec = psVectorAlloc(value->n, PS_TYPE_F32);
+    psVector *newValueVec = psVectorAlloc(value->n, PS_TYPE_F32);
+
+    float currChi2 = 0.0;
+    float newChi2 = 0.0;
+    float lamda = 0.00005;
+
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            "min->maxIter is %d\n", min->maxIter);
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            "min->tol is %f\n", min->tol);
+
+    for (p=0;p<numParams;p++) {
+        origParams->data.F32[p] = params->data.F32[p];
+    }
+
+    min->lastDelta = HUGE;
+    min->lastDelta = 12345.0;
+    min->iter = 0;
+
+    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "------------------------------------------------------\n");
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
+
+        //
+        // Calculate the current values and chi-squared of the function.
+        //
+        currChi2 = 0.0;
+        for (n=0;n<numData;n++) {
+            currValueVec->data.F32[n] = func(deriv[n],
+                                             params,
+                                             (psVector *) coords->data[n]);
+            currChi2+= (currValueVec->data.F32[n] * currValueVec->data.F32[n]);
+        }
+
+        for (p=0;p<numParams;p++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);
+        }
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                "Current chi-squared is (%f)\n", currChi2);
+
+        //
+        // Mask elements of the derivative for each data point.
+        //
+        for (p=0;p<numParams;p++) {
+            if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
+                for (n=0;n<numData;n++) {
+                    (deriv[n])->data.F32[p] = 0.0;
+                }
+            }
+        }
+
+        //
+        // Calculate the BETA vector.
+        //
+        for (p=0;p<numParams;p++) {
+            beta->data.F64[p] = 0.0;
+            for (n=0;n<numData;n++) {
+                (beta->data.F64[p])+=
+                    (value->data.F32[n] - currValueVec->data.F32[n]) *
+                    (deriv[n])->data.F32[p];
+            }
+            // XXX: multiple by -1 here?
+            (beta->data.F64[p])*= -1.0;
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);
+        }
+
+        //
+        // Calculate the ALPHA matrix.
+        //
+        for (k=0;k<numParams;k++) {
+            for (l=0;l<numParams;l++) {
+                alpha->data.F32[k][l] = 0.0;
+                for (n=0;n<numData;n++) {
+                    alpha->data.F32[k][l]+= (deriv[n])->data.F32[k] *
+                                            (deriv[n])->data.F32[l];
+                }
+            }
+        }
+
+        //
+        // Calculate the matrix A.
+        //
+        for (j=0;j<numParams;j++) {
+            for (k=0;k<numParams;k++) {
+                if (j == k) {
+                    A->data.F64[j][k] =
+                        (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
+                } else {
+                    A->data.F64[j][k] = (double) alpha->data.F32[j][k];
+                }
+            }
+        }
+        for (j=0;j<numParams;j++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
+            for (k=0;k<numParams;k++) {
+                psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "%f ", A->data.F64[j][k]);
+            }
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
+        }
+
+        //
+        // Solve A * alpha = Beta
+        //
+        aOut = psMatrixLUD(aOut, perm, A);
+        paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);
+
+        //
+        // Mask any masked parameters.
+        //
+        for (i=0;i<numParams;i++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);
+            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
+                newParams->data.F32[i] = origParams->data.F32[i];
+            } else {
+                newParams->data.F32[i] = params->data.F32[i] -
+                                         (float) paramDeltasF64->data.F64[i];
+            }
+        }
+
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                "Calling func() with new parameters:\n");
+        for (i=0;i<numParams;i++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);
+        }
+
+
+        //
+        // Calculate new function values.
+        //
+        newChi2 = 0.0;
+        for (n=0;n<numData;n++) {
+            newValueVec->data.F32[n] = func(deriv[n], newParams,
+                                            (psVector *) coords->data[n]);
+            newChi2+= (newValueVec->data.F32[n] * newValueVec->data.F32[n]);
+        }
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);
+
+        //
+        // If the new chi-squared is lower, then keep it.
+        //
+        if (currChi2 > newChi2) {
+            min->lastDelta = currChi2 - newChi2;
+            min->value = newChi2;
+
+            // We already masked params.
+            for (i=0;i<numParams;i++) {
+                params->data.F32[i] = (float) newParams->data.F32[i];
+            }
+            lamda*= 0.1;
+        } else {
+            lamda*= 10.0;
+        }
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "lamda is %f\n", lamda);
+        min->iter++;
+    }
+    psFree(beta);
+    psFree(perm);
+    psFree(paramDeltasF64);
+    psFree(origParams);
+    psFree(newParams);
+    psFree(alpha);
+    psFree(A);
+    psFree(aOut);
+    for (i=0;i<numData;i++) {
+        psFree(deriv[i]);
+    }
+    psFree(deriv);
+    psFree(currValueVec);
+    psFree(newValueVec);
+
+    if ((min->iter < min->maxIter) ||
+            (min->lastDelta <= min->tol)) {
+        return(true);
+    }
+
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+            "---- psMinimizeLMChi2() end (false) ----\n");
     return(false);
 }
Index: /trunk/psLib/src/math/psMinimize.h
===================================================================
--- /trunk/psLib/src/math/psMinimize.h	(revision 1830)
+++ /trunk/psLib/src/math/psMinimize.h	(revision 1831)
@@ -7,6 +7,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-08 23:32:03 $
+ *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-18 01:50:45 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -137,4 +137,8 @@
                                    const psArray *coords);
 
+typedef float (*psMinimizeLMChi2Func) (psVector *deriv,
+                                       const psVector *params,
+                                       const psVector *coords);
+
 bool psMinimizeLM(psMinimization *min,
                   psImage *covar,
@@ -143,4 +147,12 @@
                   const psArray *coords,
                   psMinimizeLMFunc func);
+
+bool psMinimizeLMChi2(psMinimization *min,
+                      psVector *params,
+                      const psVector *paramMask,
+                      psImage *covar,
+                      const psArray *coords,
+                      const psVector *value,
+                      psMinimizeLMChi2Func func);
 
 
Index: /trunk/psLib/src/math/psPolynomial.c
===================================================================
--- /trunk/psLib/src/math/psPolynomial.c	(revision 1830)
+++ /trunk/psLib/src/math/psPolynomial.c	(revision 1831)
@@ -7,6 +7,6 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-17 02:14:52 $
+ *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-18 01:50:45 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1809,10 +1809,10 @@
 p_psInterpolate1D(): This routine will take as input n-element floating
 point arrays domain and range, and the x value, assumed to lie with the
-domain vector.  It produces as output the n-order LaGrange interpolated
+domain vector.  It produces as output the (n-1)-order LaGrange interpolated
 value of x.
  
 XXX: do we error check for non-distinct domain values?
  *****************************************************************************/
-float p_psFullInterpolate1DF32(float *domain,
+float p_ps1DFullInterpolateF32(float *domain,
                                float *range,
                                int n,
@@ -1822,21 +1822,39 @@
     int m;
     static psVector *p = NULL;
-    psVectorRecycle(p, n, PS_TYPE_F32);
+    p = psVectorRecycle(p, n, PS_TYPE_F32);
     p_psMemSetPersistent(p, true);
+    p_psMemSetPersistent(p->data.F32, true);
+
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
+            "---- p_ps1DFullInterpolateF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n);
+
+    for (i=0;i<n;i++) {
+        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                "domain/range is (%f %f)\n", domain[i], range[i]);
+    }
+
+    for (i=0;i<n;i++) {
+        p->data.F32[i] = range[i];
+        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
+
+    }
 
     // From NR, during each iteration of the m loop, we are computing the
     // p_{i ... i+m} terms.
-    for (m=0;m<n;m++) {
+    for (m=1;m<n;m++) {
         for (i=0;i<n-m;i++) {
-            if (m == 0) {
-                p->data.F32[i] = range[i];
-            } else {
-                // From NR: we are computing P_{i ... i+m}
-                p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) +
-                                  ((range[i]-x) * p->data.F32[i+1])) /
-                                 (domain[i] - domain[i+m]);
-            }
-        }
-    }
+            // From NR: we are computing P_{i ... i+m}
+            p->data.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) +
+                              ((domain[i]-x) * p->data.F32[i+1])) /
+                             (domain[i] - domain[i+m]);
+            //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]);
+            psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                    "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
+        }
+    }
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
+            "---- p_ps1DFullInterpolateF32() end ----\n");
+
     return(p->data.F32[0]);
 }
@@ -1844,5 +1862,5 @@
 
 // This is a
-float p_psInterpolate1DF32(float *domain,
+float p_ps1DInterpolateF32(float *domain,
                            float *range,
                            int n,
@@ -1854,5 +1872,9 @@
     int origin;
 
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
+            "---- p_ps1DInterpolateF32() begin ----\n");
+
     binNum = VectorBinDisectF32(domain, n, x);
+
     if (0 == numIntPoints%2) {
         origin = binNum - ((numIntPoints/2) - 1);
@@ -1871,5 +1893,7 @@
     }
 
-    return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x));
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
+            "---- p_ps1DInterpolateF32() end ----\n");
+    return(p_ps1DFullInterpolateF32(&domain[origin], &range[origin], order+1, x));
 }
 
@@ -1882,9 +1906,12 @@
 XXX: This stuff does not work with a mask.
  *****************************************************************************/
-float p_psInterpolate1D(psVector *domain,
-                        psVector *range,
-                        int order,
-                        psScalar *x)
-{
+psScalar *p_psVectorInterpolate(psVector *domain,
+                                psVector *range,
+                                int order,
+                                psScalar *x)
+{
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+            "---- p_psVectorInterpolate() begin ----\n");
+
     if (domain->type.type != range->type.type != x->type.type) {
         // XXX psError
@@ -1898,11 +1925,18 @@
 
     if (x->type.type == PS_TYPE_F32) {
-        return(p_psInterpolate1DF32(domain->data.F32, range->data.F32,
-                                    domain->n, order, x->data.F32));
+        psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+                "---- p_psVectorInterpolate() end ----\n");
+        return(psScalarAlloc(p_ps1DInterpolateF32(domain->data.F32,
+                             range->data.F32,
+                             domain->n,
+                             order,
+                             x->data.F32), PS_TYPE_F32));
     } else {
         // XXX psError: type not supported
     }
 
-    return(-1.0);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+            "---- p_psVectorInterpolate() end ----\n");
+    return(NULL);
 }
 
Index: /trunk/psLib/src/math/psPolynomial.h
===================================================================
--- /trunk/psLib/src/math/psPolynomial.h	(revision 1830)
+++ /trunk/psLib/src/math/psPolynomial.h	(revision 1831)
@@ -12,6 +12,6 @@
 *  @author George Gusciora, MHPCC
 *
-*  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-09-17 02:14:52 $
+*  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-09-18 01:50:45 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -409,4 +409,9 @@
                         psScalar *x);
 
+psScalar *p_psVectorInterpolate(psVector *domain,
+                                psVector *range,
+                                int order,
+                                psScalar *x);
+
 /* \} */// End of MathGroup Functions
 
Index: /trunk/psLib/src/math/psSpline.c
===================================================================
--- /trunk/psLib/src/math/psSpline.c	(revision 1830)
+++ /trunk/psLib/src/math/psSpline.c	(revision 1831)
@@ -7,6 +7,6 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-17 02:14:52 $
+ *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-18 01:50:45 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1809,10 +1809,10 @@
 p_psInterpolate1D(): This routine will take as input n-element floating
 point arrays domain and range, and the x value, assumed to lie with the
-domain vector.  It produces as output the n-order LaGrange interpolated
+domain vector.  It produces as output the (n-1)-order LaGrange interpolated
 value of x.
  
 XXX: do we error check for non-distinct domain values?
  *****************************************************************************/
-float p_psFullInterpolate1DF32(float *domain,
+float p_ps1DFullInterpolateF32(float *domain,
                                float *range,
                                int n,
@@ -1822,21 +1822,39 @@
     int m;
     static psVector *p = NULL;
-    psVectorRecycle(p, n, PS_TYPE_F32);
+    p = psVectorRecycle(p, n, PS_TYPE_F32);
     p_psMemSetPersistent(p, true);
+    p_psMemSetPersistent(p->data.F32, true);
+
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
+            "---- p_ps1DFullInterpolateF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n);
+
+    for (i=0;i<n;i++) {
+        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                "domain/range is (%f %f)\n", domain[i], range[i]);
+    }
+
+    for (i=0;i<n;i++) {
+        p->data.F32[i] = range[i];
+        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
+
+    }
 
     // From NR, during each iteration of the m loop, we are computing the
     // p_{i ... i+m} terms.
-    for (m=0;m<n;m++) {
+    for (m=1;m<n;m++) {
         for (i=0;i<n-m;i++) {
-            if (m == 0) {
-                p->data.F32[i] = range[i];
-            } else {
-                // From NR: we are computing P_{i ... i+m}
-                p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) +
-                                  ((range[i]-x) * p->data.F32[i+1])) /
-                                 (domain[i] - domain[i+m]);
-            }
-        }
-    }
+            // From NR: we are computing P_{i ... i+m}
+            p->data.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) +
+                              ((domain[i]-x) * p->data.F32[i+1])) /
+                             (domain[i] - domain[i+m]);
+            //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]);
+            psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                    "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
+        }
+    }
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
+            "---- p_ps1DFullInterpolateF32() end ----\n");
+
     return(p->data.F32[0]);
 }
@@ -1844,5 +1862,5 @@
 
 // This is a
-float p_psInterpolate1DF32(float *domain,
+float p_ps1DInterpolateF32(float *domain,
                            float *range,
                            int n,
@@ -1854,5 +1872,9 @@
     int origin;
 
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
+            "---- p_ps1DInterpolateF32() begin ----\n");
+
     binNum = VectorBinDisectF32(domain, n, x);
+
     if (0 == numIntPoints%2) {
         origin = binNum - ((numIntPoints/2) - 1);
@@ -1871,5 +1893,7 @@
     }
 
-    return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x));
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
+            "---- p_ps1DInterpolateF32() end ----\n");
+    return(p_ps1DFullInterpolateF32(&domain[origin], &range[origin], order+1, x));
 }
 
@@ -1882,9 +1906,12 @@
 XXX: This stuff does not work with a mask.
  *****************************************************************************/
-float p_psInterpolate1D(psVector *domain,
-                        psVector *range,
-                        int order,
-                        psScalar *x)
-{
+psScalar *p_psVectorInterpolate(psVector *domain,
+                                psVector *range,
+                                int order,
+                                psScalar *x)
+{
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+            "---- p_psVectorInterpolate() begin ----\n");
+
     if (domain->type.type != range->type.type != x->type.type) {
         // XXX psError
@@ -1898,11 +1925,18 @@
 
     if (x->type.type == PS_TYPE_F32) {
-        return(p_psInterpolate1DF32(domain->data.F32, range->data.F32,
-                                    domain->n, order, x->data.F32));
+        psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+                "---- p_psVectorInterpolate() end ----\n");
+        return(psScalarAlloc(p_ps1DInterpolateF32(domain->data.F32,
+                             range->data.F32,
+                             domain->n,
+                             order,
+                             x->data.F32), PS_TYPE_F32));
     } else {
         // XXX psError: type not supported
     }
 
-    return(-1.0);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+            "---- p_psVectorInterpolate() end ----\n");
+    return(NULL);
 }
 
Index: /trunk/psLib/src/math/psSpline.h
===================================================================
--- /trunk/psLib/src/math/psSpline.h	(revision 1830)
+++ /trunk/psLib/src/math/psSpline.h	(revision 1831)
@@ -12,6 +12,6 @@
 *  @author George Gusciora, MHPCC
 *
-*  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-09-17 02:14:52 $
+*  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-09-18 01:50:45 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -409,4 +409,9 @@
                         psScalar *x);
 
+psScalar *p_psVectorInterpolate(psVector *domain,
+                                psVector *range,
+                                int order,
+                                psScalar *x);
+
 /* \} */// End of MathGroup Functions