Index: trunk/psLib/src/math/Makefile.am
===================================================================
--- trunk/psLib/src/math/Makefile.am	(revision 36180)
+++ trunk/psLib/src/math/Makefile.am	(revision 36222)
@@ -7,4 +7,5 @@
 	psUnaryOp.c \
 	psBinaryOp.c \
+	psBinomialCoeff.c \
 	psClip.c \
 	psCompare.c \
@@ -36,4 +37,5 @@
 	psUnaryOp.h \
 	psBinaryOp.h \
+	psBinomialCoeff.h \
 	psClip.h \
 	psCompare.h \
Index: trunk/psLib/src/math/psBinomialCoeff.c
===================================================================
--- trunk/psLib/src/math/psBinomialCoeff.c	(revision 36222)
+++ trunk/psLib/src/math/psBinomialCoeff.c	(revision 36222)
@@ -0,0 +1,30 @@
+/** @file psBinomialCoeff.c
+ *
+ *  @brief Calculate binomial coefficient
+ *
+ *  @author Chris Waters, IfA
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include "psError.h"
+
+#include "psBinomialCoeff.h"
+
+psS64 psBinomialCoeff( psS64 N,
+		       psS64 k) {
+  psS64 i;
+  psS64 b = 1;
+  if (N == k) { return (b); }
+  if (k == 0) { return (b); }
+
+  for (i = 1; i <= k; i++) {
+    b *= (N - (k - i))/i;
+  }
+
+  return(b);
+}
+
Index: trunk/psLib/src/math/psBinomialCoeff.h
===================================================================
--- trunk/psLib/src/math/psBinomialCoeff.h	(revision 36222)
+++ trunk/psLib/src/math/psBinomialCoeff.h	(revision 36222)
@@ -0,0 +1,8 @@
+#ifndef PS_BINCOEFF_H
+#define PS_BINCOEFF_H
+
+// Calculate the binomial coefficient
+psS64 psBinomialCoeff( psS64 N,
+		       psS64 k);
+
+#endif
Index: trunk/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- trunk/psLib/src/math/psMinimizePolyFit.c	(revision 36180)
+++ trunk/psLib/src/math/psMinimizePolyFit.c	(revision 36222)
@@ -42,7 +42,9 @@
 #include "psLogMsg.h"
 #include "psMathUtils.h"
+#include "psBinomialCoeff.h"
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
+#define CZW 0
 
 # define USE_GAUSS_JORDAN 1
@@ -603,5 +605,26 @@
     bool status = false;
     if (USE_GAUSS_JORDAN) {
+
+#if (CZW)
+      	printf("CZW: about to do GJ: %d\n",status);
+	for (psS32 k = 0; k < nTerm; k++) {
+	  printf("CZW: %d %f \t",k,B->data.F64[k]);
+	  for (psS32 kk = 0; kk < nTerm; kk++) {
+	    printf(" %f ",(A->data.F64[k][kk]));
+	  }
+	  printf("\n");
+	}
+#endif
         status = psMatrixGJSolve(A, B);
+#if (CZW)
+	printf("CZW: just did GJ: %d\n",status);
+	for (psS32 k = 0; k < nTerm; k++) {
+	  printf("CZW: %d %f \t",k,B->data.F64[k]);
+	  for (psS32 kk = 0; kk < nTerm; kk++) {
+	    printf(" %f ",(A->data.F64[k][kk]));
+	  }
+	  printf("\n");
+	}
+#endif
     } else {
         status = psMatrixLUSolve(A, B);
@@ -676,10 +699,77 @@
     bool result = true;
 
+    // Define values that may be used by PS_POLYNOMIAL_ORD form.
+    bool scale = false;
+    double median = 0.0;
+    double sigma = 1.0;
+    psVector *sorted = NULL;
+    psVector *z64 = NULL;
+    
     switch (poly->type) {
     case PS_POLYNOMIAL_ORD:
-        result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
-        if (!result) {
-            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
-        }
+      if ((f64->n < 10000)&&(poly->nX > 1)) {
+	scale = true;
+	sorted = psVectorSort(NULL,x64);
+	median = sorted->data.F64[sorted->n / 2];
+	// CZW: I'm not bothering to scale this because it doesn't really matter.
+	sigma  = (sorted->data.F64[3 * sorted->n / 4] - sorted->data.F64[sorted->n / 4]); 
+	psFree(sorted);
+	// I can't see a way to not clobber x if it's already F64, so make a copy.x
+	z64 = psVectorCopy(NULL,x64,PS_TYPE_F64);
+	psBinaryOp(z64,z64,"-",psScalarAlloc(median,PS_TYPE_F64));
+	psBinaryOp(z64,z64,"/",psScalarAlloc(sigma,PS_TYPE_F64));
+
+#if (CZW)
+	printf("poly1d: Scale parameters: %f %f\n",median,sigma);
+#endif
+
+	
+	result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, z64);
+      }
+      else {
+	result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
+      }
+        
+      if (!result) {
+	psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning NULL.\n");
+      }
+      
+      if (scale) {
+	psFree(z64); // Done with this.
+
+	// Undo scaling in the polynomial values.
+	psF64 *Zcoeff = psAlloc((1 + poly->nX) * sizeof(psF64 *));
+	psF64 *ZcoeffErr = psAlloc((1 + poly->nX) * sizeof(psF64 *));
+	
+	for (psS32 i = 0; i <= poly->nX; i++) {
+	  Zcoeff[i] = poly->coeff[i];
+	  ZcoeffErr[i] = poly->coeffErr[i];
+#if (CZW)
+	  printf("poly1d: fit parameters: %d %f %f\n",
+		 i,poly->coeff[i],poly->coeffErr[i]);
+#endif
+	}
+	for (psS32 i = 0; i <= poly->nX; i++) {
+	  poly->coeff[i] = 0.0;
+	  poly->coeffErr[i] = 0.0;
+	  
+	  for (psS32 j = 0; j <= poly->nX; j++) {
+#if (CZW)
+	    printf("        %d %d %f %f %f %f => %f\n",
+		   i,j,Zcoeff[j],pow(1.0 / sigma,j) * pow(-1,j - i),pow(median,j - i),1.0 * psBinomialCoeff(j,i),
+		   Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j  -i) * pow(median,j - i) * 1.0 * psBinomialCoeff(j,i)
+		   );
+#endif
+	    poly->coeff[i] += Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - i) * psBinomialCoeff(j,i);
+	    poly->coeffErr[i] += pow(ZcoeffErr[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - 1) * psBinomialCoeff(j,i),2);
+	  }
+	  poly->coeffErr[i] = sqrt(poly->coeffErr[i]);
+#if (CZW)
+	  printf("poly1d: unscaled parameters: %d %f %f\n",
+		 i,poly->coeff[i], poly->coeffErr[i]);
+#endif
+	}
+      }
+	
         break;
     case PS_POLYNOMIAL_CHEB:
Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 36180)
+++ trunk/psLib/src/math/psStats.c	(revision 36222)
@@ -141,5 +141,5 @@
 
 // Debug information
-#define CZW 0
+#define CZW 1
 
 /*****************************************************************************/
