From c27ad1c07f28fc39b632144666789331bfc61e51 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 23 Nov 2021 16:04:30 +0100
Subject: Swiched to using Jacobi polynomials instead of Legendre.

---
 stokes.hpp | 37 ++++++++++++++++++++++++++++++++++---
 1 file changed, 34 insertions(+), 3 deletions(-)

diff --git a/stokes.hpp b/stokes.hpp
index d86b479..76b950f 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -5,6 +5,38 @@
 
 #include <iostream>
 
+template <class Real>
+Real jacobi(unsigned α, unsigned β, unsigned n, Real x) {
+  Real f = 0;
+
+  Real y = (x - 1) / 2;
+  Real z = (x + 1) / 2;
+
+  for (unsigned s = 0; s <= n; s++) {
+    f += pow(y, n - s) * pow(z, s) / (tgamma(s + 1) * tgamma(n + α - s + 1) * tgamma(β + s + 1) * tgamma(n - s + 1));
+  }
+
+  return tgamma(n + α + 1) * tgamma(n + β + 1) * f;
+}
+
+template <class Real>
+Real dJacobi(unsigned α, unsigned β, unsigned n, Real x) {
+  Real f = 0;
+
+  Real y = (x - 1) / 2;
+  Real z = (x + 1) / 2;
+
+  for (unsigned s = 0; s < n; s++) {
+    f += pow(y, n - s - 1) * pow(z, s) / (tgamma(s + 1) * tgamma(n + α - s + 1) * tgamma(β + s + 1) * tgamma(n - s - 1 + 1));
+  }
+
+  for (unsigned s = 1; s <= n; s++) {
+    f += pow(y, n - s) * pow(z, s - 1) / (tgamma(s - 1 + 1) * tgamma(n + α - s + 1) * tgamma(β + s + 1) * tgamma(n - s + 1));
+  }
+
+  return 0.5 * tgamma(n + α + 1) * tgamma(n + β + 1) * f;
+}
+
 template <class Scalar, int ...ps>
 class Cord {
 public:
@@ -28,12 +60,11 @@ public:
   }
 
   Real gCoeff(unsigned i, Real t) const {
-    return (1 - t) * t * std::legendre(i, 2 * t - 1);
+    return (1 - t) * t * jacobi(1, 1, i, 2 * t - 1);
   }
 
   Real dgCoeff(unsigned i, Real t) const {
-    Real x = 2 * t - 1;
-    return ((i - 1.0) * x * std::legendre(i, x) - (i + 1) * std::legendre(i + 1, x)) / 2.0;
+    return (1 - 2 * t) * jacobi(1, 1, i, 2 * t - 1) + 2 * (1 - t) * t * dJacobi(1, 1, i, 2 * t - 1);
   }
 
   Vector<Scalar> f(Real t) const {
-- 
cgit v1.2.3-54-g00ecf