Index: trunk/tests/test_ClassicalOrthoPolys1D.cc
===================================================================
--- trunk/tests/test_ClassicalOrthoPolys1D.cc	(revision 573)
+++ trunk/tests/test_ClassicalOrthoPolys1D.cc	(revision 574)
@@ -1,170 +1,171 @@
 #include <cmath>
 
 #include "UnitTest++.h"
 #include "test_utils.hh"
 
 #include "npstat/nm/ClassicalOrthoPolys1D.hh"
 #include "npstat/nm/GaussLegendreQuadrature.hh"
+#include "npstat/nm/FejerQuadrature.hh"
 
 #include "npstat/rng/MersenneTwister.hh"
 
 using namespace npstat;
 using namespace std;
 
 namespace {
     TEST(LegendreOrthoPoly1D_values)
     {
         const double eps = 1.0e-15;
 
         LegendreOrthoPoly1D leg;
 
         CHECK(leg.weight(0.123) == 0.5);
         CHECK(leg.weight(10.0) == 0.0);
 
         unsigned deg = 0;
         CHECK_CLOSE(1.0, leg.poly(deg, 0.3), eps);
 
         deg = 3;
         CHECK_CLOSE(-71/343.0L*sqrtl(2*deg+1.0L), leg.poly(deg, 1/7.0L), eps);
 
         deg = 7;
         CHECK_CLOSE(-326569/1771561.0L*sqrtl(2*deg+1.0L), leg.poly(deg, 1/11.0L), eps);
 
         GaussLegendreQuadrature glq(10);
         for (unsigned i=0; i<5; ++i)
             for (unsigned j=0; j<5; ++j)
                 CHECK_CLOSE((i == j ? 1.0 : 0.0), leg.empiricalKroneckerDelta(glq, i, j), eps);
 
         CHECK_CLOSE(3374/715.0L, leg.integratePoly(3, 4), eps);
         unsigned degs[3] = {3, 4, 5};
         CHECK_CLOSE(1.051943782704350297, leg.jointIntegral(degs, 3), eps);
     }
 
     TEST(JacobiOrthoPoly1D_values)
     {
         const double eps = 1.0e-15;
 
         JacobiOrthoPoly1D jac(3, 4);
 
         long double x = 1/3.0L;
         CHECK_CLOSE(35/32.0L*powl(1-x,3)*powl(1+x,4), jac.weight(x), eps);
 
         GaussLegendreQuadrature glq(10);
         for (unsigned i=0; i<5; ++i)
             for (unsigned j=0; j<5; ++j)
             {
                 CHECK_CLOSE((i == j ? 1.0 : 0.0), jac.empiricalKroneckerDelta(glq, i, j), eps);
                 unsigned degs[2];
                 degs[0] = i;
                 degs[1] = j;
                 CHECK_CLOSE((i == j ? 1.0 : 0.0), jac.jointAverage(glq, degs, 2), eps);
             }
 
         JacobiOrthoPoly1D jac2(2, 5);
         for (unsigned i=0; i<5; ++i)
             for (unsigned j=0; j<5; ++j)
             {
                 CHECK_CLOSE((i == j ? 1.0 : 0.0), jac2.empiricalKroneckerDelta(glq, i, j), eps);
                 unsigned degs[2];
                 degs[0] = i;
                 degs[1] = j;
                 CHECK_CLOSE((i == j ? 1.0 : 0.0), jac2.jointAverage(glq, degs, 2), eps);
             }
 
         JacobiOrthoPoly1D jac3(0, 0);
         CHECK_CLOSE(0.5, jac3.weight(0.123), eps);
         CHECK(jac3.weight(10.0) == 0.0);
 
         unsigned deg = 7;
         CHECK_CLOSE(-326569/1771561.0L*sqrtl(2*deg+1.0L), jac3.poly(deg, 1/11.0L), eps);
 
         MersenneTwister rng;
         const unsigned npoints = 64;
         const unsigned maxdeg = 7;
         const unsigned ntries = 5;
         std::vector<double> points(npoints);
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], npoints, npoints);
             for (unsigned ipt=0; ipt<npoints; ++ipt)
             {
                 points[ipt] *= 2.0;
                 points[ipt] -= 1.0;
             }
             const Matrix<double> avemat(jac.sampleProductAverages(&points[0], npoints, maxdeg));
             std::vector<double> averages(maxdeg + 1);
             jac.sampleAverages(&points[0], npoints, &averages[0], maxdeg);
             CHECK_EQUAL(maxdeg + 1, avemat.nRows());
             CHECK_EQUAL(maxdeg + 1, avemat.nColumns());
             for (unsigned i=0; i<=maxdeg; ++i)
             {
                 CHECK_CLOSE(averages[i], avemat[i][0], eps);
                 CHECK_CLOSE(averages[i], avemat[0][i], eps);
                 for (unsigned j=0; j<=maxdeg; ++j)
                 {
                     unsigned degs[2];
                     degs[0] = i;
                     degs[1] = j;
                     const double ave = jac.jointSampleAverage(&points[0], npoints, degs, 2);
                     CHECK_CLOSE(ave, avemat[i][j], eps);
                     CHECK_CLOSE(ave, avemat[j][i], eps);
                 }
             }
 
             // Test "twoJointSampleAverages" method
             const double degmax = maxdeg + 0.99;
             unsigned degs0[3];
             for (unsigned i=0; i<3; ++i)
                 degs0[i] = rng()*degmax;
             unsigned degs1[4];
             for (unsigned i=0; i<4; ++i)
                 degs1[i] = rng()*degmax;
             const double ave0 = jac.jointSampleAverage(&points[0], npoints, degs0, 3);
             const double ave1 = jac.jointSampleAverage(&points[0], npoints, degs1, 4);
             const std::pair<double, double>& ave2 = jac.twoJointSampleAverages(
                 &points[0], npoints, degs0, 3, degs1, 4);
             CHECK_CLOSE(ave0, ave2.first, eps);
             CHECK_CLOSE(ave1, ave2.second, eps);
         }
     }
 
     TEST(ChebyshevOrthoPoly2nd_values)
     {
         const double epsVal = 1.0e-15;
         const double epsKron = 1.0e-7;
         ChebyshevOrthoPoly2nd cheb;
 
         unsigned deg = 3;
         CHECK_CLOSE(-0.54810495626822157434L, cheb.poly(deg, 1/7.0L), epsVal);
 
         deg = 7;
         CHECK_CLOSE(-0.66835314371696127673L, cheb.poly(deg, 1/11.0L), epsVal);
 
         GaussLegendreQuadrature glq(1024);
         for (unsigned i=0; i<5; ++i)
             for (unsigned j=0; j<5; ++j)
                 CHECK_CLOSE((i == j ? 1.0 : 0.0), cheb.empiricalKroneckerDelta(glq, i, j), epsKron);
     }
 
     TEST(ChebyshevOrthoPoly1st_values)
     {
         const double epsVal = 1.0e-15;
-        const double epsKron = 0.005;
+        const double epsKron = 1.0e-3;
         const double norm = 1.0/sqrt(2.0);
         ChebyshevOrthoPoly1st cheb;
 
         unsigned deg = 0;
         CHECK_CLOSE(1.0, cheb.poly(deg, 0.3), epsVal);
 
         deg = 3;
         CHECK_CLOSE(-0.41690962099125364431L/norm, cheb.poly(deg, 1/7.0L), epsVal);
 
         deg = 7;
         CHECK_CLOSE(-0.59498215518301758629L/norm, cheb.poly(deg, 1/11.0L), epsVal);
 
-        GaussLegendreQuadrature glq(1024);
+        FejerQuadrature fej(1000);
         for (unsigned i=0; i<5; ++i)
             for (unsigned j=0; j<5; ++j)
-                CHECK_CLOSE((i == j ? 1.0 : 0.0), cheb.empiricalKroneckerDelta(glq, i, j), epsKron);
+                CHECK_CLOSE((i == j ? 1.0 : 0.0), cheb.empiricalKroneckerDelta(fej, i, j), epsKron);
     }
 }