Review request -- Mathieu fcns

include/xsf/mathieu/LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2025 Stuart Brorson
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

include/xsf/mathieu/Makefile

@@ -0,0 +1,36 @@
+CC := g++
+# CC := gcc
+CFLAGS := -O3 -std=c++17 -fno-fast-math -march=x86-64 -m128bit-long-double
+#CFLAGS := -O3 -std=c++17 -fno-fast-math -mfpmath=sse -march=x86-64 -m128bit-long-double
+
+SRCS := $(wildcard *.c) 
+OBJS := $(SRCS:.c=.o)
+EXES := main
+TARGETS = main
+INCLUDES := -I. -I/usr/include -I/usr/local/include
+LDFLAGS := -L/usr/lib/x86_64-linux-gnu/lapack/
+#LDFLAGS := -L/usr/lib/x86_64-linux-gnu/openblas-pthread/ -L/usr/lib -L/usr/local/lib 
+#INCLUDES := -I. -I/usr/include/x86_64-linux-gnu/
+#LDFLAGS := -L/usr/lib/x86_64-linux-gnu
+LIBS := -llapacke -llapack -lblas -lm 
+
+
+#=================================================
+all: $(TARGETS)
+
+# Link all object files into the executable
+$(TARGETS): $(OBJS)
+	$(CC) $(OBJS) -o $@ $(CFLAGS) $(LDFLAGS) $(LIBS) 
+
+# Compile each .c into a .o
+%.o: %.c %.h
+	$(CC) $(CFLAGS) $(INCLUDES) -c $< -o $@
+
+# Also handle .c files without headers
+%.o: %.c
+	$(CC) $(CFLAGS) $(INCLUDES) -c $< -o $@
+
+#-------------------------------
+# Clean up directory -- remove executables and intermediate files.
+clean:
+	-rm -f *~ *.o *.obj *.out *.map *.h.gch $(TARGETS) $(OBJS)

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,85 @@
+#ifndef BESSELJYD_H
+#define BESSELJYD_H
+
+#include "../config.h"
+#include "../bessel.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 xsf
+} // namespace mathieu
+
+#endif  // #ifndef BESSELJYD_H