@@ -232,4 +232,5 @@
         }
         count++;
+
     }
     if (errors) {
@@ -424,4 +425,5 @@
     for (long i = 0; i < myVector->n; i++) {
         // Check if the data is with the specified range
+
         if (!isfinite(data[i]))
             continue;
@@ -1043,6 +1045,6 @@
 
 
-#if (CZW)
-	printf("CZW (%d): median %f sigma %f delta: %f \t %f %f %f %f %f %f %f \t %f %f %f %f %f %f %f\n",
+#if (CZW && 0)
+	printf("CZW (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n",
 	       iterate,
 	   stats->robustMedian,stats->robustStdev,
@@ -1060,8 +1062,8 @@
 	   cumulative->nums->data.F32[binMedian+1],
 	   cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]);
-	PS_VECTOR_PRINT_F32(histogram->bounds);
-	PS_VECTOR_PRINT_F32(histogram->nums);
-	PS_VECTOR_PRINT_F32(cumulative->bounds);
-	PS_VECTOR_PRINT_F32(cumulative->nums);
+	//	PS_VECTOR_PRINT_F32(histogram->bounds);
+	//	PS_VECTOR_PRINT_F32(histogram->nums);
+	//	PS_VECTOR_PRINT_F32(cumulative->bounds);
+	//	PS_VECTOR_PRINT_F32(cumulative->nums);
 #endif
 
@@ -1276,5 +1278,5 @@
             return true;
         }
-	//	printf("BINS: %ld %f %f %f %f\n",stats->robustN50,guessStdev,binSize,min,max);
+
         // Calculate the number of bins.
         // XXX can we calculate the binMin, binMax **before** building this histogram?
@@ -1286,5 +1288,4 @@
         psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
 
-	//	printf("BINS2: %ld %f %f %f %f %ld\n",stats->robustN50,guessStdev,binSize,min,max,numBins);
 	
         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
@@ -1385,5 +1386,5 @@
             }
             y->n = x->n = j;
-
+	    
             // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
             // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
@@ -1397,4 +1398,5 @@
                 psErrorClear();
                 COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
+
                 psFree(poly);
                 psFree(histogram);
@@ -1478,4 +1480,12 @@
             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
             bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
+#if (CZW && 0)
+	    for (long i = 0; i < x->n; i++) {
+	      printf("CZW: Dcheck: %ld %f %f %f\n",
+		     i,x->data.F32[i],y->data.F32[i],
+		     poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
+		     poly->coeff[2] * pow(x->data.F32[i],2));
+	    }
+#endif
             psFree(x);
             psFree(y);
@@ -1493,4 +1503,5 @@
             fullfitStdev = sqrt(-0.5/poly->coeff[2]);
             fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev);
+
 #ifndef PS_NO_TRACE
             psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
@@ -1515,4 +1526,5 @@
             }
 
+	    
             psFree (poly);
         }
@@ -1537,7 +1549,16 @@
 	    done = true;
 	}
