Review request -- Mathieu fcns

include/xsf/mathieu/README.md

@@ -0,0 +1,39 @@
+This is an implementation of the Mathieu fcns in C/C++.  The
+implementation follows the prototype algos created in Matlab and
+maintained on GitHub at
+https://github.com/brorson/MathieuFcnsFourier.  This impl is a
+header-only library for compatability with Scipy's xsf library.
+
+The following Mathieu fcns are implemented:
+
+*  Angular fcn ce(n,q,v)
+*  Angular fcn se(n,q,v)
+*  Radial (modified) fcn of first kind mc1(n,q,u)
+*  Radial (modified) fcn of first kind ms1(n,q,u)
+*  Radial (modified) fcn of second kind mc2(n,q,u)
+*  Radial (modified) fcn of second kind ms2(n,q,u)
+
+Here, n = fcn order, q = frequency (geometry) parmeter, v = angular
+coord (radians), u = radial coord (au).
+
+I also provide the following utility fcns:
+
+*  Eigenvalue a_n(q)
+*  Eigenvalue b_n(q)
+*  Fourier coeffs A_n^k(q) for ce fcns
+*  Fourier coeffs B_n^k(q) for se fcns
+
+The goal is to provide a replacement of the Mathieu fcn suite used by
+Scipy. 
+
+These programs may be built the usual way on a Linux system using the
+usual GNU build tools.  The main() function runs some simple sanity
+checks on the functions.  In particular, it verifies some output
+values against those computed by the Matlab programs.  I did a lot of
+verification and accuracy testing on the Matlab implementations.
+Therefore, tests run here just make sure the C implementation's
+outputs match those from Matlab.  The code in main() also shows how to
+invoke the various fcns.
+
+Summer 2025, SDB
+

include/xsf/mathieu/besseljyd.h

@@ -0,0 +1,84 @@
+#ifndef BESSELJYD_H
+#define BESSELJYD_H
+
+#include "../bessel.h"
+#include "../config.h"
+
+/*
+ *
+ * This is part of the Mathieu function suite -- a reimplementation
+ * of the Mathieu functions for Scipy.  This file holds helpers
+ * to the Bessel J and Y functions and also returns derivatives
+ * of those fcns.
+ *
+ */
+
+namespace xsf {
+namespace mathieu {
+
+    //==================================================================
+    double besselj(int k, double z) {
+        // This is just a thin wrapper around the Bessel impl in the
+        // std library.
+        double v = (double)k;
+        return xsf::cyl_bessel_j(v, z);
+    }
+
+    //==================================================================
+    double bessely(int k, double z) {
+        // This is just a thin wrapper around the Bessel impl in the
+        // std library.
+        double v = (double)k;
+        return xsf::cyl_bessel_y(v, z);
+    }
+
+    //==================================================================
+    double besseljd(int k, double z) {
+        // This returns the derivative of besselj.  The deriv is
+        // computed using common identities.
+        double y;
+
+        if (k == 0) {
+            double v = 1.0;
+            y = -besselj(v, z);
+        } else {
+            double kp1 = (double)(k + 1);
+            double km1 = (double)(k - 1);
+            y = (besselj(km1, z) - besselj(kp1, z)) / 2.0;
+        }
+
+        // Must flip sign for negative k and odd k.
+        if (k < 0 && ((k % 2) != 0)) {
+            y = -y;
+        }
+
+        return y;
+    }
+
+    //==================================================================
+    double besselyd(int k, double z) {
+        // This returns the derivative of besselj.  The deriv is
+        // computed using common identities.
+        double y;
+
+        if (k == 0) {
+            double v = 1.0;
+            y = -bessely(v, z);
+        } else {
+            double kp1 = (double)(k + 1);
+            double km1 = (double)(k - 1);
+            y = (bessely(km1, z) - bessely(kp1, z)) / 2.0;
+        }
+
+        // Must flip sign for negative k and odd k.
+        if (k < 0 && ((k % 2) != 0)) {
+            y = -y;
+        }
+
+        return y;
+    }
+
+} // namespace mathieu
+} // namespace xsf
+
+#endif // #ifndef BESSELJYD_H

include/xsf/mathieu/make_matrix.h

