diff --git a/include/xsf/mathieu.h b/include/xsf/mathieu.h index 2431dde44..c5ce88155 100644 --- a/include/xsf/mathieu.h +++ b/include/xsf/mathieu.h @@ -1,228 +1,352 @@ -#pragma once - -#include "error.h" -#include "specfun/specfun.h" - -namespace xsf { - -template -T sem_cva(T m, T q); - -template -void sem(T m, T q, T x, T &csf, T &csd); - -/* Mathieu functions */ -/* Characteristic values */ -template -T cem_cva(T m, T q) { - int int_m, kd = 1; - - if ((m < 0) || (m != floor(m))) { - set_error("mathieu_a", SF_ERROR_DOMAIN, NULL); - return std::numeric_limits::quiet_NaN(); - } - int_m = (int)m; - if (q < 0) { - /* https://dlmf.nist.gov/28.2#E26 */ - if (int_m % 2 == 0) { - return cem_cva(m, -q); - } else { - return sem_cva(m, -q); - } - } - - if (int_m % 2) { - kd = 2; - } - return specfun::cva2(kd, int_m, q); -} - -template -T sem_cva(T m, T q) { - int int_m, kd = 4; - - if ((m <= 0) || (m != floor(m))) { - set_error("mathieu_b", SF_ERROR_DOMAIN, NULL); - return std::numeric_limits::quiet_NaN(); - } - int_m = (int)m; - if (q < 0) { - /* https://dlmf.nist.gov/28.2#E26 */ - if (int_m % 2 == 0) { - return sem_cva(m, -q); - } else { - return cem_cva(m, -q); - } - } - if (int_m % 2) { - kd = 3; - } - return specfun::cva2(kd, int_m, q); -} - -/* Mathieu functions */ -template -void cem(T m, T q, T x, T &csf, T &csd) { - int int_m, kf = 1, sgn; - T f = 0.0, d = 0.0; - if ((m < 0) || (m != floor(m))) { - csf = std::numeric_limits::quiet_NaN(); - csd = std::numeric_limits::quiet_NaN(); - set_error("mathieu_cem", SF_ERROR_DOMAIN, NULL); - } else { - int_m = (int)m; - if (q < 0) { - /* https://dlmf.nist.gov/28.2#E34 */ - if (int_m % 2 == 0) { - sgn = ((int_m / 2) % 2 == 0) ? 1 : -1; - cem(m, -q, 90 - x, f, d); - csf = sgn * f; - csd = -sgn * d; - - } else { - sgn = ((int_m / 2) % 2 == 0) ? 1 : -1; - sem(m, -q, 90 - x, f, d); - csf = sgn * f; - csd = -sgn * d; - } - } else { - using specfun::Status; - Status status = specfun::mtu0(kf, int_m, q, x, &csf, &csd); - if (status != Status::OK) { - csf = std::numeric_limits::quiet_NaN(); - csd = std::numeric_limits::quiet_NaN(); - sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; - set_error("mathieu_cem", sf_error, NULL); - } - } - } -} - -template -void sem(T m, T q, T x, T &csf, T &csd) { - int int_m, kf = 2, sgn; - T f = 0.0, d = 0.0; - if ((m < 0) || (m != floor(m))) { - csf = std::numeric_limits::quiet_NaN(); - csd = std::numeric_limits::quiet_NaN(); - set_error("mathieu_sem", SF_ERROR_DOMAIN, NULL); - } else { - int_m = (int)m; - if (int_m == 0) { - csf = 0; - csd = 0; - } else if (q < 0) { - /* https://dlmf.nist.gov/28.2#E34 */ - if (int_m % 2 == 0) { - sgn = ((int_m / 2) % 2 == 0) ? -1 : 1; - sem(m, -q, 90 - x, f, d); - csf = sgn * f; - csd = -sgn * d; - } else { - sgn = ((int_m / 2) % 2 == 0) ? 1 : -1; - cem(m, -q, 90 - x, f, d); - csf = sgn * f; - csd = -sgn * d; - } - } else { - using specfun::Status; - Status status = specfun::mtu0(kf, int_m, q, x, &csf, &csd); - if (status != Status::OK) { - csf = std::numeric_limits::quiet_NaN(); - csd = std::numeric_limits::quiet_NaN(); - sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; - set_error("mathieu_sem", sf_error, NULL); - } - } - } -} - -template -void mcm1(T m, T q, T x, T &f1r, T &d1r) { - int int_m, kf = 1, kc = 1; - T f2r = 0.0, d2r = 0.0; - - if ((m < 0) || (m != floor(m)) || (q < 0)) { - f1r = std::numeric_limits::quiet_NaN(); - d1r = std::numeric_limits::quiet_NaN(); - set_error("mathieu_modcem1", SF_ERROR_DOMAIN, NULL); - } else { - using specfun::Status; - int_m = (int)m; - Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r); - if (status != Status::OK) { - f1r = std::numeric_limits::quiet_NaN(); - d1r = std::numeric_limits::quiet_NaN(); - sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; - set_error("mathieu_modcem1", sf_error, NULL); - } - } -} - -template -void msm1(T m, T q, T x, T &f1r, T &d1r) { - int int_m, kf = 2, kc = 1; - T f2r = 0.0, d2r = 0.0; - - if ((m < 1) || (m != floor(m)) || (q < 0)) { - f1r = std::numeric_limits::quiet_NaN(); - d1r = std::numeric_limits::quiet_NaN(); - set_error("mathieu_modsem1", SF_ERROR_DOMAIN, NULL); - } else { - using specfun::Status; - int_m = (int)m; - Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r); - if (status != Status::OK) { - f1r = std::numeric_limits::quiet_NaN(); - d1r = std::numeric_limits::quiet_NaN(); - sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; - set_error("mathieu_modsem1", sf_error, NULL); - } - } -} - -template -void mcm2(T m, T q, T x, T &f2r, T &d2r) { - int int_m, kf = 1, kc = 2; - T f1r = 0.0, d1r = 0.0; - - if ((m < 0) || (m != floor(m)) || (q < 0)) { - f2r = std::numeric_limits::quiet_NaN(); - d2r = std::numeric_limits::quiet_NaN(); - set_error("mathieu_modcem2", SF_ERROR_DOMAIN, NULL); - } else { - using specfun::Status; - int_m = (int)m; - Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r); - if (status != Status::OK) { - f2r = std::numeric_limits::quiet_NaN(); - d2r = std::numeric_limits::quiet_NaN(); - sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; - set_error("mathieu_modcem2", sf_error, NULL); - } - } -} - -template -void msm2(T m, T q, T x, T &f2r, T &d2r) { - int int_m, kf = 2, kc = 2; - T f1r = 0.0, d1r = 0.0; - - if ((m < 1) || (m != floor(m)) || (q < 0)) { - f2r = std::numeric_limits::quiet_NaN(); - d2r = std::numeric_limits::quiet_NaN(); - set_error("mathieu_modsem2", SF_ERROR_DOMAIN, NULL); - } else { - using specfun::Status; - int_m = (int)m; - Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r); - if (status != Status::OK) { - f2r = std::numeric_limits::quiet_NaN(); - d2r = std::numeric_limits::quiet_NaN(); - sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; - set_error("mathieu_modsem2", sf_error, NULL); - } - } -} - -} // namespace xsf +#pragma once + +#include "error.h" +#include "mathieu/besseljyd.h" +#include "mathieu/make_matrix.h" +#include "mathieu/mathieu_coeffs.h" +#include "mathieu/mathieu_eigs.h" +#include "mathieu/mathieu_fcns.h" +#include "mathieu/matrix_utils.h" + +/* + * + * This is part of the Mathieu function suite -- a reimplementation + * of the Mathieu functions for Scipy. This file #includes all the + * fcn impls in the mathieu/ subdirectory, and provides translation + * from the Scipy call to the calling signature I implemented in my + * reimplementation. + * + * Stuart Brorson, Summer 2025. + * + */ + +namespace xsf { + +/* Characteristic values */ +//------------------------------------------------------------- +/** + * Mathieu characteristic values (eigenvalues) for even parity functions. + * + * Even parity characteristic values a. + * + * @param m Eigenvalue order. Must be positive integer less than 500. + * @param q Mathieu parameter q. Real number. + * @return Mathieu eigenvalue a. + */ +template +T cem_cva(T m, T q) { + // This returns the even Mathieu characteristic value (eigenvalue) a. + + // Check for invalid Mathieu order. + if ((m < 0) || (m != floor(m))) { + set_error("mathieu_a", SF_ERROR_DOMAIN, NULL); + return std::numeric_limits::quiet_NaN(); + } + + // Must cast to correct types prior to fcn call. + int im = static_cast(m); + double dq = static_cast(q); + double da; + + int retcode = xsf::mathieu::mathieu_a(im, dq, &da); + if (retcode != SF_ERROR_OK) { + set_error("mathieu_a", SF_ERROR_NO_RESULT, NULL); + return std::numeric_limits::quiet_NaN(); + } + + // Now cast back. + T a = static_cast(da); + return a; +} + +//------------------------------------------------------------- +template +/** + * Mathieu characteristic values (eigenvalues) b for odd functions. + * + * + * Odd parity characteristic values b. + * + * @param m Eigenvalue order. Must be positive integer less than 500. + * @param q Mathieu parameter q. Real number. + * @return Mathieu eigenvalue b. + */ +T sem_cva(T m, T q) { + // This returns the odd Mathieu characteristic value (eigenvalue) b. + + // Check for invalid Mathieu order. + if ((m < 1) || (m != floor(m))) { + set_error("mathieu_b", SF_ERROR_DOMAIN, NULL); + return std::numeric_limits::quiet_NaN(); + } + + // Must cast to correct types prior to fcn call. + int im = static_cast(m); + double dq = static_cast(q); + double db; + + int retcode = xsf::mathieu::mathieu_b(im, dq, &db); + if (retcode != SF_ERROR_OK) { + set_error("mathieu_b", SF_ERROR_NO_RESULT, NULL); + return std::numeric_limits::quiet_NaN(); + } + + // Now cast back. + T b = static_cast(db); + return b; +} + +//--------------------------------------------------------------- +/* Mathieu functions */ +/** + * Even parity Mathieu angular function ce(m, q, x) + * + * This implementation of ce follows the definitions on the + * DLMF, https://dlmf.nist.gov/28 + * + * @param m Function order. Must be positive integer less than 500. + * @param q Parameter q. Real number. + * @param x Angular coordinate x (radians). Real number. + * @param csf Value of function. Real number. + * @param csd Value of derivative w.r.t. x. Real number. + */ +template +void cem(T m, T q, T x, T &csf, T &csd) { + + if ((m < 0) || (m != floor(m))) { + csf = std::numeric_limits::quiet_NaN(); + csd = std::numeric_limits::quiet_NaN(); + set_error("mathieu_ce", SF_ERROR_DOMAIN, NULL); + } else { + + // Must cast to correct types prior to fcn call. + int im = static_cast(m); + double dq = static_cast(q); + double dx = static_cast(x); + double dcsf; + double dcsd; + + // Call fcn and cast back. + int retcode = xsf::mathieu::mathieu_ce(im, dq, dx, &dcsf, &dcsd); + if (retcode != SF_ERROR_OK) { + csf = std::numeric_limits::quiet_NaN(); + csd = std::numeric_limits::quiet_NaN(); + set_error("mathieu_ce", (sf_error_t)retcode, NULL); + } else { + csf = static_cast(dcsf); + csd = static_cast(dcsd); + } + } +} + +//--------------------------------------------------------------- +/** + * Odd parity Mathieu angular function se(m, q, x) + * + * This implementation of ce follows the definitions on the + * DLMF, https://dlmf.nist.gov/28 + * + * @param m Function order. Must be positive integer less than 500. + * @param q Parameter q. Real number. + * @param x Angular coordinate x (radians). Real number. + * @param ssf Value of function. Real number. + * @param ssd Value of derivative w.r.t. x. Real number + */ +template +void sem(T m, T q, T x, T &ssf, T &ssd) { + + if ((m < 1) || (m != floor(m))) { + ssf = std::numeric_limits::quiet_NaN(); + ssd = std::numeric_limits::quiet_NaN(); + set_error("mathieu_sem", SF_ERROR_DOMAIN, NULL); + } else { + + // Must cast to correct types prior to fcn call. + int im = static_cast(m); + double dq = static_cast(q); + double dx = static_cast(x); + double dssf; + double dssd; + + // Call fcn and cast back. + int retcode = xsf::mathieu::mathieu_se(im, dq, dx, &dssf, &dssd); + if (retcode != SF_ERROR_OK) { + ssf = std::numeric_limits::quiet_NaN(); + ssd = std::numeric_limits::quiet_NaN(); + set_error("mathieu_sem", (sf_error_t)retcode, NULL); + } else { + ssf = static_cast(dssf); + ssd = static_cast(dssd); + } + } +} + +//--------------------------------------------------------------- +/** + * Even parity modified (radial) Mathieu function of first kind Mc1(m, q, x) + * + * This implementation of ce follows the definitions on the + * DLMF, https://dlmf.nist.gov/28 + * + * @param m Function order. Must be positive integer less than 500. + * @param q Parameter q. Positive real number + * @param x Radial coordinate x. Positive real number. + * @param f1r Value of function. Real number. + * @param d1r Value of derivative w.r.t. x. Real number + */ +template +void mcm1(T m, T q, T x, T &f1r, T &d1r) { + + if ((m < 0) || (m != floor(m))) { + f1r = std::numeric_limits::quiet_NaN(); + d1r = std::numeric_limits::quiet_NaN(); + set_error("mathieu_mcm1", SF_ERROR_DOMAIN, NULL); + } else { + + // Must cast to correct types prior to fcn call. + int im = static_cast(m); + double dq = static_cast(q); + double dx = static_cast(x); + double df1r; + double dd1r; + + // Call fcn and cast back. + int retcode = xsf::mathieu::mathieu_modmc1(im, dq, dx, &df1r, &dd1r); + if (retcode != SF_ERROR_OK) { + f1r = std::numeric_limits::quiet_NaN(); + d1r = std::numeric_limits::quiet_NaN(); + set_error("mathieu_mcm1", (sf_error_t)retcode, NULL); + } else { + f1r = static_cast(df1r); + d1r = static_cast(dd1r); + } + } +} + +//--------------------------------------------------------------- +/** + * Odd parity modified (radial) Mathieu function of first kind Ms1(m, q, x) + * + * This implementation of ce follows the definitions on the + * DLMF, https://dlmf.nist.gov/28 + * + * @param m Function order. Must be positive integer less than 500. + * @param q Parameter q. Positive real number + * @param x Radial coordinate x. Positive real number. + * @param f1r Value of function. Real number. + * @param d1r Value of derivative w.r.t. x. Real number + */ +template +void msm1(T m, T q, T x, T &f1r, T &d1r) { + + if ((m < 1) || (m != floor(m))) { + f1r = std::numeric_limits::quiet_NaN(); + d1r = std::numeric_limits::quiet_NaN(); + set_error("mathieu_msm1", SF_ERROR_DOMAIN, NULL); + } else { + + // Must cast to correct types prior to fcn call. + int im = static_cast(m); + double dq = static_cast(q); + double dx = static_cast(x); + double df1r; + double dd1r; + + // Call fcn and cast back. + int retcode = xsf::mathieu::mathieu_modms1(im, dq, dx, &df1r, &dd1r); + if (retcode != SF_ERROR_OK) { + f1r = std::numeric_limits::quiet_NaN(); + d1r = std::numeric_limits::quiet_NaN(); + set_error("mathieu_msm1", (sf_error_t)retcode, NULL); + } else { + f1r = static_cast(df1r); + d1r = static_cast(dd1r); + } + } +} + +//--------------------------------------------------------------- +/** + * Even parity modified (radial) Mathieu function of second kind Mc2(m, q, x) + * + * This implementation of ce follows the definitions on the + * DLMF, https://dlmf.nist.gov/28 + * + * @param m Function order. Must be positive integer less than 500. + * @param q Parameter q. Positive real number + * @param x Radial coordinate x. Positive real number. + * @param f2r Value of function. Real number. + * @param d2r Value of derivative w.r.t. x. Real number + */ +template +void mcm2(T m, T q, T x, T &f2r, T &d2r) { + + if ((m < 0) || (m != floor(m))) { + f2r = std::numeric_limits::quiet_NaN(); + d2r = std::numeric_limits::quiet_NaN(); + set_error("mathieu_mcm2", SF_ERROR_DOMAIN, NULL); + } else { + + // Must cast to correct types prior to fcn call. + int im = static_cast(m); + double dq = static_cast(q); + double dx = static_cast(x); + double df2r; + double dd2r; + + // Call fcn and cast back. + int retcode = xsf::mathieu::mathieu_modmc2(im, dq, dx, &df2r, &dd2r); + if (retcode != SF_ERROR_OK) { + f2r = std::numeric_limits::quiet_NaN(); + d2r = std::numeric_limits::quiet_NaN(); + set_error("mathieu_mcm2", (sf_error_t)retcode, NULL); + } else { + f2r = static_cast(df2r); + d2r = static_cast(dd2r); + } + } +} + +//--------------------------------------------------------------- +/** + * Odd parity modified (radial) Mathieu function of second kind Ms2(m, q, x) + * + * This implementation of ce follows the definitions on the + * DLMF, https://dlmf.nist.gov/28 + * + * @param m Function order. Must be positive integer less than 500. + * @param q Parameter q. Positive real number + * @param x Radial coordinate x. Positive real number. + * @param f2r Value of function. Real number. + * @param d2r Value of derivative w.r.t. x. Real number + */ +template +void msm2(T m, T q, T x, T &f2r, T &d2r) { + + if ((m < 1) || (m != floor(m))) { + f2r = std::numeric_limits::quiet_NaN(); + d2r = std::numeric_limits::quiet_NaN(); + set_error("mathieu_msm2", SF_ERROR_DOMAIN, NULL); + } else { + + // Must cast to correct types prior to fcn call. + int im = static_cast(m); + double dq = static_cast(q); + double dx = static_cast(x); + double df2r; + double dd2r; + + // Call fcn and cast back. + int retcode = xsf::mathieu::mathieu_modms2(im, dq, dx, &df2r, &dd2r); + if (retcode != SF_ERROR_OK) { + f2r = std::numeric_limits::quiet_NaN(); + d2r = std::numeric_limits::quiet_NaN(); + set_error("mathieu_msm2", (sf_error_t)retcode, NULL); + } else { + f2r = static_cast(df2r); + d2r = static_cast(dd2r); + } + } +} + +} // namespace xsf diff --git a/include/xsf/mathieu/README.md b/include/xsf/mathieu/README.md new file mode 100644 index 000000000..5e22ce075 --- /dev/null +++ b/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 + diff --git a/include/xsf/mathieu/besseljyd.h b/include/xsf/mathieu/besseljyd.h new file mode 100644 index 000000000..f47c4f44c --- /dev/null +++ b/include/xsf/mathieu/besseljyd.h @@ -0,0 +1,86 @@ +#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. + * + * Stuart Brorson -- Summer and Fall 2025. + * + */ + +namespace xsf { +namespace mathieu { + + //================================================================== + double besselj(int k, double z) { + // This is just a thin wrapper around the Bessel impl in the + // xsf 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 + // xsf 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 diff --git a/include/xsf/mathieu/make_matrix.h b/include/xsf/mathieu/make_matrix.h new file mode 100644 index 000000000..bfd012c4e --- /dev/null +++ b/include/xsf/mathieu/make_matrix.h @@ -0,0 +1,210 @@ +#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 = SF_ERROR_OK 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 SF_ERROR_OK; + } + + /*----------------------------------------------- + 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 = SF_ERROR_OK 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 SF_ERROR_OK; + } + + /*----------------------------------------------- + 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 = SF_ERROR_OK 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 SF_ERROR_OK; + } + + /*----------------------------------------------- + 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 = SF_ERROR_OK 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 SF_ERROR_OK; + } + +} // namespace mathieu +} // namespace xsf + +#endif // #ifndef MAKE_MATRIX_H diff --git a/include/xsf/mathieu/mathieu_coeffs.h b/include/xsf/mathieu/mathieu_coeffs.h new file mode 100644 index 000000000..3450ca0b9 --- /dev/null +++ b/include/xsf/mathieu/mathieu_coeffs.h @@ -0,0 +1,266 @@ +#ifndef MATHIEU_COEFFS_H +#define MATHIEU_COEFFS_H + +#include "../config.h" +#include "../error.h" +#include "make_matrix.h" +#include "matrix_utils.h" +#include + +#define SQRT2 1.414213562373095 + +/* + * + * This is part of the Mathieu function suite -- a reimplementation + * of the Mathieu functions for Scipy. This file holds the functions + * which return the Mathieu A and B coefficients used in the Fourier + * series computing the Mathieu fcns. + * + * Stuart Brorson, Summer 2025. + * + */ + +/* DSYEV_ prototype */ +#ifdef __cplusplus +extern "C" { +#endif +void dsyev_(char *jobz, char *uplo, int *n, double *a, int *lda, double *w, double *work, int *lwork, int *info); +#ifdef __cplusplus +} +#endif + +namespace xsf { +namespace mathieu { + + //------------------------------------------------------ + int mathieu_coeffs_ee(int N, double q, int m, double *AA) { + // Returns Fourier coeffs for the mth order ce_2n Mathieu fcn. + // Allowed value of m = 0, 2, 4, 6, ... + // Inputs: + // N = size of recursion matrix to use. + // q = frequency parameter + // m = order of Mathieu fcn desired. + // Output: + // AA = length N vector preallocated to hold coeffs. + // Returns SF_ERROR_OK if all goes well. + + int retcode = SF_ERROR_OK; + + // Bail out if m is not even. + if (m % 2 != 0) + return SF_ERROR_ARG; + + // Allocate recursion matrix + std::vector A(N * N); + + // Do EVD + retcode = make_matrix_ee(N, q, A.data()); + if (retcode != 0) { + return SF_ERROR_NO_RESULT; + } + + char V[1] = {'V'}; + char U[1] = {'U'}; + double wkopt; + /* Query and allocate the optimal workspace */ + int lwork = -1; + dsyev_(V, U, &N, A.data(), &N, AA, &wkopt, &lwork, &retcode); + lwork = (int)wkopt; + std::vector work(lwork); + + /* Solve eigenproblem */ + dsyev_(V, U, &N, A.data(), &N, AA, work.data(), &lwork, &retcode); + + // Check return code from dsyev and bail if it's not 0. + if (retcode != 0) { + return SF_ERROR_NO_RESULT; + } + + // Sort AA vector from lowest to highest + // quickSort(AA, 0, N-1); + // print_matrix(AA, N, 1); + + // Undo sqrt(2) in make_matrix by normalizing elet in first col by sqrt(2). + int idx; + int row = m / 2; + idx = MATRIX_IDX(N, row, 0); + AA[0] = A[idx] / SQRT2; + // Transfer remaining elets in correct row to coeff vector. + for (int j = 1; j < N; j++) { + idx = MATRIX_IDX(N, row, j); + AA[j] = A[idx]; + } + return retcode; + } + + //------------------------------------------------------ + int mathieu_coeffs_eo(int N, double q, int m, double *AA) { + // Returns Fourier coeffs for the mth order ce_2n+1 Mathieu fcn. + // Allowed value of m = 1, 3, 5, 7 ... + + int retcode = SF_ERROR_OK; + + // Bail out if m is not odd. + if (m % 2 != 1) + return SF_ERROR_ARG; + + // Allocate recursion matrix + std::vector A(N * N); + + // Do EVD + retcode = make_matrix_eo(N, q, A.data()); + if (retcode != 0) { + return SF_ERROR_NO_RESULT; + } + + char V[1] = {'V'}; + char U[1] = {'U'}; + double wkopt; + + /* Query and allocate the optimal workspace */ + int lwork = -1; + dsyev_(V, U, &N, A.data(), &N, AA, &wkopt, &lwork, &retcode); + lwork = (int)wkopt; + std::vector work(lwork); + + /* Solve eigenproblem */ + dsyev_(V, U, &N, A.data(), &N, AA, work.data(), &lwork, &retcode); + + // Check return code from dsyev and bail if it's not 0. + if (retcode != 0) { + return SF_ERROR_NO_RESULT; + } + + // Sort AA vector from lowest to highest + // quickSort(AA, 0, N-1); + // print_matrix(AA, N, 1); + + // Transfer correct row to coeff vector. + int idx; + int row = (m - 1) / 2; + // Transfer elets in correct row to coeff vector. + for (int j = 0; j < N; j++) { + idx = MATRIX_IDX(N, row, j); + AA[j] = A[idx]; + } + return retcode; + } + + //------------------------------------------------------ + int mathieu_coeffs_oe(int N, double q, int m, double *AA) { + // Returns Fourier coeffs for the mth order se_2n Mathieu fcn. + // Allowed value of m = 2, 4, 6, ... + // Inputs: + // N = size of recursion matrix to use. + // q = frequency parameter + // m = order of Mathieu fcn desired. + // Output: + // AA = length N vector preallocated to hold coeffs. + // Returns 0 if all goes well. Must put check on calloc + // here. + + int retcode = SF_ERROR_OK; + + // Bail out if m is not even or >= 2. + if ((m % 2 != 0) || (m < 2)) + return SF_ERROR_ARG; + + // Allocate recursion matrix + std::vector A(N * N); + + // Do EVD + retcode = make_matrix_oe(N, q, A.data()); + if (retcode != 0) { + return SF_ERROR_NO_RESULT; + } + + // Work in local scope. + char V[1] = {'V'}; + char U[1] = {'U'}; + double wkopt; + + /* Query and allocate the optimal workspace */ + int lwork = -1; + dsyev_(V, U, &N, A.data(), &N, AA, &wkopt, &lwork, &retcode); + lwork = (int)wkopt; + std::vector work(lwork); + + /* Solve eigenproblem */ + dsyev_(V, U, &N, A.data(), &N, AA, work.data(), &lwork, &retcode); + + // Bail out if dsyev doesn't return 0. + if (retcode != 0) { + return SF_ERROR_NO_RESULT; + } + + // Sort AA vector from lowest to highest + // quickSort(AA, 0, N-1); + // print_matrix(AA, N, 1); + + // Transfer remaining elets in correct row to coeff vector. + int idx; + int row = (m - 2) / 2; + for (int j = 0; j < N; j++) { + idx = MATRIX_IDX(N, row, j); + AA[j] = A[idx]; + } + return retcode; + } + + //------------------------------------------------------ + int mathieu_coeffs_oo(int N, double q, int m, double *AA) { + // Returns Fourier coeffs for the mth order se_2n+1 Mathieu fcn. + // Allowed value of m = 1, 3, 5, 7 ... + + int retcode = SF_ERROR_OK; + + // Bail out if m is not odd. + if (m % 2 != 1) + return SF_ERROR_ARG; + + // Allocate recursion matrix + std::vector A(N * N); + + // Do EVD + retcode = make_matrix_oo(N, q, A.data()); + if (retcode != 0) { + return SF_ERROR_NO_RESULT; + } + + // Work in local scope + char V[1] = {'V'}; + char U[1] = {'U'}; + double wkopt; + + /* Query and allocate the optimal workspace */ + int lwork = -1; + dsyev_(V, U, &N, A.data(), &N, AA, &wkopt, &lwork, &retcode); + lwork = (int)wkopt; + std::vector work(lwork); + /* Solve eigenproblem */ + dsyev_(V, U, &N, A.data(), &N, AA, work.data(), &lwork, &retcode); + + // Bail out if dsyev didn't return 0; + if (retcode != 0) { + return SF_ERROR_NO_RESULT; + } + + // Sort AA vector from lowest to highest + // quickSort(AA, 0, N-1); + // print_matrix(AA, N, 1); + + // Transfer correct row to coeff vector. + int idx; + int row = (m - 1) / 2; + // Transfer elets in correct row to coeff vector. + for (int j = 0; j < N; j++) { + idx = MATRIX_IDX(N, row, j); + AA[j] = A[idx]; + } + return retcode; + } + +} // namespace mathieu +} // namespace xsf + +#endif // #ifndef MATHIEU_COEFFS_H diff --git a/include/xsf/mathieu/mathieu_eigs.h b/include/xsf/mathieu/mathieu_eigs.h new file mode 100644 index 000000000..d78693e2b --- /dev/null +++ b/include/xsf/mathieu/mathieu_eigs.h @@ -0,0 +1,226 @@ +#ifndef MATHIEU_EIGS_H +#define MATHIEU_EIGS_H + +#include "../config.h" +#include "../error.h" +#include "make_matrix.h" +#include "matrix_utils.h" +#include + +/* + * + * This is part of the Mathieu function suite -- a reimplementation + * of the Mathieu functions for Scipy. This file holds the functions + * which return the Mathieu eigenvalues (characteristic values) a and + * b as a function of parameter q. + * + * Stuart Brorson, summer 2025. + * + */ + +/* DSYEV_ prototype */ +#ifdef __cplusplus +extern "C" { +#endif +void dsyev_(char *jobz, char *uplo, int *n, double *a, int *lda, double *w, double *work, int *lwork, int *info); +#ifdef __cplusplus +} +#endif + +namespace xsf { +namespace mathieu { + + //------------------------------------------------------ + // This is the Mathieu characteristic value (eigenvalue) + // a for even fcns. + int mathieu_a(int m, double q, double *a) { + // printf("--> mathieu_a, m = %d, q = %e\n", m, q); + + int N = m + 25; // Sets size of recursion matrix + int retcode = SF_ERROR_OK; + + if (m > 500) { + // Don't support absurdly larger orders for now. + *a = std::numeric_limits::quiet_NaN(); + return SF_ERROR_DOMAIN; + } + + // Allocate recursion matrix + std::vector A(N * N); + + // Allocate vector for eigenvalues + std::vector ww(N); + + // Do EVD + if (m % 2 == 0) { + // Even order m + retcode = make_matrix_ee(N, q, A.data()); + if (retcode != SF_ERROR_OK) { + *a = std::numeric_limits::quiet_NaN(); + return SF_ERROR_OTHER; // Not sure what went wrong. + } + + char V = 'V'; + char U = 'U'; + double wkopt; + + /* Query and allocate the optimal workspace */ + int lwork = -1; + dsyev_(&V, &U, &N, A.data(), &N, ww.data(), &wkopt, &lwork, &retcode); + lwork = (int)wkopt; + std::vector work(lwork); + + /* Solve eigenproblem */ + dsyev_(&V, &U, &N, A.data(), &N, ww.data(), work.data(), &lwork, &retcode); + + // Check if dsyev was successful + if (retcode != 0) { + *a = std::numeric_limits::quiet_NaN(); + return SF_ERROR_NO_RESULT; + } + + // Sort ww vector from lowest to highest + // quickSort(ww, 0, N-1); + // print_matrix(ww, N, 1); + + // Now figure out which one to return. + int idx = m / 2; + *a = ww[idx]; + + } else { + // Odd order m + retcode = make_matrix_eo(N, q, A.data()); + if (retcode != SF_ERROR_OK) { + *a = std::numeric_limits::quiet_NaN(); + return SF_ERROR_OTHER; + } + + char V = 'V'; + char U = 'U'; + double wkopt; + + /* Query and allocate the optimal workspace */ + int lwork = -1; + dsyev_(&V, &U, &N, A.data(), &N, ww.data(), &wkopt, &lwork, &retcode); + lwork = (int)wkopt; + std::vector work(lwork); + + /* Solve eigenproblem */ + dsyev_(&V, &U, &N, A.data(), &N, ww.data(), work.data(), &lwork, &retcode); + + // Check if dsyev was successful + if (retcode != 0) { + *a = std::numeric_limits::quiet_NaN(); + return SF_ERROR_NO_RESULT; + } + + // Sort ww vector from lowest to highest + // quickSort(ww, 0, N-1); + // print_matrix(ww, N, 1); + + // Now figure out which one to return. + int idx = (m - 1) / 2; + *a = ww[idx]; + } + + // printf("<-- mathieu_a\n"); + return retcode; + } + + //------------------------------------------------------ + int mathieu_b(int m, double q, double *b) { + // This computes the Mathieu characteristic value (eigenvalue) + // for odd fcns. + // printf("--> mathieu_b, m = %d, q = %e\n", m, q); + int N = m + 25; // Sets size of recursion matrix + int retcode = SF_ERROR_OK; + + if (m > 500) { + // Don't support absurdly larger orders for now. + *b = std::numeric_limits::quiet_NaN(); + return SF_ERROR_DOMAIN; + } + + // Allocate recursion matrix + std::vector B(N * N); + + // Allocate vector for eigenvalues + std::vector ww(N); + + // Do EVD + if (m % 2 == 0) { + // Even order m + retcode = make_matrix_oe(N, q, B.data()); + if (retcode != SF_ERROR_OK) { + *b = std::numeric_limits::quiet_NaN(); + return SF_ERROR_OTHER; + } + + char V = 'V'; + char U = 'U'; + double wkopt; + + /* Query and allocate the optimal workspace */ + int lwork = -1; + dsyev_(&V, &U, &N, B.data(), &N, ww.data(), &wkopt, &lwork, &retcode); + lwork = (int)wkopt; + std::vector work(lwork); + + /* Solve eigenproblem */ + dsyev_(&V, &U, &N, B.data(), &N, ww.data(), work.data(), &lwork, &retcode); + + if (retcode != 0) { + *b = std::numeric_limits::quiet_NaN(); + return SF_ERROR_NO_RESULT; + } + + // Sort ww vector from lowest to highest + // quickSort(ww, 0, N-1); + // print_matrix(ww, N, 1); + + // Now figure out which one to return. + int idx = (m - 2) / 2; + *b = ww[idx]; + + } else { + // Odd order m + retcode = make_matrix_oo(N, q, B.data()); + if (retcode != SF_ERROR_OK) { + *b = std::numeric_limits::quiet_NaN(); + return SF_ERROR_OTHER; + } + + char V = 'V'; + char U = 'U'; + double wkopt; + + /* Query and allocate the optimal workspace */ + int lwork = -1; + dsyev_(&V, &U, &N, B.data(), &N, ww.data(), &wkopt, &lwork, &retcode); + lwork = (int)wkopt; + std::vector work(lwork); + /* Solve eigenproblem */ + dsyev_(&V, &U, &N, B.data(), &N, ww.data(), work.data(), &lwork, &retcode); + + if (retcode != 0) { + *b = std::numeric_limits::quiet_NaN(); + return SF_ERROR_NO_RESULT; + } + + // Sort ww vector from lowest to highest + // quickSort(ww, 0, N-1); + // print_matrix(ww, N, 1); + + // Now figure out which one to return. + int idx = (m - 1) / 2; + *b = ww[idx]; + } + + // printf("<-- mathieu_b\n"); + return retcode; + } + +} // namespace mathieu +} // namespace xsf + +#endif // #ifndef MATHIEU_EIGS_H diff --git a/include/xsf/mathieu/mathieu_fcns.h b/include/xsf/mathieu/mathieu_fcns.h new file mode 100644 index 000000000..a6a565460 --- /dev/null +++ b/include/xsf/mathieu/mathieu_fcns.h @@ -0,0 +1,1480 @@ +#ifndef MATHIEU_FCNS_H +#define MATHIEU_FCNS_H + +#include "../config.h" +#include "../error.h" +#include "besseljyd.h" +#include "mathieu_coeffs.h" +#include "matrix_utils.h" +#include +#include + +/* + * + * This is part of the Mathieu function suite -- a reimplementation + * of the Mathieu functions for Scipy. This file holds the function + * implementations themselves. The prototype was written in Matlab + * and validated. This is a translation from Matlab to C++. + * + * Stuart Brorson, Summer 2025. + * + */ + +namespace xsf { +namespace mathieu { + + // Forward declarations + int check_angular_fcn_domain(int m, double q); + int check_modified_fcn_domain(int m, double q); + int set_adaptive_offset_c(int m, double q); + + //================================================================== + int mathieu_ce(int m, double q, double v, double *ce, double *ced) { + // This computes the Mathieu fcn ce + // Inputs: + // m = Mathieu fcn order (scalar) + // q = frequency parameter (scalar) + // v = angle in radians (scalar) + // Outputs: + // ce = value of fcn for these inputs (scalar) + // ced = value of fcn deriv w.r.t. v for these inputs (scalar) + // Return code: + // Codes in error.h. + + int retcode = SF_ERROR_OK; + + // Check inputs. Note that retcode can include SF_ERROR_LOSS + // but the program can keep going in that case. + retcode = check_angular_fcn_domain(m, q); + if (retcode == SF_ERROR_DOMAIN) { + *ce = std::numeric_limits::quiet_NaN(); + *ced = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // I find the peak Fourier coeff tracks m. Therefore, + // adjust the matrix size based on order m. Later make this + // a fcn of q also since the distribution of coeff mags allegedly + // flattens out for large q. + int N = m + 25; // N = size of recursion matrix to use. + + // Use different coeffs depending upon whether m is even or odd. + if (m % 2 == 0) { + // Even order + + // Get coeff vector for even ce + std::vector AA(N); + retcode = mathieu_coeffs_ee(N, q, m, AA.data()); + if (retcode != SF_ERROR_OK) { + *ce = std::numeric_limits::quiet_NaN(); + *ced = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Local scope variables used in summing the Fourier series. + double tt, td, cep, cem, cedp, cedm; + cem = 0.0; + cep = 0.0; + cedm = 0.0; + cedp = 0.0; + + // Sum from smallest to largest coeff. + for (int k = (N - 1); k >= 0; k--) { + tt = AA[k] * cos(2.0 * k * v); // Term for Mathieu ce + if (tt < 0) { + cem = cem + tt; // Neg running sum + } else { + cep = cep + tt; // Pos running sum + } + + td = -2.0 * k * AA[k] * sin(2.0 * k * v); // Term for deriv + if (td < 0) { + cedm = cedm + td; + } else { + cedp = cedp + td; + } + } // for (k=(N-1) ... + + // I should do a sort before doing the sum + *ce = cep + cem; + *ced = cedp + cedm; + + // This makes sure the fcn has the right overall sign for q<0. + // Someday combine this with the above sums into the same for loop. + double s = 0.0; + for (int l = 0; l < N; l++) { + s = s + AA[l]; + } + *ce = SIGN(s) * (*ce); + *ced = SIGN(s) * (*ced); + + } else { + // Odd order + + // Get coeff vector for odd ce + std::vector AA(N); + + retcode = mathieu_coeffs_eo(N, q, m, AA.data()); + if (retcode != SF_ERROR_OK) { + *ce = std::numeric_limits::quiet_NaN(); + *ced = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Local scope variables used in summing the Fourier series. + double tt, td, cep, cem, cedp, cedm; + cem = 0.0; + cep = 0.0; + cedm = 0.0; + cedp = 0.0; + + // Perform Fourier sum on k = 0, 2, 4, ... + for (int k = (N - 1); k >= 0; k--) { + tt = AA[k] * cos((2.0 * k + 1.0) * v); // Term for Mathieu ce + if (tt < 0) { + cem = cem + tt; // Neg running sum + } else { + cep = cep + tt; // Pos running sum + } + + td = -(2.0 * k + 1.0) * AA[k] * sin((2.0 * k + 1.0) * v); // Deriv. + if (td < 0) { + cedm = cedm + td; + } else { + cedp = cedp + td; + } + } // for (k=(N-1) ... + + // I should do a sort before doing the sum + *ce = cep + cem; + *ced = cedp + cedm; + + // This makes sure the fcn has the right overall sign for q<0. + // Someday combine this with the above sums into the same for loop. + double s = 0.0; + for (int l = 0; l < N; l++) { + s = s + AA[l]; + } + *ce = SIGN(s) * (*ce); + *ced = SIGN(s) * (*ced); + + } // if (m % 2 == 0) + + return retcode; + } + + //================================================================== + int mathieu_se(int m, double q, double v, double *se, double *sed) { + // This computes the Mathieu fcn se + // Inputs: + // m = Mathieu fcn order (scalar) + // q = frequency parameter (scalar) + // v = angle in radians (scalar) + // Outputs: + // se = value of fcn for these inputs (scalar) + // sed = value of fcn deriv w.r.t. v for these inputs (scalar) + // Return code: + // Success = 0 + + int retcode = SF_ERROR_OK; + + // Check inputs. If domain err (inputs out of bounds) then return. + // Note that retcode can include SF_ERROR_LOSS + // but the program can keep going in that case. + retcode = check_angular_fcn_domain(m, q); + if (retcode == SF_ERROR_DOMAIN) { + *se = std::numeric_limits::quiet_NaN(); + *sed = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // I find the peak Fourier coeff tracks m. Therefore, + // adjust the matrix size based on order m. Later make this + // a fcn of q also since the distribution of coeff mags allegedly + // flattens out for large q. + int N = m + 25; // N = size of recursion matrix to use. + + // Use different coeffs depending upon whether m is even or odd. + if (m % 2 == 0) { + // Even order. + + // Get coeff vector for even se + std::vector BB(N); + + retcode = mathieu_coeffs_oe(N, q, m, BB.data()); + if (retcode != SF_ERROR_OK) { + *se = std::numeric_limits::quiet_NaN(); + *sed = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Local scope variables used in summing the Fourier series. + double tt, td, sep, sem, sedp, sedm; + sem = 0.0; + sep = 0.0; + sedm = 0.0; + sedp = 0.0; + + // Sum from smallest to largest coeff. + for (int k = N; k >= 1; k--) { + tt = BB[k - 1] * sin(2.0 * k * v); // Mathieu se term + if (tt < 0) { + sem = sem + tt; // Neg running sum + } else { + sep = sep + tt; // Pos running sum + } + + td = 2.0 * k * BB[k - 1] * cos(2.0 * k * v); // Deriv term. + if (td < 0) { + sedm = sedm + td; + } else { + sedp = sedp + td; + } + } // for (k=(N-1) ... + + // I should do a sort before doing the sum + *se = sep + sem; + *sed = sedp + sedm; + + // This makes sure the fcn has the right overall sign for q<0. + // Someday combine this with the above sums into the same for loop. + double s = 0.0; + for (int l = 0; l < N; l++) { + s = s + BB[l]; + } + *se = SIGN(s) * (*se); + *sed = SIGN(s) * (*sed); + + } else { + // Odd order + + // Get coeff vector for odd se + std::vector BB(N); + + retcode = mathieu_coeffs_oo(N, q, m, BB.data()); + if (retcode != SF_ERROR_OK) { + *se = std::numeric_limits::quiet_NaN(); + *sed = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Local scope variables used in summing the Fourier series. + double tt, td, sep, sem, sedp, sedm; + sem = 0.0; + sep = 0.0; + sedm = 0.0; + sedp = 0.0; + + // Sum from smallest to largest coeff. + for (int k = (N - 1); k >= 0; k--) { + tt = BB[k] * sin((2.0 * k + 1.0) * v); // Mathieu se term + if (tt < 0) { + sem = sem + tt; // Neg running sum + } else { + sep = sep + tt; // Pos running sum + } + + td = (2.0 * k + 1.0) * BB[k] * cos((2.0 * k + 1.0) * v); // Deriv term. + if (td < 0) { + sedm = sedm + td; + } else { + sedp = sedp + td; + } + } // for (k=(N-1) ... + + // I should do a sort before doing the sum + *se = sep + sem; + *sed = sedp + sedm; + + // This makes sure the fcn has the right overall sign for q<0. + // Someday combine this with the above sums into the same for loop. + double s = 0.0; + for (int l = 0; l < N; l++) { + s = s + BB[l]; + } + *se = SIGN(s) * (*se); + *sed = SIGN(s) * (*sed); + } + + return retcode; + } // int mathieu_se + + //================================================================== + int mathieu_modmc1(int m, double q, double u, double *mc1, double *mc1d) { + // This computes the Mathieu fcn modmc1 + // Inputs: + // m = Mathieu fcn order (scalar) + // q = frequency parameter (scalar) + // u = radial coord (scalar) + // Outputs: + // mc1 = value of fcn for these inputs (scalar) + // mc1d = value of fcn deriv w.r.t. u for these inputs (scalar) + // Return code: + // Success = 0 + + int retcode = SF_ERROR_OK; + int c; // Offset used in adaptive computation. + + // Check inputs. Note that retcode can include SF_ERROR_LOSS + // but the program can keep going in that case. + retcode = check_modified_fcn_domain(m, q); + if (retcode == SF_ERROR_DOMAIN) { + *mc1 = std::numeric_limits::quiet_NaN(); + *mc1d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // I find the peak Fourier coeff tracks m. Therefore, + // adjust the matrix size based on order m. Later make this + // a fcn of q also since the distribution of coeff mags allegedly + // flattens out for large q. + int N = m + 25; // N = size of recursion matrix to use. + + // Utility vars. + double sqq = sqrt(q); + double exppu = exp(u); + double expmu = exp(-u); + double s = sqq * expmu; + double t = sqq * exppu; + + // Set offset c for adaptive calc. + c = set_adaptive_offset_c(m, q); + + // Use different coeffs depending upon whether m is even or odd. + if (m % 2 == 0) { + // Even order + + // Get coeff vector for even modmc1 + std::vector AA(N); + retcode = mathieu_coeffs_ee(N, q, m, AA.data()); + if (retcode != SF_ERROR_OK) { + *mc1 = std::numeric_limits::quiet_NaN(); + *mc1d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Local scope variables used in summing the Fourier series. + // These are Float128 since some of the terms are near + // equal amplitude, but different sign, and I want to + // avoid catastrophic cancellation. + _Float128 mc1p, mc1m, mc1dp, mc1dm; + mc1p = 0.0; + mc1m = 0.0; + mc1dp = 0.0; + mc1dm = 0.0; + + // Sum from smallest to largest coeff. + for (int k = (N - 1); k >= 0; k--) { + if (c == 0) { + // Non-adaptive calc + double Jks = besselj(k, s); + double Jkt = besselj(k, t); + double Jdks = besseljd(k, s); + double Jdkt = besseljd(k, t); + + _Float128 tt = AA[k] * (Jks * Jkt); + _Float128 ttd = AA[k] * (exppu * Jks * Jdkt - expmu * Jdks * Jkt); + + // Even terms have + sign, odd terms have - sign + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + mc1m = mc1m + tt; + } else { + // Pos terms + mc1p = mc1p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + mc1dm = mc1dm + ttd; + } else { + // Pos terms + mc1dp = mc1dp + ttd; + } + + } else { + // Adaptive calc + double Jkmcs = besselj(k - c, s); + double Jkpcs = besselj(k + c, s); + double Jkpct = besselj(k + c, t); + double Jkmct = besselj(k - c, t); + + double Jdkmcs = besseljd(k - c, s); + double Jdkpcs = besseljd(k + c, s); + double Jdkpct = besseljd(k + c, t); + double Jdkmct = besseljd(k - c, t); + + _Float128 tt = AA[k] * (Jkmcs * Jkpct + Jkpcs * Jkmct); + _Float128 ttd = + AA[k] * (exppu * (Jkmcs * Jdkpct + Jkpcs * Jdkmct) - expmu * (Jdkmcs * Jkpct + Jdkpcs * Jkmct)); + + // Even terms have + sign, odd terms have - sign + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + mc1m = mc1m + tt; + } else { + // Pos terms + mc1p = mc1p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + mc1dm = mc1dm + ttd; + } else { + // Pos terms + mc1dp = mc1dp + ttd; + } + + } // if (c==0) + + } // for (int k=(N-1) ... + + // Sum pos and neg terms to get final result + *mc1 = static_cast(mc1p + mc1m); + *mc1d = static_cast(mc1dp + mc1dm); + + // Do normalization. Note normalization depends upon c. + int sgn = m / 2; + if (sgn % 2 == 0) { + *mc1 = (*mc1) / AA[c]; + *mc1d = sqq * (*mc1d) / AA[c]; + } else { + *mc1 = -(*mc1) / AA[c]; + *mc1d = -sqq * (*mc1d) / AA[c]; + } + + } else { + // Odd order -- m = 1, 3, 5, 7 ... + + // Get coeff vector for odd modmc1 + std::vector AA(N); + retcode = mathieu_coeffs_eo(N, q, m, AA.data()); + if (retcode != SF_ERROR_OK) { + *mc1 = std::numeric_limits::quiet_NaN(); + *mc1d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Variables used in summing the Fourier series. + _Float128 mc1p, mc1m, mc1dp, mc1dm; + mc1p = 0.0; + mc1m = 0.0; + mc1dp = 0.0; + mc1dm = 0.0; + + // Sum from smallest to largest coeff. + for (int k = (N - 1); k >= 0; k--) { + if (c == 0) { + // Non-adaptive calc + double Jks = besselj(k, s); + double Jkp1s = besselj(k + 1, s); + double Jkt = besselj(k, t); + double Jkp1t = besselj(k + 1, t); + + double Jdks = besseljd(k, s); + double Jdkp1s = besseljd(k + 1, s); + double Jdkt = besseljd(k, t); + double Jdkp1t = besseljd(k + 1, t); + + _Float128 tt = AA[k] * (Jks * Jkp1t + Jkp1s * Jkt); + _Float128 ttd = + AA[k] * (exppu * (Jks * Jdkp1t + Jkp1s * Jdkt) - expmu * (Jdks * Jkp1t + Jdkp1s * Jkt)); + + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + mc1m = mc1m + tt; + } else { + // Pos terms + mc1p = mc1p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + mc1dm = mc1dm + ttd; + } else { + // Pos terms + mc1dp = mc1dp + ttd; + } + + } else { + // Adaptive calc + double Jkmcs = besselj(k - c, s); + double Jkpcs = besselj(k + c + 1, s); + double Jkpct = besselj(k + c + 1, t); + double Jkmct = besselj(k - c, t); + + double Jdkmcs = besseljd(k - c, s); + double Jdkpcs = besseljd(k + c + 1, s); + double Jdkpct = besseljd(k + c + 1, t); + double Jdkmct = besseljd(k - c, t); + + _Float128 tt = AA[k] * (Jkmcs * Jkpct + Jkpcs * Jkmct); + _Float128 ttd = + AA[k] * (exppu * (Jkmcs * Jdkpct + Jkpcs * Jdkmct) - expmu * (Jdkmcs * Jkpct + Jdkpcs * Jkmct)); + + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + mc1m = mc1m + tt; + } else { + // Pos terms + mc1p = mc1p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + mc1dm = mc1dm + ttd; + } else { + // Pos terms + mc1dp = mc1dp + ttd; + } + + } // if (c==0) + + } // for (int k=(N-1) ... + + // Sum pos and neg terms to get final answer + *mc1 = static_cast(mc1p + mc1m); + *mc1d = static_cast(mc1dp + mc1dm); + + // Do normalization. Note normalization depends upon c. + int sgn = (m - 1) / 2; + if (sgn % 2 == 0) { + *mc1 = (*mc1) / AA[c]; + *mc1d = sqq * (*mc1d) / AA[c]; + } else { + *mc1 = -(*mc1) / AA[c]; + *mc1d = -sqq * (*mc1d) / AA[c]; + } + } + + return retcode; + } // int mathieu_modmc1 + + //================================================================== + int mathieu_modms1(int m, double q, double u, double *ms1, double *ms1d) { + // This computes the Mathieu fcn modms1 + // Inputs: + // m = Mathieu fcn order (scalar) + // q = frequency parameter (scalar) + // u = radial coord (scalar) + // Outputs: + // ms1 = value of fcn for these inputs (scalar) + // ms1d = value of fcn deriv w.r.t. u for these inputs (scalar) + // Return code: + // Success = 0 + + int retcode = SF_ERROR_OK; + int c; // Offset used in adaptive computation. + + // Check inputs. Note that retcode can include SF_ERROR_LOSS + // but the program can keep going in that case. + retcode = check_modified_fcn_domain(m, q); + if (retcode == SF_ERROR_DOMAIN) { + *ms1 = std::numeric_limits::quiet_NaN(); + *ms1d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // I find the peak Fourier coeff tracks m. Therefore, + // adjust the matrix size based on order m. Later make this + // a fcn of q also since the distribution of coeff mags allegedly + // flattens out for large q. + int N = m + 25; // N = size of recursion matrix to use. + + // Utility vars. + double sqq = sqrt(q); + double exppu = exp(u); + double expmu = exp(-u); + double s = sqq * expmu; + double t = sqq * exppu; + + // Set offset c for adaptive calc. + c = set_adaptive_offset_c(m, q); + + // Use different coeffs depending upon whether m is even or odd. + if (m % 2 == 0) { + // Even order + + // Get coeff vector for even modms1 + std::vector BB(N); + retcode = mathieu_coeffs_oe(N, q, m, BB.data()); + if (retcode != SF_ERROR_OK) { + *ms1 = std::numeric_limits::quiet_NaN(); + *ms1d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Variables used in summing the Fourier series. + // These are Float128 since some of the terms are near + // equal amplitude, but different sign. + _Float128 ms1p, ms1m, ms1dp, ms1dm; + ms1p = 0.0; + ms1m = 0.0; + ms1dp = 0.0; + ms1dm = 0.0; + + // Sum from smallest to largest coeff. + for (int k = (N - 1); k >= 0; k--) { + if (c == 0) { + // Non-adaptive calc + double Jks = besselj(k, s); + double Jkp2t = besselj(k + 2, t); + double Jkp2s = besselj(k + 2, s); + double Jkt = besselj(k, t); + + double Jdks = besseljd(k, s); + double Jdkp2t = besseljd(k + 2, t); + double Jdkp2s = besseljd(k + 2, s); + double Jdkt = besseljd(k, t); + + _Float128 tt = BB[k] * (Jks * Jkp2t - Jkp2s * Jkt); + _Float128 ttd = + BB[k] * (exppu * (Jks * Jdkp2t - Jkp2s * Jdkt) - expmu * (Jdks * Jkp2t - Jdkp2s * Jkt)); + + // Even terms have + sign, odd terms have - sign + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + ms1m = ms1m + tt; + } else { + // Pos terms + ms1p = ms1p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + ms1dm = ms1dm + ttd; + } else { + // Pos terms + ms1dp = ms1dp + ttd; + } + + } else { + // Adaptive calc + double Jkmcs = besselj(k - c, s); + double Jkpct = besselj(k + c + 2, t); + double Jkpcs = besselj(k + c + 2, s); + double Jkmct = besselj(k - c, t); + + double Jdkmcs = besseljd(k - c, s); + double Jdkpct = besseljd(k + c + 2, t); + double Jdkpcs = besseljd(k + c + 2, s); + double Jdkmct = besseljd(k - c, t); + + _Float128 tt = BB[k] * (Jkmcs * Jkpct - Jkpcs * Jkmct); + _Float128 ttd = + BB[k] * (exppu * (Jkmcs * Jdkpct - Jkpcs * Jdkmct) - expmu * (Jdkmcs * Jkpct - Jdkpcs * Jkmct)); + + // Even terms have + sign, odd terms have - sign + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + ms1m = ms1m + tt; + } else { + // Pos terms + ms1p = ms1p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + ms1dm = ms1dm + ttd; + } else { + // Pos terms + ms1dp = ms1dp + ttd; + } + + } // if (c==0) + + } // for (int k=(N-1) ... + + // Sum pos and neg terms to get final result + *ms1 = static_cast(ms1p + ms1m); + *ms1d = static_cast(ms1dp + ms1dm); + + // Do normalization. Note normalization depends upon c. + int sgn = (m - 2) / 2; + if (sgn % 2 == 0) { + *ms1 = (*ms1) / BB[c]; + *ms1d = sqq * (*ms1d) / BB[c]; + } else { + *ms1 = -(*ms1) / BB[c]; + *ms1d = -sqq * (*ms1d) / BB[c]; + } + + } else { + // Odd order -- m = 1, 3, 5, 7 ... + + // Get coeff vector for odd modms1 + std::vector BB(N); + retcode = mathieu_coeffs_oo(N, q, m, BB.data()); + if (retcode != SF_ERROR_OK) { + *ms1 = std::numeric_limits::quiet_NaN(); + *ms1d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Variables used in summing the Fourier series. + _Float128 ms1p, ms1m, ms1dp, ms1dm; + ms1p = 0.0; + ms1m = 0.0; + ms1dp = 0.0; + ms1dm = 0.0; + + // Sum from smallest to largest coeff. + for (int k = (N - 1); k >= 0; k--) { + if (c == 0) { + // Non-adaptive calc + double Jks = besselj(k, s); + double Jkt = besselj(k, t); + double Jkp1s = besselj(k + 1, s); + double Jkp1t = besselj(k + 1, t); + + double Jdks = besseljd(k, s); + double Jdkt = besseljd(k, t); + double Jdkp1s = besseljd(k + 1, s); + double Jdkp1t = besseljd(k + 1, t); + + _Float128 tt = BB[k] * (Jks * Jkp1t - Jkp1s * Jkt); + _Float128 ttd = + BB[k] * (exppu * (Jks * Jdkp1t - Jkp1s * Jdkt) - expmu * (Jdks * Jkp1t - Jdkp1s * Jkt)); + + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + ms1m = ms1m + tt; + } else { + // Pos terms + ms1p = ms1p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + ms1dm = ms1dm + ttd; + } else { + // Pos terms + ms1dp = ms1dp + ttd; + } + + } else { + // Adaptive calc + double Jkmcs = besselj(k - c, s); + double Jkpcs = besselj(k + c + 1, s); + double Jkpct = besselj(k + c + 1, t); + double Jkmct = besselj(k - c, t); + + double Jdkmcs = besseljd(k - c, s); + double Jdkpcs = besseljd(k + c + 1, s); + double Jdkpct = besseljd(k + c + 1, t); + double Jdkmct = besseljd(k - c, t); + + _Float128 tt = BB[k] * (Jkmcs * Jkpct - Jkpcs * Jkmct); + _Float128 ttd = + BB[k] * (exppu * (Jkmcs * Jdkpct - Jkpcs * Jdkmct) - expmu * (Jdkmcs * Jkpct - Jdkpcs * Jkmct)); + + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + ms1m = ms1m + tt; + } else { + // Pos terms + ms1p = ms1p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + ms1dm = ms1dm + ttd; + } else { + // Pos terms + ms1dp = ms1dp + ttd; + } + + } // if (c==0) + + } // for (int k=(N-1) ... + + // Sum pos and neg terms to get final answer + *ms1 = static_cast(ms1p + ms1m); + *ms1d = static_cast(ms1dp + ms1dm); + + // Do normalization. Note normalization depends upon c. + int sgn = (m - 1) / 2; + if (sgn % 2 == 0) { + *ms1 = (*ms1) / BB[c]; + *ms1d = sqq * (*ms1d) / BB[c]; + } else { + *ms1 = -(*ms1) / BB[c]; + *ms1d = -sqq * (*ms1d) / BB[c]; + } + } + + return retcode; + } // int mathieu_modms1 + + //================================================================== + int mathieu_modmc2(int m, double q, double u, double *mc2, double *mc2d) { + // This computes the Mathieu fcn modmc2 + // Inputs: + // m = Mathieu fcn order (scalar) + // q = frequency parameter (scalar) + // u = radial coord (scalar) + // Outputs: + // mc2 = value of fcn for these inputs (scalar) + // mc2d = value of fcn deriv w.r.t. u for these inputs (scalar) + // Return code: + // Success = 0 + + int retcode = SF_ERROR_OK; + int c; // Offset used in adaptive computation. + + // Check inputs. Note that retcode can include SF_ERROR_LOSS + // but the program can keep going in that case. + retcode = check_modified_fcn_domain(m, q); + if (retcode == SF_ERROR_DOMAIN) { + *mc2 = std::numeric_limits::quiet_NaN(); + *mc2d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // I find the peak Fourier coeff tracks m. Therefore, + // adjust the matrix size based on order m. Later make this + // a fcn of q also since the distribution of coeff mags allegedly + // flattens out for large q. + int N = m + 25; // N = size of recursion matrix to use. + + // Utility vars. + double sqq = sqrt(q); + double exppu = exp(u); + double expmu = exp(-u); + double s = sqq * expmu; + double t = sqq * exppu; + + // Set offset c for adaptive calc. + c = set_adaptive_offset_c(m, q); + + // Use different coeffs depending upon whether m is even or odd. + if (m % 2 == 0) { + // Even order + + // Get coeff vector for even modmc2 + std::vector AA(N); + retcode = mathieu_coeffs_ee(N, q, m, AA.data()); + if (retcode != SF_ERROR_OK) { + *mc2 = std::numeric_limits::quiet_NaN(); + *mc2d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Variables used in summing the Fourier series. + // These are Float128 since some of the terms are near + // equal amplitude, but different sign. + _Float128 mc2p, mc2m, mc2dp, mc2dm; + + // Sum from smallest to largest coeff. + mc2p = 0.0; + mc2m = 0.0; + mc2dp = 0.0; + mc2dm = 0.0; + for (int k = (N - 1); k >= 0; k--) { + if (c == 0) { + // Non-adaptive calc + double Jks = besselj(k, s); + double Ykt = bessely(k, t); + double Jdks = besseljd(k, s); + double Ydkt = besselyd(k, t); + + _Float128 tt = AA[k] * Jks * Ykt; + _Float128 ttd = AA[k] * (exppu * Jks * Ydkt - expmu * Jdks * Ykt); + + // Even terms have + sign, odd terms have - sign + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + mc2m = mc2m + tt; + } else { + // Pos terms + mc2p = mc2p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + mc2dm = mc2dm + ttd; + } else { + // Pos terms + mc2dp = mc2dp + ttd; + } + + } else { + // Adaptive calc + double Jkmcs = besselj(k - c, s); + double Jkpcs = besselj(k + c, s); + double Ykpct = bessely(k + c, t); + double Ykmct = bessely(k - c, t); + + double Jdkmcs = besseljd(k - c, s); + double Jdkpcs = besseljd(k + c, s); + double Ydkpct = besselyd(k + c, t); + double Ydkmct = besselyd(k - c, t); + + _Float128 tt = AA[k] * (Jkmcs * Ykpct + Jkpcs * Ykmct); + _Float128 ttd = + AA[k] * (exppu * (Jkmcs * Ydkpct + Jkpcs * Ydkmct) - expmu * (Jdkmcs * Ykpct + Jdkpcs * Ykmct)); + + // Even terms have + sign, odd terms have - sign + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + mc2m = mc2m + tt; + } else { + // Pos terms + mc2p = mc2p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + mc2dm = mc2dm + ttd; + } else { + // Pos terms + mc2dp = mc2dp + ttd; + } + + } // if (c==0) + + } // for (int k=(N-1) ... + + // Sum pos and neg terms to get final result + *mc2 = static_cast(mc2p + mc2m); + *mc2d = static_cast(mc2dp + mc2dm); + + // Do normalization. Note normalization depends upon c. + int sgn = m / 2; + if (sgn % 2 == 0) { + *mc2 = (*mc2) / AA[c]; + *mc2d = sqq * (*mc2d) / AA[c]; + } else { + *mc2 = -(*mc2) / AA[c]; + *mc2d = -sqq * (*mc2d) / AA[c]; + } + + } else { + // Odd order -- m = 1, 3, 5, 7 ... + + // Get coeff vector for odd mc2 + std::vector AA(N); + retcode = mathieu_coeffs_eo(N, q, m, AA.data()); + if (retcode != SF_ERROR_OK) { + *mc2 = std::numeric_limits::quiet_NaN(); + *mc2d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Variables used in summing the Fourier series. + _Float128 mc2p, mc2m, mc2dp, mc2dm; + mc2p = 0.0; + mc2m = 0.0; + mc2dp = 0.0; + mc2dm = 0.0; + + // Sum from smallest to largest coeff. + for (int k = (N - 1); k >= 0; k--) { + if (c == 0) { + // Non-adaptive calc + double Jks = besselj(k, s); + double Ykt = bessely(k, t); + double Jkp1s = besselj(k + 1, s); + double Ykp1t = bessely(k + 1, t); + + double Jdks = besseljd(k, s); + double Ydkt = besselyd(k, t); + double Jdkp1s = besseljd(k + 1, s); + double Ydkp1t = besselyd(k + 1, t); + + _Float128 tt = AA[k] * (Jks * Ykp1t + Jkp1s * Ykt); + _Float128 ttd = + AA[k] * (exppu * (Jks * Ydkp1t + Jkp1s * Ydkt) - expmu * (Jdks * Ykp1t + Jdkp1s * Ykt)); + + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + mc2m = mc2m + tt; + } else { + // Pos terms + mc2p = mc2p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + mc2dm = mc2dm + ttd; + } else { + // Pos terms + mc2dp = mc2dp + ttd; + } + + } else { + // Adaptive calc + double Jkmcs = besselj(k - c, s); + double Jkpcs = besselj(k + c + 1, s); + double Ykpct = bessely(k + c + 1, t); + double Ykmct = bessely(k - c, t); + + double Jdkmcs = besseljd(k - c, s); + double Jdkpcs = besseljd(k + c + 1, s); + double Ydkpct = besselyd(k + c + 1, t); + double Ydkmct = besselyd(k - c, t); + + _Float128 tt = AA[k] * (Jkmcs * Ykpct + Jkpcs * Ykmct); + _Float128 ttd = + AA[k] * (exppu * (Jkmcs * Ydkpct + Jkpcs * Ydkmct) - expmu * (Jdkmcs * Ykpct + Jdkpcs * Ykmct)); + + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + mc2m = mc2m + tt; + } else { + // Pos terms + mc2p = mc2p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + mc2dm = mc2dm + ttd; + } else { + // Pos terms + mc2dp = mc2dp + ttd; + } + + } // if (c==0) + + } // for (int k=(N-1) ... + + // Sum pos and neg terms to get final answer + *mc2 = static_cast(mc2p + mc2m); + *mc2d = static_cast(mc2dp + mc2dm); + + // Do normalization. Note normalization depends upon c. + int sgn = (m - 1) / 2; + if (sgn % 2 == 0) { + *mc2 = (*mc2) / AA[c]; + *mc2d = sqq * (*mc2d) / AA[c]; + } else { + *mc2 = -(*mc2) / AA[c]; + *mc2d = -sqq * (*mc2d) / AA[c]; + } + } + + return retcode; + } // int mathieu_modmc2 + + //================================================================== + int mathieu_modms2(int m, double q, double u, double *ms2, double *ms2d) { + // This computes the Mathieu fcn modms2 + // Inputs: + // m = Mathieu fcn order (scalar) + // q = frequency parameter (scalar) + // u = radial coord (scalar) + // Outputs: + // ms2 = value of fcn for these inputs (scalar) + // ms2d = value of fcn deriv w.r.t. u for these inputs (scalar) + // Return code: + // Success = 0 + + int retcode = SF_ERROR_OK; + int c; // Offset used in adaptive computation. + + // Check inputs. Note that retcode can include SF_ERROR_LOSS + // but the program can keep going in that case. + retcode = check_modified_fcn_domain(m, q); + if (retcode == SF_ERROR_DOMAIN) { + *ms2 = std::numeric_limits::quiet_NaN(); + *ms2d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // I find the peak Fourier coeff tracks m. Therefore, + // adjust the matrix size based on order m. Later make this + // a fcn of q also since the distribution of coeff mags allegedly + // flattens out for large q. + int N = m + 25; // N = size of recursion matrix to use. + + // Utility vars. + double sqq = sqrt(q); + double exppu = exp(u); + double expmu = exp(-u); + double s = sqq * expmu; + double t = sqq * exppu; + + // Set offset c for adaptive calc. + // c = set_adaptive_offset_c(m, q); + c = 0; // Turn off adaptive c in modms2 for now ... + + // Use different coeffs depending upon whether m is even or odd. + if (m % 2 == 0) { + // Even order m. + + // Get coeff vector for even modms2 + std::vector BB(N); + retcode = mathieu_coeffs_oe(N, q, m, BB.data()); + if (retcode != SF_ERROR_OK) { + *ms2 = std::numeric_limits::quiet_NaN(); + *ms2d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Variables used in summing the Fourier series. + // These are Float128 since some of the terms are near + // equal amplitude, but different sign. + _Float128 ms2p, ms2m, ms2dp, ms2dm; + ms2p = 0.0; + ms2m = 0.0; + ms2dp = 0.0; + ms2dm = 0.0; + + // Sum from smallest to largest coeff. + for (int k = (N - 1); k >= 0; k--) { + if (c == 0) { + // Non-adaptive calc + double Jks = besselj(k, s); + double Ykp2t = bessely(k + 2, t); + double Jkp2s = besselj(k + 2, s); + double Ykt = bessely(k, t); + + double Jdks = besseljd(k, s); + double Ydkp2t = besselyd(k + 2, t); + double Jdkp2s = besseljd(k + 2, s); + double Ydkt = besselyd(k, t); + + _Float128 tt = BB[k] * (Jks * Ykp2t - Jkp2s * Ykt); + _Float128 ttd = + BB[k] * (exppu * (Jks * Ydkp2t - Jkp2s * Ydkt) - expmu * (Jdks * Ykp2t - Jdkp2s * Ykt)); + + // Even terms have + sign, odd terms have - sign + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + ms2m = ms2m + tt; + } else { + // Pos terms + ms2p = ms2p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + ms2dm = ms2dm + ttd; + } else { + // Pos terms + ms2dp = ms2dp + ttd; + } + + } else { + // Adaptive calc + double Jkmcs = besselj(k - c, s); + double Ykpct = bessely(k + c + 2, t); + double Jkpcs = besselj(k + c + 2, s); + double Ykmct = bessely(k - c, t); + + double Jdkmcs = besseljd(k - c, s); + double Ydkpct = besselyd(k + c + 2, t); + double Jdkpcs = besseljd(k + c + 2, s); + double Ydkmct = besselyd(k - c, t); + + _Float128 tt = BB[k] * (Jkmcs * Ykpct - Jkpcs * Ykmct); + _Float128 ttd = + BB[k] * (exppu * (Jkmcs * Ydkpct - Jkpcs * Ydkmct) - expmu * (Jdkmcs * Ykpct - Jdkpcs * Ykmct)); + + // Even terms have + sign, odd terms have - sign + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + ms2m = ms2m + tt; + } else { + // Pos terms + ms2p = ms2p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + ms2dm = ms2dm + ttd; + } else { + // Pos terms + ms2dp = ms2dp + ttd; + } + + } // if (c==0) + + } // for (int k=(N-1) ... + + // Sum pos and neg terms to get final result + *ms2 = static_cast(ms2p + ms2m); + *ms2d = static_cast(ms2dp + ms2dm); + + // Do normalization. Note normalization depends upon c. + int sgn = (m - 2) / 2; + if (sgn % 2 == 0) { + *ms2 = (*ms2) / BB[c]; + *ms2d = sqq * (*ms2d) / BB[c]; + } else { + *ms2 = -(*ms2) / BB[c]; + *ms2d = -sqq * (*ms2d) / BB[c]; + } + + } else { + // Odd order -- m = 1, 3, 5, 7 ... + + // Get coeff vector for odd modms2 + std::vector BB(N); + retcode = mathieu_coeffs_oo(N, q, m, BB.data()); + if (retcode != SF_ERROR_OK) { + *ms2 = std::numeric_limits::quiet_NaN(); + *ms2d = std::numeric_limits::quiet_NaN(); + return retcode; + } + + // Variables used in summing the Fourier series. + _Float128 ms2p, ms2m, ms2dp, ms2dm; + ms2p = 0.0; + ms2m = 0.0; + ms2dp = 0.0; + ms2dm = 0.0; + + // Sum from smallest to largest coeff. + for (int k = (N - 1); k >= 0; k--) { + if (c == 0) { + // Non-adaptive calc + double Jks = besselj(k, s); + double Ykt = bessely(k, t); + double Jkp1s = besselj(k + 1, s); + double Ykp1t = bessely(k + 1, t); + + double Jdks = besseljd(k, s); + double Ydkt = besselyd(k, t); + double Jdkp1s = besseljd(k + 1, s); + double Ydkp1t = besselyd(k + 1, t); + + _Float128 tt = BB[k] * (Jks * Ykp1t - Jkp1s * Ykt); + _Float128 ttd = + BB[k] * (exppu * (Jks * Ydkp1t - Jkp1s * Ydkt) - expmu * (Jdks * Ykp1t - Jdkp1s * Ykt)); + + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + ms2m = ms2m + tt; + } else { + // Pos terms + ms2p = ms2p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + ms2dm = ms2dm + ttd; + } else { + // Pos terms + ms2dp = ms2dp + ttd; + } + + } else { + // Adaptive calc + double Jkmcs = besselj(k - c, s); + double Jkpcs = besselj(k + c + 1, s); + double Ykpct = bessely(k + c + 1, t); + double Ykmct = bessely(k - c, t); + + double Jdkmcs = besseljd(k - c, s); + double Jdkpcs = besseljd(k + c + 1, s); + double Ydkpct = besselyd(k + c + 1, t); + double Ydkmct = besselyd(k - c, t); + + _Float128 tt = BB[k] * (Jkmcs * Ykpct - Jkpcs * Ykmct); + _Float128 ttd = + BB[k] * (exppu * (Jkmcs * Ydkpct - Jkpcs * Ydkmct) - expmu * (Jdkmcs * Ykpct - Jdkpcs * Ykmct)); + + int sgn = (k % 2 == 0) ? 1 : -1; + + // Do sum using separate sums for + and - + tt = sgn * tt; + if (tt < 0) { + // Neg terms + ms2m = ms2m + tt; + } else { + // Pos terms + ms2p = ms2p + tt; + } + + // Do sum using separate sums for + and - + ttd = sgn * ttd; + if (ttd < 0) { + // Neg terms + ms2dm = ms2dm + ttd; + } else { + // Pos terms + ms2dp = ms2dp + ttd; + } + + } // if (c==0) + + } // for (int k=(N-1) ... + + // Sum pos and neg terms to get final answer + *ms2 = static_cast(ms2p + ms2m); + *ms2d = static_cast(ms2dp + ms2dm); + + // Do normalization. Note normalization depends upon c. + int sgn = (m - 1) / 2; + if (sgn % 2 == 0) { + *ms2 = (*ms2) / BB[c]; + *ms2d = sqq * (*ms2d) / BB[c]; + } else { + *ms2 = -(*ms2) / BB[c]; + *ms2d = -sqq * (*ms2d) / BB[c]; + } + } + + return retcode; + } // int mathieu_modms2 + + //================================================================ + // Helper fcns -- these help reduce the amount of redundant code. + + int check_angular_fcn_domain(int m, double q) { + // Verify inputs are OK. If not indicate err. + int retcode = SF_ERROR_OK; + + if (m > 500) { + // Don't support absurdly larger orders for now. + return SF_ERROR_DOMAIN; + } + + // abs(q) > 1000 leads to low accuracy. + if (abs(q) > 1.0e3d) + retcode = SF_ERROR_LOSS; + + return retcode; + } + + //--------------------------------------------------- + int check_modified_fcn_domain(int m, double q) { + int retcode = SF_ERROR_OK; + + // Check input domain and flag any problems + if (m > 500) { + return SF_ERROR_DOMAIN; + } + + if (q < 0) { + return SF_ERROR_DOMAIN; // q<0 is unimplemented + } + + // Don't need to bail out of main computation for these, just set retcode. + if (abs(q) > 1.0e3d) + retcode = SF_ERROR_LOSS; // q>1000 is inaccurate + if (m > 15 && q > 0.1d) + retcode = SF_ERROR_LOSS; + + return retcode; + } + + //--------------------------------------------------- + int set_adaptive_offset_c(int m, double q) { + // This is used to set the c used in the adaptive computation. + // I set the offset used in Bessel fcn depending upon order m + // and shape/frequency parameter q. This improves the accuracy + // for larger values of m. + // The idea comes from the book "Accurate Computation of Mathieu Functions", + // Malcolm M. Bibby & Andrew F. Peterson. Also used in the paper + // "Accurate calculation of the modified Mathieu functions of + // integer order", Van Buren & Boisvert. The values I use here + // were found from experiment using my Matlab prototype. However, + // better values are likely -- finding them is a future project. + int c; + + if ((m > 5 && q < .001) || (m > 7 && q < .01) || (m > 10 && q < .1) || (m > 15 && q < 1) || + (m > 20 && q < 10) || (m > 30 && q < 100)) { + c = m / 2; + } else { + c = 0; + } + + return c; + } + +} // namespace mathieu +} // namespace xsf + +#endif // #ifndef MATHIEU_FCNS_H diff --git a/include/xsf/mathieu/matrix_utils.h b/include/xsf/mathieu/matrix_utils.h new file mode 100644 index 000000000..ee151ba58 --- /dev/null +++ b/include/xsf/mathieu/matrix_utils.h @@ -0,0 +1,88 @@ +#ifndef MATRIX_UTILS_H +#define MATRIX_UTILS_H + +#include "matrix_utils.h" +#include +#include + +// These fcns are meant to make it easier to deal with +// matrices in C/C++. We use col major format since that's +// what underlies Lapack. + +// returns +/-1 depending upon sign of x +#define SIGN(x) (((x) > 0) - ((x) < 0)) + +// Macros to extract matrix index and element. +// Matrix is NxN, i = row idx, j = col idx. +// MATRIX_IDX is where col major format is enforced. +#define MATRIX_IDX(N, I, J) (((N) * (I)) + (J)) +#define MATRIX_ELEMENT(A, m, n, i, j) A[MATRIX_IDX(n, i, j)] + +// Min and max macros for scalars. +#define MIN(a, b) (((a) < (b)) ? (a) : (b)) +#define MAX(a, b) (((a) > (b)) ? (a) : (b)) + +//=========================================================== +// This file holds utility functions for dealing with vectors +// and matrices. The idea is to be able to reuse common matrix +// operations. I will name the utils analogously to their names +// in Matlab. +// Note that C matrices are row-major. + +namespace xsf { +namespace mathieu { + + //----------------------------------------------------- + void print_matrix(const double *A, int m, int n) { + // prints matrix as 2-dimensional tablei -- this is how we + // usually think of matrices. + int i, j; + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf("% 10.4e ", MATRIX_ELEMENT(A, m, n, i, j)); + } + printf("\n"); + } + } + + //----------------------------------------------------- + // Stuff to sort a vector. + // Function to swap two elements + void swap(double *a, double *b) { + double temp = *a; + *a = *b; + *b = temp; + } + + // Partition function for quicksort + int partition(double *arr, int low, int high) { + double pivot = arr[high]; // Choose last element as pivot + int i = (low - 1); // Index of smaller element + + for (int j = low; j <= high - 1; j++) { + // If current element is smaller than or equal to pivot + if (arr[j] <= pivot) { + i++; + swap(&arr[i], &arr[j]); + } + } + swap(&arr[i + 1], &arr[high]); + return (i + 1); + } + + // Quicksort function + void quickSort(double *arr, int low, int high) { + if (low < high) { + // Partition the array and get pivot index + int pivotIndex = partition(arr, low, high); + + // Recursively sort elements before and after partition + quickSort(arr, low, pivotIndex - 1); + quickSort(arr, pivotIndex + 1, high); + } + } + +} // namespace mathieu +} // namespace xsf + +#endif // #ifndef MATRIX_UTILS_H diff --git a/tests/scipy_special_tests/test_mathieu.cpp b/tests/scipy_special_tests/test_mathieu.cpp new file mode 100644 index 000000000..3c37f144e --- /dev/null +++ b/tests/scipy_special_tests/test_mathieu.cpp @@ -0,0 +1,1130 @@ +#include "../testing_utils.h" + +#include + +#include +#include + +/* + * + * The goal of these tests are to verify that my C/C++ + * impl of the Mathieu fcns has been carried over from + * my Matlab impl correctly. Therefore, main just calls + * a bunch of golden value tests. The GVs were generated + * using the Matlab impl. I did fairly extensive + * validation of the Matlab impl. Therefore, if the tests + * here show the C impl matches the Matlab impl, then that + * should serve as verification of the C impl's correctness. + * + * A secondary goal is to show how to call the various fcns + * in my API. + * + */ + +namespace xsf { +namespace mathieu { + + //------------------------------------------------------------- + extern "C" int main() { + int N = 6; + int pass = 0; + int fail = 0; + + //******************************************************* + // First print out the recursion matrices. These are + // private -- not public fcns. But I want to see they + // are correct. + //******************************************************* + { + double *A = (double *)calloc(N * N, sizeof(double)); + double q = 2.0d; + + make_matrix_ee(N, q, A); + print_matrix(A, N, N); + printf("----------------------------------------------\n"); + + make_matrix_eo(N, q, A); + print_matrix(A, N, N); + printf("----------------------------------------------\n"); + + make_matrix_oe(N, q, A); + print_matrix(A, N, N); + printf("----------------------------------------------\n"); + + make_matrix_oo(N, q, A); + print_matrix(A, N, N); + printf("----------------------------------------------\n"); + + free(A); + } + + //******************************************************* + // Try computing the eigenvalues. These are public fcns. + //******************************************************* + ///* + printf("==============================================\n"); + printf("Test a eigenvalues\n"); + double a; + { + double q = 0.001; + double tol = 1e-13; + + // Golden values from Matlab. + double a_true[6] = {-4.999999453125127e-07, 1.000999874984374, 4.000000416666611, + 9.000000062515628, 16.000000033333333, 25.000000020833337}; + printf("q = %f\n", q); + for (int m = 0; m < N; m++) { + mathieu_a(m, q, &a); + printf("Order m = %d, eigenvalue a = %8.5f\n", m, a); + printf("Checking, abs tol = %e ... ", tol); + if (abs(a - a_true[m]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + printf("----------------------------------------------\n"); + { + double q = 1.0; + double tol = 1e-13; + double a_true[6] = {-0.455138604107414, 1.859108072514364, 4.371300982735085, + 9.078368847203102, 16.033832340359510, 25.020854345448583}; + printf("q = %f\n", q); + for (int m = 0; m < N; m++) { + mathieu_a(m, q, &a); + printf("Order m = %d, eigenvalue a = %8.5f\n", m, a); + printf("Checking, abs tol = %e ... ", tol); + if (abs(a - a_true[m]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + printf("----------------------------------------------\n"); + { + double q = 100.0; + double tol = 1e-13; + double a_true[6] = {-1.802532491522514e+02, -1.412800568086200e+02, -1.033705070649672e+02, + -66.574389967485317, -30.950104000623185, 3.432562047992776}; + q = 100.0; + printf("q = %f\n", q); + for (int m = 0; m < N; m++) { + mathieu_a(m, q, &a); + printf("Order m = %d, eigenvalue a = %8.5f\n", m, a); + printf("Checking, abs tol = %e ... ", tol); + if (abs(a - a_true[m]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + //*/ + + ///* + printf("==============================================\n"); + double b; + printf("Test b eigenvalues\n"); + { + double q = 0.001; + double tol = 1e-13; + double b_true[6] = {0.998999875015624, 3.999999916666667, 9.000000062484375, + 16.000000033333333, 25.000000020833333, 36.000000014285725}; + printf("q = %f\n", q); + for (int m = 1; m < N; m++) { + mathieu_b(m, q, &b); + printf("Order m = %d, eigenvalue b = %8.5f\n", m, b); + printf("Checking, abs tol = %e ... ", tol); + if (abs(b - b_true[m - 1]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + printf("----------------------------------------------\n"); + { + double q = 1.0; + double tol = 1e-13; + double b_true[6] = {-0.110248816992095, 3.917024772998471, 9.047739259809376, + 16.032970081405796, 25.020840823289770, 36.014289910628236}; + printf("q = %f\n", q); + for (int m = 1; m < N; m++) { + mathieu_b(m, q, &b); + printf("Order m = %d, eigenvalue b = %8.5f\n", m, b); + printf("Checking, abs tol = %e ... ", tol); + if (abs(b - b_true[m - 1]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + printf("----------------------------------------------\n"); + { + double q = 100.0; + double tol = 1e-13; + double b_true[6] = {-1.802532491522514e+02, -1.412800568086194e+02, -1.033705070649312e+02, + -66.574389965842357, -30.950103947238059, 3.432563359324842}; + printf("q = %f\n", q); + for (int m = 1; m < N; m++) { + mathieu_b(m, q, &b); + printf("Order m = %d, eigenvalue b = %8.5f\n", m, b); + printf("Checking, abs tol = %e ... ", tol); + if (abs(b - b_true[m - 1]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + //*/ + + //******************************************************* + // Compute the Fourier coeffs -- these are public fcns. + //******************************************************* + printf("==============================================\n"); + printf("Test A Fourier coeffs\n"); + int retcode; + int m; + double *AA = (double *)calloc(N, sizeof(double)); + N = 5; + { + double q = 1.0d; + double tol = 1e-13; + double AA_true[5] = { + 6.72989672316633203e-01, -3.06303580036475842e-01, 1.86455593633495648e-02, -5.11683638542392003e-04, + 7.93860116701059744e-06 + }; + + m = 0; + printf("ee, N = %d, q = %f, m = %d\n", N, q, m); + retcode = mathieu_coeffs_ee(N, q, m, AA); + if (retcode == 0) { + pass++; + printf("Function call succeeded!\n"); + } else { + printf("Call to mathieu_coeffs_ee returned failure!\n"); + fail++; + } + + printf("Checking values ... \n"); + for (int j = 0; j < N; j++) { + printf("AA[%d] = % 9.6e ...", j, AA[j]); + // I test the relative err here since the coeffs vary widely + if (abs((AA[j] - AA_true[j]) / AA[j]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! AA_true = % 9.6e\n", AA_true[j]); + } + } + } + printf("----------------------------------------------\n"); + { + double q = 1.0d; + double tol = 1e-13; + m = 1; + double AA_true[5] = { + 9.90202059407947477e-01, -1.39511476750210750e-01, 6.03431870922933183e-03, -1.28040356069675752e-04, + 1.61787860802725383e-06 + }; + printf("eo, N = %d, q = %f, m = %d\n", N, q, m); + retcode = mathieu_coeffs_eo(N, q, m, AA); + if (retcode == 0) { + pass++; + printf("Function call succeeded!\n"); + } else { + printf("Call to mathieu_coeffs_eo returned failure!\n"); + fail++; + } + + printf("Checking values ... \n"); + for (int j = 0; j < N; j++) { + printf("AA[%d] = % 9.6e ...", j, AA[j]); + if (abs((AA[j] - AA_true[j]) / AA[j]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! AA_true = % 9.6e\n", AA_true[j]); + } + } + } + printf("----------------------------------------------\n"); + { + double q = 1.0d; + double tol = 1e-13; + m = 2; + double AA_true[5] = { + 2.16927946751974854e-01, 9.48257346823881631e-01, -8.17670087238141774e-02, 2.58658716581847224e-03, + -4.33782257276887011e-05 + }; + printf("N = %d, q = %f, m = %d\n", N, q, m); + retcode = mathieu_coeffs_ee(N, q, m, AA); + if (retcode == 0) { + pass++; + printf("Function call succeeded!\n"); + } else { + printf("Call to mathieu_coeffs_ee returned failure!\n"); + fail++; + } + + printf("Checking values ... \n"); + for (int j = 0; j < N; j++) { + printf("AA[%d] = % 9.6e ...", j, AA[j]); + if (abs((AA[j] - AA_true[j]) / AA[j]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! AA_true = % 9.6e\n", AA_true[j]); + } + } + } + + printf("----------------------------------------------\n"); + { + double q = 1.0d; + double tol = 1e-13; + m = 3; + double AA_true[5] = { + 1.39615654604830663e-01, 9.88251100137287120e-01, -6.21675551357310924e-02, 1.55778240472689887e-03, + -2.16594420865885830e-05 + }; + + printf("N = %d, q = %f, m = %d\n", N, q, m); + retcode = mathieu_coeffs_eo(N, q, m, AA); + if (retcode == 0) { + pass++; + printf("Function call succeeded!\n"); + } else { + printf("Call to mathieu_coeffs_eo returned failure!\n"); + fail++; + } + + printf("Checking values ... \n"); + for (int j = 0; j < N; j++) { + printf("AA[%d] = % 9.6e ...", j, AA[j]); + if (abs((AA[j] - AA_true[j]) / AA[j]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! AA_true = % 9.6e\n", AA_true[j]); + } + } + } + printf("----------------------------------------------\n"); + { + printf("This should fail -- calling coeffs_eo with even m\n"); + double q = 1.0d; + m = 4; + printf("N = %d, q = %f, m = %d\n", N, q, m); + retcode = mathieu_coeffs_eo(N, q, m, AA); + if (retcode != 0) { + pass++; + printf("Function call correctly returned error!\n"); + } else { + printf("Bad call to mathieu_coeffs_eo returned success -- bad!\n"); + fail++; + } + } + + //===================================================== + printf("==============================================\n"); + printf("Test B Fourier coeffs\n"); + { + printf("This should fail -- calling coeffs_oe with m=0\n"); + double q = 1.0d; + m = 0; + printf("N = %d, q = %f, m = %d\n", N, q, m); + retcode = mathieu_coeffs_oe(N, q, m, AA); + if (retcode != 0) { + pass++; + printf("Function call correctly returned error!\n"); + } else { + printf("Bad call to mathieu_coeffs_eo returned success -- bad!\n"); + fail++; + } + } + + printf("----------------------------------------------\n"); + { + double q = 1.0d; + double tol = 1e-13; + m = 1; + double AA_true[5] = { + 9.93967961398935729e-01, -1.09583791872267272e-01, 4.36764886689439396e-03, -8.89579207045385789e-05, + 1.09675314774650835e-06 + }; + printf("N = %d, q = %f, m = %d\n", N, q, m); + retcode = mathieu_coeffs_oo(N, q, m, AA); + if (retcode == 0) { + pass++; + printf("Function call succeeded!\n"); + } else { + printf("Call to mathieu_coeffs_oo returned failure!\n"); + fail++; + } + + printf("Checking values ... \n"); + for (int j = 0; j < N; j++) { + printf("AA[%d] = % 9.6e ...", j, AA[j]); + if (abs((AA[j] - AA_true[j]) / AA[j]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! AA_true = % 9.6e\n", AA_true[j]); + } + } + } + printf("----------------------------------------------\n"); + { + double q = 1.0d; + double tol = 1e-13; + m = 2; + double AA_true[5] = { + 9.96571915618007398e-01, -8.26907809217511391e-02, 2.57874176092240956e-03, + -4.29271107568478840e-05, 4.46771248032548085e-07, + }; + + printf("N = %d, q = %f, m = %d\n", N, q, m); + retcode = mathieu_coeffs_oe(N, q, m, AA); + if (retcode == 0) { + pass++; + printf("Function call succeeded!\n"); + } else { + printf("Call to mathieu_coeffs_oe returned failure!\n"); + fail++; + } + + printf("Checking values ... \n"); + for (int j = 0; j < N; j++) { + printf("AA[%d] = % 9.6e ...", j, AA[j]); + if (abs((AA[j] - AA_true[j]) / AA[j]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! AA_true = % 9.6e\n", AA_true[j]); + } + } + } + printf("----------------------------------------------\n"); + { + double q = 1.0d; + double tol = 1e-13; + m = 3; + double AA_true[5] = { + 1.09642473301236554e-01, 9.92016510230659843e-01, -6.22843393799822065e-02, 1.55951158907772385e-03, + -2.16742541934712708e-05 + }; + + printf("N = %d, q = %f, m = %d\n", N, q, m); + retcode = mathieu_coeffs_oo(N, q, m, AA); + if (retcode == 0) { + pass++; + printf("Function call succeeded!\n"); + } else { + printf("Call to mathieu_coeffs_oo returned failure!\n"); + fail++; + } + + printf("Checking values ... \n"); + for (int j = 0; j < N; j++) { + printf("AA[%d] = % 9.6e ...", j, AA[j]); + if (abs((AA[j] - AA_true[j]) / AA[j]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! AA_true = % 9.6e\n", AA_true[j]); + } + } + } + printf("----------------------------------------------\n"); + { + double q = 1.0d; + double tol = 1e-13; + m = 4; + double AA_true[5] = { + 8.27162782989298156e-02, 9.95322502016357302e-01, -4.99004143812368725e-02, 1.04056488398845915e-03, + -1.23925412747991493e-05 + }; + + printf("N = %d, q = %f, m = %d\n", N, q, m); + retcode = mathieu_coeffs_oe(N, q, m, AA); + if (retcode == 0) { + pass++; + printf("Function call succeeded!\n"); + } else { + printf("Call to mathieu_coeffs_oe returned failure!\n"); + fail++; + } + + printf("Checking values ... \n"); + for (int j = 0; j < N; j++) { + printf("AA[%d] = % 9.6e ...", j, AA[j]); + if (abs((AA[j] - AA_true[j]) / AA[j]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! AA_true = % 9.6e\n", AA_true[j]); + } + } + } + + //******************************************************* + // Compute the Mathieu fcns & compare the results to + // those from the Matlab impls. These are + // the publically-accessible Mathieu fcn impls. + //******************************************************* + printf("==============================================\n"); + printf("Test Mathieu ce\n"); + double ce; + double ced; + N = 6; + { + double v = 1.0d; + double q = 0.001; + double tol = 1e-13; + // Golden values obtained from Matlab for m = 0, 1, 2, 3, 4, 5. + double ce_true[6] = {0.707253852647939, 0.540426067653929, -0.415842336334926, + -0.989942668408815, -0.653726300137563, 0.283568898267299}; + double ced_true[6] = {6.430371575928823e-04, -0.841418026637729, -1.818846996799653, + -0.423764888085533, 3.026974584470965, 4.794786515927561}; + + printf("q = %f\n", q); + for (int m = 0; m < N; m++) { + retcode = mathieu_ce(m, q, v, &ce, &ced); + printf("Order m = %d, v = %f, ce = %8.5f, ced = %8.5f\n", m, v, ce, ced); + printf("Checking, abs tol = %e ... ", tol); + if (abs(ce - ce_true[m]) < tol && abs(ced - ced_true[m]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + printf("----------------------------------------------\n"); + { + double v = 1.0d; + double q = 100; + double tol = 1e-13; + double ce_true[6] = {0.086725459442521, 0.315928563892497, 0.737195938153019, + 1.233155548880747, 1.476584096867507, 1.114354874812387}; + double ced_true[6] = { + 0.924256230828739, 2.780839256567336, 4.896513216077115, 4.882886873958746, + 0.264149721572143, -7.943564372635874 + + }; + printf("q = %f\n", q); + for (int m = 0; m < N; m++) { + retcode = mathieu_ce(m, q, v, &ce, &ced); + printf("Order m = %d, v = %f, ce = %8.5f, ced = %8.5f\n", m, v, ce, ced); + printf("Checking, abs tol = %e ... ", tol); + if (abs(ce - ce_true[m]) < tol && abs(ced - ced_true[m]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + + printf("==============================================\n"); + printf("Test Mathieu se\n"); + double se; + double sed; + N = 6; + { + double v = 1.0d; + double q = 0.001; + double tol = 1e-13; + // Golden values obtained from Matlab for m = 1, 2, 3, 4, 5, 6. + double se_true[6] = {0.841453335445950, 0.909360489814828, 0.141285111200371, + -0.756712745143668, -0.958942823901047, -0.279488670902052}; + double sed_true[6] = {0.540673509812678, -0.832075773996413, -2.969998567644164, + -2.614931881211228, 1.417905406870398, 5.760932545758129}; + + printf("q = %f\n", q); + for (int m = 1; m <= N; m++) { + retcode = mathieu_se(m, q, v, &se, &sed); + printf("Order m = %d, v = %f, se = %8.5f, sed = %8.5f\n", m, v, se, sed); + printf("Checking, abs tol = %e ... ", tol); + if (abs(se - se_true[m - 1]) < tol && abs(sed - sed_true[m - 1]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + printf("----------------------------------------------\n"); + { + double v = 1.0d; + double q = 100; + double tol = 1e-13; + // Golden values obtained from Matlab for m = 1, 2, 3, 4, 5, 6. + double se_true[6] = {0.086725459442519, 0.315928563892424, 0.737195938151456, + 1.233155548910162, 1.476584099964093, 1.114354958780249}; + double sed_true[6] = {0.924256230828753, 2.780839256567937, 4.896513216103436, + 4.882886874949449, 0.264149739778980, -7.943564564030449}; + + printf("q = %f\n", q); + for (int m = 1; m <= N; m++) { + retcode = mathieu_se(m, q, v, &se, &sed); + printf("Order m = %d, v = %f, se = %8.5f, sed = %8.5f\n", m, v, se, sed); + printf("Checking, abs tol = %e ... ", tol); + if (abs(se - se_true[m - 1]) < tol && abs(sed - sed_true[m - 1]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + + //******************************************************* + // Check that my Bessel fcns return the right thing + //******************************************************* + printf("==============================================\n"); + printf("Test my BesselJ implementation ... \n"); + double J; + double Y; + { + // Golden values obtained from Matlab for z = 1 and m = -3, -2, -1, 0, 1, 2, 3.. + double z = 1.0d; + double tol = 1e-13; + double J_true[7] = {-0.019563353982668, 0.114903484931901, -0.440050585744934, 0.765197686557967, + 0.440050585744934, 0.114903484931901, 0.019563353982668}; + + for (int m = -3; m <= 3; m++) { + J = besselj(m, z); + printf("Order m = %d, z = %f, J = %20.17f, J_true = %20.17e\n", m, z, J, J_true[m + 3]); + printf("Checking, abs tol = %e ... ", tol); + if (abs(J - J_true[m + 3]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + + printf("------------------------------------------------\n"); + printf("Test my BesselY implementation, z = 1 ... \n"); + { + // Golden values obtained from Matlab for z = 1 and m = -3, -2, -1, 0, 1, 2, 3.. + double z = 1.0d; + double tol = 1e-13; + double Y_true[7] = {5.821517605964731, -1.650682606816255, 0.781212821300289, 0.088256964215677, + -0.781212821300289, -1.650682606816255, -5.821517605964731}; + + for (int m = -3; m <= 3; m++) { + Y = bessely(m, z); + printf("Order m = %d, z = %f, Y = %20.17e, Y_true = %20.17e\n", m, z, Y, Y_true[m + 3]); + printf("Checking, abs tol = %e ... ", tol); + if (abs(Y - Y_true[m + 3]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + + printf("------------------------------------------------\n"); + printf("Test my BesselY implementation, z = .01 ... \n"); + { + // Golden values obtained from Matlab for z = 1 and m = -3, -2, -1, 0, 1, 2, 3.. + double z = 0.01d; + double tol = 1e-13; + double Y_true[7] = {5.093021841713738e+06, -1.273271380077505e+04, 63.678596282060660, + -3.005455637083646, -63.678596282060660, -1.273271380077505e+04, + -5.093021841713738e+06}; + + for (int m = -3; m <= 3; m++) { + Y = bessely(m, z); + printf("Order m = %d, z = %f, Y = %20.17e, Y_true = %20.17e\n", m, z, Y, Y_true[m + 3]); + printf("Checking, abs tol = %e ... ", tol); + if (abs(Y - Y_true[m + 3]) < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed!\n"); + } + } + } + + printf("==============================================\n"); + printf("Test Mathieu modmc1\n"); + double mc1; + double mc1d; + N = 6; + { + double q = 0.01; + double v = 4.0d; + double tol = 1e-12; + double rtol = 1e-11; + // Golden values obtained using Matlab. I get values + // for m = 1, 4, 7, 10, 13, 16. + // I use widely spaced orders to verify adaption works OK. + double mc1_true[6] = {-3.43416995264759106e-01, 3.97873130174148326e-01, 8.35567375774955712e-02, + 3.15168992451903213e-03, 4.36682920825898135e-05, 2.91722697679619490e-07}; + double mc1d_true[6] = {2.32628352485249795e-01, -1.37343635827296445e-01, 4.09262446440566696e-01, + 2.69766257109354118e-02, 5.19428246083634753e-04, 4.40526164453159552e-06}; + + printf("q = %f\n", q); + for (int i = 0; i < N; i++) { + int m = 3 * i + 1; + retcode = mathieu_modmc1(m, q, v, &mc1, &mc1d); + printf( + "Iteration i = %d, order m = %d, v = %f, mc1 = %20.17e, mc1d = %20.17e\n", i, m, v, mc1, mc1d + ); + printf( + " mc1_true = %20.17e, mc1d_true = %20.17e\n", mc1_true[i], + mc1d_true[i] + ); + // For small 1 modmc1's output is very small. To adequately test, I check + // both the absolute err as well as the rel error. I get one failure + // which suggests some minute difference between Matlab and C imps. + // The difference occurs because of catastrophic cancellation, but + // I haven't found why Matlab and C are different yet. + printf("Checking abs err, tol = %e ... ", tol); + double diff = abs((mc1 - mc1_true[i])); + double diffd = abs((mc1d - mc1d_true[i])); + // printf("diff = %e, diffd = %e\n", diff, diffd); + if (diff < tol && diffd < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! tol = %e, diff = %e, diffd = %e\n", tol, diff, diffd); + } + printf("Checking rel err, rtol = %e... ", rtol); + double rdiff = abs((mc1 - mc1_true[i]) / mc1); + double rdiffd = abs((mc1d - mc1d_true[i]) / mc1d); + // printf("rdiff = %e, rdiffd = %e\n", rdiff, rdiffd); + if (rdiff < rtol && rdiffd < rtol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! rtol = %e, rdiff = %e, rdiffd = %e\n", rtol, rdiff, rdiffd); + } + printf("--------------------------------\n"); + } + } + printf("----------------------------------------------\n"); + { + double q = 100; + double v = 4.0d; + double tol = 1e-14; + double rtol = 1e-13; + // Golden values obtained using Matlab. I get values + // for m = 1, 4, 7, 10, 13, 16. + // I use widely spaced orders to verify adaption works OK. + double mc1_true[6] = {-3.41428322016207667e-02, 3.44948654167566600e-03, 3.35219048366175096e-02, + -9.08021166319475555e-03, -3.24223338845559816e-02, 1.27646351014092524e-02}; + double mc1d_true[6] = {2.00691110743774284e-02, 1.85469526520107095e+01, -3.57255502109930578e+00, + -1.79656659177568585e+01, 5.87447048444966935e+00, 1.72822185411075182e+01}; + + printf("q = %f\n", q); + for (int i = 0; i < N; i++) { + int m = 3 * i + 1; + retcode = mathieu_modmc1(m, q, v, &mc1, &mc1d); + printf( + "Iteration i = %d, order m = %d, v = %f, mc1 = %20.17e, mc1d = %20.17e\n", i, m, v, mc1, mc1d + ); + printf( + " mc1_true = %20.17e, mc1d_true = %20.17e\n", mc1_true[i], + mc1d_true[i] + ); + printf("Checking abs err, tol = %e ... ", tol); + double diff = abs((mc1 - mc1_true[i])); + double diffd = abs((mc1d - mc1d_true[i])); + // printf("diff = %e, diffd = %e\n", diff, diffd); + if (diff < tol && diffd < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! tol = %e, diff = %e, diffd = %e\n", tol, diff, diffd); + } + printf("Checking rel err, rtol = %e... ", rtol); + double rdiff = abs((mc1 - mc1_true[i]) / mc1); + double rdiffd = abs((mc1d - mc1d_true[i]) / mc1d); + // printf("rdiff = %e, rdiffd = %e\n", rdiff, rdiffd); + if (rdiff < rtol && rdiffd < rtol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! rtol = %e, rdiff = %e, rdiffd = %e\n", rtol, rdiff, rdiffd); + } + + printf("--------------------------------\n"); + } + } + + printf("==============================================\n"); + printf("Test Mathieu modms1\n"); + double ms1; + double ms1d; + N = 6; + { + double q = 0.01; + double v = 4.0d; + double tol = 1e-6; + double rtol = 1e-11; // I use different tol for relative err test. + // Golden values obtained using Matlab. I used + // m = 1, 4, 7, 10, 13, 16 + double ms1_true[6] = {-3.43379253985569288e-01, 3.97873130174058565e-01, 8.35567375774849824e-02, + 3.15168992462108635e-03, 4.36682920825898135e-05, 2.91722697679619490e-07}; + double ms1d_true[6] = {2.29156633917084712e-01, -1.37343635825750598e-01, 4.09262446440550431e-01, + 2.69766257110715633e-02, 5.19428246083634753e-04, 4.40526164453159637e-06}; + + printf("q = %f\n", q); + for (int i = 0; i < N; i++) { + int m = 3 * i + 1; + retcode = mathieu_modms1(m, q, v, &ms1, &ms1d); + printf( + "Iteration i = %d, order m = %d, v = %f, ms1 = %20.17e, ms1d = %20.17e\n", i, m, v, ms1, ms1d + ); + printf( + " ms1_true = %20.17e, ms1d_true = %20.17e\n", ms1_true[i], + ms1d_true[i] + ); + // For small 1 modmc1's output is very small. To adequately test, I check + // both the absolute err as well as the rel error. I get one failure + // which suggests some minute difference between Matlab and C imps. + // The difference occurs because of catastrophic cancellation, but + // I haven't found why Matlab and C are different yet. + printf("Checking abs err, tol = %e ... ", tol); + double diff = abs((ms1 - ms1_true[i])); + double diffd = abs((ms1d - ms1d_true[i])); + // printf("diff = %e, diffd = %e\n", diff, diffd); + if (diff < tol && diffd < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! tol = %e, diff = %e, diffd = %e\n", tol, diff, diffd); + } + printf("Checking rel err, rtol = %e... ", rtol); + double rdiff = abs((ms1 - ms1_true[i]) / ms1); + double rdiffd = abs((ms1d - ms1d_true[i]) / ms1d); + // printf("rdiff = %e, rdiffd = %e\n", rdiff, rdiffd); + if (rdiff < rtol && rdiffd < rtol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! rtol = %e, rdiff = %e, rdiffd = %e\n", rtol, rdiff, rdiffd); + } + printf("--------------------------------\n"); + } + } + printf("----------------------------------------------\n"); + { + double q = 100; + double v = 4.0d; + double tol = 1e-14; + double rtol = 1e-14; + // Golden values obtained using Matlab. I used + // m = 1, 4, 7, 10, 13, 16 + double ms1_true[6] = {-3.41201739333336501e-02, 2.33952744720816123e-03, 3.36955475281936756e-02, + -8.29453993546177236e-03, -3.24933710575481122e-02, 1.27631677399401809e-02}; + double ms1d_true[6] = {-6.45251275907079758e-01, 1.85996271072809449e+01, -3.04135814724805620e+00, + -1.80792391857584605e+01, 5.75576271594593436e+00, 1.72825425245382398e+01}; + + printf("q = %f\n", q); + for (int i = 0; i < N; i++) { + int m = 3 * i + 1; + retcode = mathieu_modms1(m, q, v, &ms1, &ms1d); + printf( + "Iteration i = %d, order m = %d, v = %f, ms1 = %20.17e, ms1d = %20.17e\n", i, m, v, ms1, ms1d + ); + printf( + " ms1_true = %20.17e, ms1d_true = %20.17e\n", ms1_true[i], + ms1d_true[i] + ); + printf("Checking abs err, tol = %e ... ", tol); + double diff = abs((ms1 - ms1_true[i])); + double diffd = abs((ms1d - ms1d_true[i])); + // printf("diff = %e, diffd = %e\n", diff, diffd); + if (diff < tol && diffd < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! tol = %e, diff = %e, diffd = %e\n", tol, diff, diffd); + } + printf("Checking rel err, rtol = %e... ", rtol); + double rdiff = abs((ms1 - ms1_true[i]) / ms1); + double rdiffd = abs((ms1d - ms1d_true[i]) / ms1d); + // printf("rdiff = %e, rdiffd = %e\n", rdiff, rdiffd); + if (rdiff < rtol && rdiffd < rtol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! rtol = %e, rdiff = %e, rdiffd = %e\n", rtol, rdiff, rdiffd); + } + printf("--------------------------------\n"); + } + } + + printf("==============================================\n"); + printf("Test Mathieu modmc2\n"); + double mc2; + double mc2d; + N = 6; + { + double v = 4.0d; // Test at v=4 since the value is extremely large for + // smaller v. + double q = 0.01; + double tol = 1e-9; + double rtol = 1e-12; // I use different tol for relative err test. + + // Golden values obtained using Matlab. I get values + // for m = 1, 4, 7, 10, 13, 16. + // I use widely spaced orders to verify adaption works OK. + double mc2_true[6] = {-1.05364049460497881e-02, -5.43273917516489124e-02, -8.98079478272420850e-01, + -1.21204986957549021e+01, -6.18466504094926336e+02, -7.25758230510118883e+04}; + double mc2d_true[6] = {-1.84664333620107901e+00, 1.61881073394890906e+00, 3.22020193455027437e+00, + 9.82487563907863120e+01, 7.22196427512374794e+03, 1.08632040247457405e+06}; + + printf("q = %f\n", q); + for (int i = 0; i < N; i++) { + int m = 3 * i + 1; + retcode = mathieu_modmc2(m, q, v, &mc2, &mc2d); + printf( + "Iteration i = %d, order m = %d, v = %f, mc2 = %20.17e, mc2d = %20.17e\n", i, m, v, mc2, mc2d + ); + printf( + " mc2_true = %20.17e, mc2d_true = %20.17e\n", mc2_true[i], + mc2d_true[i] + ); + // This fcns output varies widely. To adequately test, I check + // both the absolute err as well as the rel error. + printf("Checking abs err, tol = %e ... ", tol); + double diff = abs((mc2 - mc2_true[i])); + double diffd = abs((mc2d - mc2d_true[i])); + // printf("diff = %e, diffd = %e\n", diff, diffd); + if (diff < tol && diffd < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! tol = %e, diff = %e, diffd = %e\n", tol, diff, diffd); + } + printf("Checking rel err, rtol = %e... ", rtol); + double rdiff = abs((mc2 - mc2_true[i]) / mc2); + double rdiffd = abs((mc2d - mc2d_true[i]) / mc2d); + // printf("rdiff = %e, rdiffd = %e\n", rdiff, rdiffd); + if (rdiff < rtol && rdiffd < rtol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! rtol = %e, rdiff = %e, rdiffd = %e\n", rtol, rdiff, rdiffd); + } + printf("--------------------------------\n"); + } + } + printf("----------------------------------------------\n"); + { + double q = 100; + double v = 4.0d; // Test at v=4 since the value is extremely large for smaller v. + double tol = 1e-14; + double rtol = 1e-12; + // Golden values obtained using Matlab. I get values + // for m = 1, 4, 7, 10, 13, 16. + // I use widely spaced orders to verify adaption works OK. + double mc2_true[6] = {-5.50400810104661117e-06, -3.39713077181867523e-02, 6.51339805952160303e-03, + 3.29219862870713670e-02, -1.07334892119392861e-02, -3.16798833654585860e-02}; + double mc2d_true[6] = {-1.86457777769476074e+01, 1.90043837315601305e+00, 1.82970001977819763e+01, + -4.97283180832788130e+00, -1.76904663568990372e+01, 6.98187640938441945e+00}; + + printf("q = %f\n", q); + for (int i = 0; i < N; i++) { + int m = 3 * i + 1; + retcode = mathieu_modmc2(m, q, v, &mc2, &mc2d); + printf( + "Iteration i = %d, order m = %d, v = %f, mc2 = %20.17e, mc2d = %20.17e\n", i, m, v, mc2, mc2d + ); + printf( + " mc2_true = %20.17e, mc2d_true = %20.17e\n", mc2_true[i], + mc2d_true[i] + ); + printf("Checking abs err, tol = %e ... ", tol); + double diff = abs((mc2 - mc2_true[i])); + double diffd = abs((mc2d - mc2d_true[i])); + // printf("diff = %e, diffd = %e\n", diff, diffd); + if (diff < tol && diffd < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! diff = %e, diffd = %e\n", diff, diffd); + } + printf("Checking rel err, rtol = %e... ", rtol); + double rdiff = abs((mc2 - mc2_true[i]) / mc2); + double rdiffd = abs((mc2d - mc2d_true[i]) / mc2d); + // printf("rdiff = %e, rdiffd = %e\n", rdiff, rdiffd); + if (rdiff < rtol && rdiffd < rtol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! rdiff = %e, rdiffd = %e\n", rdiff, rdiffd); + } + printf("--------------------------------\n"); + } + } + + printf("==============================================\n"); + printf("Test Mathieu modms2\n"); + double ms2; + double ms2d; + N = 6; + { + double q = 0.01; + double v = 4.0d; // Test at v=4 since the value is extremely large for smaller v. + double tol = 1e-2; + double rtol = 1e-8; + // Golden values obtained using Matlab. I used + // m = 1, 4, 7, 10, 13, 16 + double ms2_true[6] = { + -9.91337437572497975e-03, -5.43273917519763866e-02, -8.98079478272391984e-01, + -1.21204986957547884e+01, -6.18466504094926336e+02, -7.25758230510119029e+04, + }; + double ms2d_true[6] = { + -1.84736861502861638e+00, 1.61881073394918618e+00, 3.22020193455016246e+00, + 9.82487563907851893e+01, 7.22196427512374976e+03, 1.08632040247457428e+06, + }; + + printf("q = %f\n", q); + for (int i = 0; i < N; i++) { + int m = 3 * i + 1; + retcode = mathieu_modms2(m, q, v, &ms2, &ms2d); + printf( + "Iteration i = %d, order m = %d, v = %f, ms2 = %20.17e, ms2d = %20.17e\n", i, m, v, ms2, ms2d + ); + printf( + " ms2_true = %20.17e, ms2d_true = %20.17e\n", ms2_true[i], + ms2d_true[i] + ); + // To adequately test, I check + // both the absolute err as well as the rel error. + printf("Checking abs err, tol = %e ... ", tol); + double diff = abs((ms2 - ms2_true[i])); + double diffd = abs((ms2d - ms2d_true[i])); + // printf("diff = %e, diffd = %e\n", diff, diffd); + if (diff < tol && diffd < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! tol = %e, diff = %e, diffd = %e\n", tol, diff, diffd); + } + printf("Checking rel err, rtol = %e... ", rtol); + double rdiff = abs((ms2 - ms2_true[i]) / ms2); + double rdiffd = abs((ms2d - ms2d_true[i]) / ms2d); + // printf("rdiff = %e, rdiffd = %e\n", rdiff, rdiffd); + if (rdiff < rtol && rdiffd < rtol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! rtol = %e, rdiff = %e, rdiffd = %e\n", rtol, rdiff, rdiffd); + } + printf("--------------------------------\n"); + } + } + printf("----------------------------------------------\n"); + { + double q = 100; + double v = 4.0d; // Test at v=4 since the value is extremely large for smaller v. + double tol = 1e-14; + double rtol = 1e-14; + // Golden values obtained using Matlab. I used + // m = 1, 4, 7, 10, 13, 16 + double ms2_true[6] = {1.21267949732755992e-03, -3.40647275296700053e-02, 5.53991666556525025e-03, + 3.31278916985982691e-02, -1.05157990917375058e-02, -3.16804730008901705e-02}; + double ms2d_true[6] = {-1.86352300142727003e+01, 1.29451048253596168e+00, 1.83932580755102961e+01, + -4.54427791318153673e+00, -1.77295648089783491e+01, 6.98107651232898974e+00}; + + printf("q = %f\n", q); + for (int i = 0; i < N; i++) { + int m = 3 * i + 1; + retcode = mathieu_modms2(m, q, v, &ms2, &ms2d); + printf( + "Iteration i = %d, order m = %d, v = %f, ms2 = %20.17e, ms2d = %20.17e\n", i, m, v, ms2, ms2d + ); + printf( + " ms2_true = %20.17e, ms2d_true = %20.17e\n", ms2_true[i], + ms2d_true[i] + ); + printf("Checking abs err, tol = %e ... ", tol); + double diff = abs((ms2 - ms2_true[i])); + double diffd = abs((ms2d - ms2d_true[i])); + // printf("diff = %e, diffd = %e\n", diff, diffd); + if (diff < tol && diffd < tol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! tol = %e, diff = %e, diffd = %e\n", tol, diff, diffd); + } + printf("Checking rel err, rtol = %e... ", rtol); + double rdiff = abs((ms2 - ms2_true[i]) / ms2); + double rdiffd = abs((ms2d - ms2d_true[i]) / ms2d); + // printf("rdiff = %e, rdiffd = %e\n", rdiff, rdiffd); + if (rdiff < rtol && rdiffd < rtol) { + pass++; + printf("Passed!\n"); + } else { + fail++; + printf("Failed! rtol = %e, rdiff = %e, rdiffd = %e\n", rtol, rdiff, rdiffd); + } + printf("--------------------------------\n"); + } + } + + printf("===========================================\n"); + printf("At end of run, pass = %d, fail = %d\n", pass, fail); + + return 0; + } + +} // namespace mathieu +} // namespace xsf