Restore the order-39 Ahrens-Beylkin rule: it is 572 points, not 552

include/integratorxx/quadratures/s2/ahrens_beylkin.hpp

@@ -3,6 +3,7 @@
 #include <integratorxx/quadrature.hpp>
 #include <integratorxx/quadratures/s2/ahrens_beylkin/ahrens_beylkin_grids.hpp>
 #include <integratorxx/util/copy_grid.hpp>
+#include <integratorxx/util/unsupported_grid.hpp>
 #include <vector>
 
 namespace IntegratorXX {
@@ -64,8 +65,14 @@ template <typename RealType> struct quadrature_traits<AhrensBeylkin<RealType>> {
       detail::copy_grid<ahrens_beylkin_432<RealType>>(points, weights);
     else if (npts == 492)
       detail::copy_grid<ahrens_beylkin_492<RealType>>(points, weights);
-    else if (npts == 552)
-      detail::copy_grid<ahrens_beylkin_552<RealType>>(points, weights);
+    // NOTE: the order-39 rule is 572 points, not 552. The 552-point table --
+    // and the qsph1-39-552DP.dat it came from -- is missing its 20-point
+    // (dodecahedral vertex) orbit, which is why its weights sum to 0.9632 * 4pi
+    // and no weighting of those points integrates to order 39 at all. The
+    // orbit restored, the same points and weights form a valid rule. 552
+    // remains unavailable; ask for 572.
+    else if (npts == 572)
+      detail::copy_grid<ahrens_beylkin_572<RealType>>(points, weights);
     else if (npts == 612)
       detail::copy_grid<ahrens_beylkin_612<RealType>>(points, weights);
     else if (npts == 672)
@@ -162,6 +169,8 @@ template <typename RealType> struct quadrature_traits<AhrensBeylkin<RealType>> {
       detail::copy_grid<ahrens_beylkin_7512<RealType>>(points, weights);
     else if (npts == 15012)
       detail::copy_grid<ahrens_beylkin_15012<RealType>>(points, weights);
+    else
+      detail::throw_unsupported_grid_size("AhrensBeylkin", npts);
 
     return std::make_tuple(points, weights);
   }
@@ -183,7 +192,7 @@ inline static int64_t npts_by_algebraic_order(int64_t order) {
   case 37:
     return 492;
   case 39:
-    return 552;
+    return 572;
   case 41:
     return 612;
   case 44:
@@ -301,7 +310,7 @@ inline static int64_t algebraic_order_by_npts(int64_t npts) {
     return 35;
   case 492:
     return 37;
-  case 552:
+  case 572:
     return 39;
   case 612:
     return 41;

include/integratorxx/quadratures/s2/ahrens_beylkin/ahrens_beylkin_grids.hpp

@@ -10,6 +10,7 @@
 #include "ahrens_beylkin_432.hpp"
 #include "ahrens_beylkin_492.hpp"
 #include "ahrens_beylkin_552.hpp"
+#include "ahrens_beylkin_572.hpp"
 #include "ahrens_beylkin_612.hpp"
 #include "ahrens_beylkin_672.hpp"
 #include "ahrens_beylkin_732.hpp"

include/integratorxx/quadratures/s2/delley.hpp

@@ -2,6 +2,7 @@
 
 #include <integratorxx/quadrature.hpp>
 #include <integratorxx/util/copy_grid.hpp>
+#include <integratorxx/util/unsupported_grid.hpp>
 #include <integratorxx/quadratures/s2/delley/delley_grids.hpp>
 #include <vector>
 
@@ -85,6 +86,8 @@ struct quadrature_traits<Delley<RealType>> {
       detail::copy_grid<delley_3074<RealType>>(points, weights);
     else if(npts == 3470)
       detail::copy_grid<delley_3470<RealType>>(points, weights);
+    else
+      detail::throw_unsupported_grid_size("Delley", npts);
 
     // Pretabulated weights are missing 4 pi
     for(auto i=0; i < npts; i++)

include/integratorxx/quadratures/s2/lebedev_laikov.hpp

@@ -3,6 +3,7 @@
 #include <integratorxx/quadrature.hpp>
 #include <integratorxx/quadratures/s2/lebedev_laikov/lebedev_laikov_grids.hpp>
 #include <integratorxx/util/copy_grid.hpp>
+#include <integratorxx/util/unsupported_grid.hpp>
 #include <vector>
 
 namespace IntegratorXX {
@@ -105,8 +106,12 @@ struct quadrature_traits< LebedevLaikov<RealType> > {
       detail::copy_grid<lebedev_laikov_4334<RealType>>( points, weights );
     else if( npts == 4802 ) 
       detail::copy_grid<lebedev_laikov_4802<RealType>>( points, weights );
+    else if( npts == 5294 ) 
+      detail::copy_grid<lebedev_laikov_5294<RealType>>( points, weights );
     else if( npts == 5810 ) 
       detail::copy_grid<lebedev_laikov_5810<RealType>>( points, weights );
+    else
+      detail::throw_unsupported_grid_size("LebedevLaikov", npts);
 
     // Pretabulated weights are missing 4 pi
     for(size_t i=0; i < npts; i++)
@@ -149,6 +154,7 @@ struct quadrature_traits< LebedevLaikov<RealType> > {
       case 107: return 3890 ; 
       case 113: return 4334 ; 
       case 119: return 4802 ; 
+      case 125: return 5294 ; 
       case 131: return 5810 ; 
       default:  return -1;
     }
@@ -188,6 +194,7 @@ struct quadrature_traits< LebedevLaikov<RealType> > {
       case 3890: return 107 ; 
       case 4334: return 113 ; 
       case 4802: return 119 ; 
+      case 5294: return 125 ; 
       case 5810: return 131 ; 
       default:  return -1;
     }
@@ -227,6 +234,7 @@ struct quadrature_traits< LebedevLaikov<RealType> > {
      else if( order <= 107) return 107; 
      else if( order <= 113) return 113; 
      else if( order <= 119) return 119; 
+     else if( order <= 125) return 125; 
      else                   return 131; 
 
   }

include/integratorxx/quadratures/s2/womersley.hpp

@@ -2,6 +2,7 @@
 
 #include <integratorxx/quadrature.hpp>
 #include <integratorxx/util/copy_grid.hpp>
+#include <integratorxx/util/unsupported_grid.hpp>
 #include <integratorxx/util/create_array.hpp>
 #include <integratorxx/quadratures/s2/womersley/womersley_grids.hpp>
 #include <vector>
@@ -184,6 +185,8 @@ struct quadrature_traits<Womersley<RealType>> {
       detail::copy_grid<womersley_7814<RealType>>(points, weights);
     else if(npts == 7939)
       detail::copy_grid<womersley_7939<RealType>>(points, weights);
+    else
+      detail::throw_unsupported_grid_size("Womersley", npts);
   }
 
   inline static std::tuple<point_container, weight_container> generate(
@@ -319,6 +322,8 @@ struct quadrature_traits<Womersley<RealType>> {
       detail::copy_grid<womersley_1923<RealType>>(points, weights);
     else if(npts == 1986)
       detail::copy_grid<womersley_1986<RealType>>(points, weights);
+    else
+      detail::throw_unsupported_grid_size("Womersley", npts);
     return std::make_tuple(points, weights);
   }
 

include/integratorxx/util/unsupported_grid.hpp

@@ -0,0 +1,33 @@
+#pragma once
+
+#include <cstddef>
+#include <stdexcept>
+#include <string>
+
+namespace IntegratorXX {
+namespace detail {
+
+/**
+ *  @brief Report a request for an angular grid size that is not tabulated.
+ *
+ *  The angular quadratures are only defined for the specific point counts
+ *  that integrate spherical harmonics exactly up to a given algebraic
+ *  order. Requesting any other size is a programming error rather than a
+ *  recoverable condition, but it is diagnosed at runtime because the size
+ *  is typically read from user input.
+ *
+ *  @param[in] family Name of the angular quadrature family
+ *  @param[in] npts   Unsupported number of points that was requested
+ */
+[[noreturn]] inline void throw_unsupported_grid_size(const char* family,
+                                                     size_t npts) {
+  throw std::runtime_error(
+      std::string("IntegratorXX: ") + family + " does not tabulate a " +
+      std::to_string(npts) +
+      "-point grid. Supported sizes are given by "
+      "quadrature_traits<" +
+      family + ">::npts_by_algebraic_order(order).");
+}
+
+}  // namespace detail
+}  // namespace IntegratorXX

test/spherical_generator.cxx

@@ -1,5 +1,7 @@
 #include "catch2/catch_all.hpp"
 #include <iostream>
+#include <algorithm>
+#include <cmath>
 
 #include <integratorxx/quadratures/radial.hpp>
 #include <integratorxx/quadratures/s2.hpp>
@@ -290,6 +292,57 @@ TEMPLATE_LIST_TEST_CASE("S2 Generator", "[sph-gen]", s2_test_types) {
 
 }
 
+/// Exact value of \int_{S^2} x^{2a} y^{2b} z^{2c} d\Omega
+static double exact_even_monomial(int a, int b, int c) {
+  auto dfact = [](int n) { double r = 1.0; for(int k = n; k > 1; k -= 2) r *= k; return r; };
+  return 4.0 * M_PI * dfact(2*a-1) * dfact(2*b-1) * dfact(2*c-1) /
+         dfact(2*(a+b+c) + 1);
+}
+
+TEMPLATE_LIST_TEST_CASE("S2 Grid Validity", "[sph-gen]", s2_test_types) {
+  using namespace IntegratorXX;
+  using angular_type   = TestType;
+  using angular_traits = quadrature_traits<angular_type>;
+
+  // Sizes that no angular scheme tabulates must be diagnosed rather than
+  // silently yielding a zero grid.
+  REQUIRE_THROWS_AS( angular_type(1), std::runtime_error );
+  REQUIRE_THROWS_AS( angular_type(7), std::runtime_error );
+
+  // Every size advertised by the algebraic order tables must be dispatched
+  // and must integrate polynomials exactly up to its claimed order. A
+  // missing dispatch branch would return zero weights; a corrupt table
+  // would integrate low-order monomials incorrectly.
+  for( int64_t order = 0; order <= 200; ++order ) {
+    const auto npts = angular_traits::npts_by_algebraic_order(order);
+    if( npts < 0 ) continue;
+
+    INFO( "algebraic order " << order << ", " << npts << " points" );
+    angular_type aq(npts);
+    REQUIRE( aq.npts() == static_cast<size_t>(npts) );
+
+    const auto& p = aq.points();
+    const auto& w = aq.weights();
+
+    // Total degree 2*(a+b+c) must not exceed the claimed algebraic order.
+    // Capped for run time; degree 8 already exercises every table.
+    const int dmax = static_cast<int>(std::min<int64_t>(order / 2, 4));
+    for( int a = 0;         a <= dmax; ++a )
+    for( int b = 0;     a + b <= dmax; ++b )
+    for( int c = 0; a + b + c <= dmax; ++c ) {
+      double integral = 0.0;
+      for( size_t i = 0; i < aq.npts(); ++i )
+        integral += w[i] * std::pow(p[i][0], 2*a)
+                         * std::pow(p[i][1], 2*b)
+                         * std::pow(p[i][2], 2*c);
+
+      INFO( "monomial x^" << 2*a << " y^" << 2*b << " z^" << 2*c );
+      REQUIRE_THAT( integral,
+        Catch::Matchers::WithinRel(exact_even_monomial(a,b,c), 1e-12) );
+    }
+  }
+}
+
 using sph_test_types = std::tuple<
   std::tuple<bk_type, ah_type>,
   std::tuple<bk_type, de_type>,