Index: trunk/tests/test_OrthoPoly1D.cc
===================================================================
--- trunk/tests/test_OrthoPoly1D.cc	(revision 875)
+++ trunk/tests/test_OrthoPoly1D.cc	(revision 876)
@@ -1,63 +1,93 @@
 #include <vector>
 
 #include "UnitTest++.h"
 
 #include "test_utils.hh"
 
 #include "npstat/nm/Matrix.hh"
 #include "npstat/nm/ArrayND.hh"
 #include "npstat/nm/OrthoPoly1D.hh"
+#include "npstat/nm/MathUtils.hh"
 
 using namespace npstat;
 using namespace std;
 
 namespace {
     TEST(OrthoPoly1D_constructors)
     {
         ArrayND<double> weight(makeShape(50));
         weight.constFill(1.0);
 
         OrthoPoly1D ortho1(3, weight.data(), weight.length(), 0.5);
         OrthoPoly1D ortho2(ortho1);
         CHECK(ortho1 == ortho2);
 
         weight(25) = 0.0;
         OrthoPoly1D ortho3(3, weight.data(), weight.length(), 0.5);
         CHECK(ortho1 != ortho3);
 
         ortho3 = ortho1;
         CHECK(ortho1 == ortho3);
     }
 
     TEST(OrthoPoly1D_globalFilter)
     {
         const unsigned n = 128;
         const double step = 0.3;
         const double coeffs[] = {0.8, 1.1, 0.4, 0.7, 0.3};
         const unsigned maxdeg = sizeof(coeffs)/sizeof(coeffs[0])-1;
 
         std::vector<double> weight(n, 1.0);
         weight[1] = 0.5;
         weight[2] = 0.1;
         weight[n - 1] = 0.2345;
         weight[n - 2] = 0.7;
 
         OrthoPoly1D poly(maxdeg, &weight[0], n, step);
         Matrix<double> gf(n, n);
         poly.globalFilter(coeffs, maxdeg, gf.data(), gf.length());
 
         Matrix<double> ref(n, n, 0);
         for (unsigned i=0; i<n; ++i)
             for (unsigned j=0; j<n; ++j)
                 for (unsigned deg=0; deg<=maxdeg; ++deg)
                     ref[i][j] += coeffs[deg]*poly.poly(deg,i)*poly.poly(deg,j)*poly.weight(j);
         ref *= step;
 
         CHECK_CLOSE(0.0, (gf - ref).maxAbsValue(), 1.0e-12);
 
         std::vector<double> diag(n);
         poly.globalFilterDiag(coeffs, maxdeg, &diag[0], diag.size());
         for (unsigned i=0; i<n; ++i)
             CHECK_CLOSE(gf[i][i], diag[i], 1.0e-12);
     }
+
+    TEST(OrthoPoly1D_derivative_fit)
+    {
+        const unsigned n = 6;
+        const unsigned maxdeg = 2;
+
+        std::vector<double> weight(n, 1.0);
+        OrthoPoly1D poly(maxdeg, &weight[0], n, 1.0);
+
+        double coeffs[maxdeg+1] = {3.0, -1.0, 0.7};
+        const double expectedDer2 = 2.0*coeffs[2];
+        std::vector<double> data(n);
+        for (unsigned i=0; i<n; ++i)
+            data[i] = polySeriesSum(coeffs, maxdeg, i);
+
+        std::vector<double> series(n);
+        poly.calculateCoeffs(&data[0], n, coeffs, maxdeg);
+        for (unsigned i=0; i<n; ++i)
+            series[i] = poly.series(coeffs, maxdeg, i);
+        
+        for (unsigned i=1; i<n-1; ++i)
+        {
+            const double secondDeriv = (data[i+1]-data[i]) - (data[i]-data[i-1]);
+            CHECK_CLOSE(expectedDer2, secondDeriv, 1.0e-14);
+
+            const double secondDer2 = (series[i+1]-series[i]) - (series[i]-series[i-1]);
+            CHECK_CLOSE(expectedDer2, secondDer2, 1.0e-14);
+        }
+    }
 }