@@ -0,0 +1,198 @@
+#ifndef MAKE_MATRIX_H
+#define MAKE_MATRIX_H
+
+#include "../config.h"
+#include "../error.h"
+#include "matrix_utils.h"
+
+/*
+ *
+ * This is part of the Mathieu function suite -- a reimplementation
+ * of the Mathieu functions for Scipy.  This file holds the functions
+ * which make the recursion matrices.
+ *
+ * Stuart Brorson, Summer 2025.
+ *
+ */
+
+#define SQRT2 1.414213562373095
+
+namespace xsf {
+namespace mathieu {
+
+    /*-----------------------------------------------
+    This creates the recurrence relation matrix for
+    the even-even Mathieu fcns (ce_2n).
+    Inputs:
+    N = matrix size (related to max order desired).
+    q = shape parameter.
+    Output:
+    A = recurrence matrix (must be calloc'ed in caller).
+    Return:
+    return code = 0 if OK.
+    -------------------------------------------------*/
+    int make_matrix_ee(int N, double q, double *A) {
+        int j;
+        int i;
+
+        // Symmetrize matrix here, then fix in caller.
+        i = MATRIX_IDX(N, 0, 1);
+        A[i] = SQRT2 * q;
+        i = MATRIX_IDX(N, 1, 0);
+        A[i] = SQRT2 * q;
+        i = MATRIX_IDX(N, 1, 1);
+        A[i] = 4.0;
+        i = MATRIX_IDX(N, 1, 2);
+        A[i] = q;
+
+        for (j = 2; j <= N - 2; j++) {
+            i = MATRIX_IDX(N, j, j - 1);
+            A[i] = q;
+            i = MATRIX_IDX(N, j, j);
+            A[i] = (2.0 * j) * (2.0 * j);
+            i = MATRIX_IDX(N, j, j + 1);
+            A[i] = q;
+        }
+
+        i = MATRIX_IDX(N, N - 1, N - 2);
+        A[i] = q;
+        i = MATRIX_IDX(N, N - 1, N - 1);
+        A[i] = (2.0 * (N - 1)) * (2.0 * (N - 1));
+
+        return 0;
+    }
+
+    /*-----------------------------------------------
+    This creates the recurrence relation matrix for
+    the even-odd Mathieu fcns (ce_2n+1).
+    Inputs:
+    N = matrix size (related to max order desired).
+    q = shape parameter.
+    Output:
+    A = recurrence matrix (calloc in caller).
+    Return:
+    return code = 0 if OK.
+    -------------------------------------------------*/
+    int make_matrix_eo(int N, double q, double *A) {
+        int j;
+        int i;
+
+        i = MATRIX_IDX(N, 0, 0);
+        A[i] = 1.0 + q;
+        i = MATRIX_IDX(N, 0, 1);
+        A[i] = q;
+        i = MATRIX_IDX(N, 1, 0);
+        A[i] = q;
+        i = MATRIX_IDX(N, 1, 1);
+        A[i] = 9.0;
+        i = MATRIX_IDX(N, 1, 2);
+        A[i] = q;
+
+        for (j = 2; j <= N - 2; j++) {
+            i = MATRIX_IDX(N, j, j - 1);
+            A[i] = q;
+            i = MATRIX_IDX(N, j, j);
+            A[i] = (2.0 * j + 1.0) * (2.0 * j + 1.0);
+            i = MATRIX_IDX(N, j, j + 1);
+            A[i] = q;
+        }
+
+        i = MATRIX_IDX(N, N - 1, N - 2);
+        A[i] = q;
+        i = MATRIX_IDX(N, N - 1, N - 1);
+        A[i] = (2.0 * (N - 1) + 1.0) * (2.0 * (N - 1) + 1.0);
+
+        return 0;
+    }
+
+    /*-----------------------------------------------
+    This creates the recurrence relation matrix for
+    the odd-even Mathieu fcns (se_2n) -- sometimes called
+    se_2n+2.
+    Inputs:
+    N = matrix size (related to max order desired).
+    q = shape parameter.
+    Output:
+    A = recurrence matrix (calloc in caller).
+    Return:
+    return code = 0 if OK.
+    -------------------------------------------------*/
+    int make_matrix_oe(int N, double q, double *A) {
+        int j;
+        int i;
+
+        i = MATRIX_IDX(N, 0, 0);
+        A[i] = 4.0;
+        i = MATRIX_IDX(N, 0, 1);
+        A[i] = q;
+        i = MATRIX_IDX(N, 1, 0);
+        A[i] = q;
+        i = MATRIX_IDX(N, 1, 1);
+        A[i] = 16.0;
+        i = MATRIX_IDX(N, 1, 2);
+        A[i] = q;
+
+        for (j = 2; j <= N - 2; j++) {
+            i = MATRIX_IDX(N, j, j - 1);
+            A[i] = q;
+            i = MATRIX_IDX(N, j, j);
+            A[i] = (2.0 * (j + 1)) * (2.0 * (j + 1));
+            i = MATRIX_IDX(N, j, j + 1);
+            A[i] = q;
+        }
+
+        i = MATRIX_IDX(N, N - 1, N - 2);
+        A[i] = q;
+        i = MATRIX_IDX(N, N - 1, N - 1);
+        A[i] = (2.0 * N) * (2.0 * N);
+
+        return 0;
+    }
+
+    /*-----------------------------------------------
+    This creates the recurrence relation matrix for
+    the odd-odd Mathieu fcns (se_2n+1).
+    Inputs:
+    N = matrix size (related to max order desired).
+    q = shape parameter.
+    Output:
+    A = recurrence matrix (calloc in caller).
+    Return:
+    return code = 0 if OK.
+    -------------------------------------------------*/
+    int make_matrix_oo(int N, double q, double *A) {
+        int j;
+        int i;
+
+        i = MATRIX_IDX(N, 0, 0);
+        A[i] = 1.0 - q;
+        i = MATRIX_IDX(N, 0, 1);
+        A[i] = q;
+        i = MATRIX_IDX(N, 1, 0);
+        A[i] = q;
+        i = MATRIX_IDX(N, 1, 1);
+        A[i] = 9.0;
+        i = MATRIX_IDX(N, 1, 2);
+        A[i] = q;
+
+        for (j = 2; j <= N - 2; j++) {
+            i = MATRIX_IDX(N, j, j - 1);
+            A[i] = q;
+            i = MATRIX_IDX(N, j, j);
+            A[i] = (2.0 * j + 1.0) * (2.0 * j + 1.0);
+            i = MATRIX_IDX(N, j, j + 1);
+            A[i] = q;
+        }
+
+        i = MATRIX_IDX(N, N - 1, N - 2);
+        A[i] = q;
+        i = MATRIX_IDX(N, N - 1, N - 1);
+        A[i] = (2.0 * N - 1.0) * (2.0 * N - 1.0);
+
+        return 0;
+    }
+
+} // namespace mathieu
+} // namespace xsf
+
+#endif // #ifndef MAKE_MATRIX_H