-#if (CZW) 
-	printf("CZW IN FITTED? low %f %f full %f %f robust %f %f final %f %f\n",
-	       lowfitMean,lowfitStdev,fullfitMean,fullfitStdev,stats->robustMedian,stats->robustStdev,
+
+	
+#if (CZW && 0)
+	printf("CZW IN FITTED: iter   %d %f \n"
+	       "               low    %f %f \n"
+	       "               full   %f %f \n"
+	       "               robust %f %f \n"
+	       "               final  %f %f\n",
+	       iteration,minValueSym,
+	       lowfitMean,lowfitStdev,
+	       fullfitMean,fullfitStdev,
+	       stats->robustMedian,stats->robustStdev,
 	       guessMean,guessStdev);
 #endif
@@ -2423,12 +2444,4 @@
         psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
 
-#if (CZW)
-	printf("  polyin: %f %f %f == %f %f %f\n",
-	       x->data.F64[0], x->data.F64[1], x->data.F64[2],
-	       y->data.F64[0], y->data.F64[1], y->data.F64[2]);
-	printf("  rawpolyin: %f %f %f == %f %f %f\n",
-	       xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum + 1],
-	       y->data.F64[0], y->data.F64[1], y->data.F64[2]);
-#endif
 
         // Ensure that the y value lies within range of the y values.
@@ -2469,11 +2482,4 @@
                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
-#if (CZW)
-	printf("  poly: %f %f %f fit: %f %f %f\n",
-	       myPoly->coeff[0],myPoly->coeff[1],myPoly->coeff[2],
-	       (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
-	       (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
-	       (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
-#endif
 
         psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
@@ -2498,8 +2504,4 @@
 	//        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
 	tmpFloat = binValue;
-		
-#if (CZW)
-	printf("   internal median: %f %f\n",tmpFloat,binValue);
-#endif
 	
     } else {
@@ -2556,36 +2558,39 @@
     )
 {
+
 # if (1)
-    double Sx = (yVec->data.F32[binNum - 2] + yVec->data.F32[binNum - 1] +
-		 yVec->data.F32[binNum - 0] +
-		 yVec->data.F32[binNum + 1] + yVec->data.F32[binNum + 2]);
-
-    double YM = 0.5*(xVec->data.F32[binNum - 2] + xVec->data.F32[binNum - 1]);
-    double Ym = 0.5*(xVec->data.F32[binNum - 1] + xVec->data.F32[binNum - 0]);
-    double Yo = 0.5*(xVec->data.F32[binNum - 0] + xVec->data.F32[binNum + 1]);
-    double Yp = 0.5*(xVec->data.F32[binNum + 1] + xVec->data.F32[binNum + 2]);
-    double YP = 0.5*(xVec->data.F32[binNum + 2] + xVec->data.F32[binNum + 3]);
-
-    double Sy = YM + Ym + Yo + Yp + YP;
-
-    double Sxx = (PS_SQR(yVec->data.F32[binNum - 2]) + PS_SQR(yVec->data.F32[binNum - 1]) +
-		  PS_SQR(yVec->data.F32[binNum - 0]) +
-		  PS_SQR(yVec->data.F32[binNum + 1]) + PS_SQR(yVec->data.F32[binNum + 2]));
-
-    double Sxy = (yVec->data.F32[binNum - 2] * YM +
-		  yVec->data.F32[binNum - 1] * Ym +
-		  yVec->data.F32[binNum - 0] * Yo +
-		  yVec->data.F32[binNum + 1] * Yp +
-		  yVec->data.F32[binNum + 2] * YP);
-
-    double Det = 5*Sxx - Sx*Sx;
-    if (Det == 0.0) return NAN;
-
-    double C0 = (Sy*Sxx - Sx*Sxy) / Det;
-    double C1 = (Sxy*5 - Sx*Sy) / Det;
-
-    double value = C0 + yVal*C1;
-
-    return value;
+# define HALF_SIZE 2
+  double Sx = 0.0;
+
+  double Sy = 0.0;
+  double Sxx = 0.0;
+  double Sxy = 0.0;
+  double deltaY = 0.0;
+  int N = 0;
+
+  for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++) {
+    if ((u >= 0)&&(u < yVec->n)) {
+      if (u+1 < xVec->n) {
+	Sx += yVec->data.F32[u];
+	Sxx += PS_SQR(yVec->data.F32[u]);
+
+	deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]);
+	Sy += deltaY;
+	Sxy += yVec->data.F32[u] * deltaY;
+	N += 1;
+      }
+    }
+  }
+  double Det = N * Sxx - Sx * Sx;
+  if (Det == 0.0) return NAN;
+  if (N == 0) return NAN;
+
+  double C0 = (Sy*Sxx - Sx*Sxy) / Det;
+  double C1 = (Sxy*N - Sx*Sy) / Det;
+  
+  double value = C0 + yVal*C1;
+  return value;
+  
+  
 # else
     psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
@@ -2682,7 +2687,4 @@
 	tmpFloat = binValue;
 		
-#if (CZW)
-	printf("   internal median: %f %f\n",tmpFloat,binValue);
-#endif
 	
     } else {