summaryrefslogtreecommitdiff
path: root/gsl-1.9/specfunc/mathieu_radfunc.c
diff options
context:
space:
mode:
Diffstat (limited to 'gsl-1.9/specfunc/mathieu_radfunc.c')
-rw-r--r--gsl-1.9/specfunc/mathieu_radfunc.c459
1 files changed, 459 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/mathieu_radfunc.c b/gsl-1.9/specfunc/mathieu_radfunc.c
new file mode 100644
index 0000000..844a0f4
--- /dev/null
+++ b/gsl-1.9/specfunc/mathieu_radfunc.c
@@ -0,0 +1,459 @@
+/* specfunc/mathieu_radfunc.c
+ *
+ * Copyright (C) 2002 Lowell Johnson
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ */
+
+/* Author: L. Johnson */
+
+#include <config.h>
+#include <stdlib.h>
+#include <math.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_sf_mathieu.h>
+
+
+int gsl_sf_mathieu_Mc(int kind, int order, double qq, double zz,
+ gsl_sf_result *result)
+{
+ int even_odd, kk, mm, status;
+ double maxerr = 1e-14, amax, pi = M_PI, fn, factor;
+ double coeff[GSL_SF_MATHIEU_COEFF], fc;
+ double j1c, z2c, j1pc, z2pc;
+ double u1, u2;
+ gsl_sf_result aa;
+
+
+ /* Check for out of bounds parameters. */
+ if (qq <= 0.0)
+ {
+ GSL_ERROR("q must be greater than zero", GSL_EINVAL);
+ }
+ if (kind < 1 || kind > 2)
+ {
+ GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
+ }
+
+ mm = 0;
+ amax = 0.0;
+ fn = 0.0;
+ u1 = sqrt(qq)*exp(-1.0*zz);
+ u2 = sqrt(qq)*exp(zz);
+
+ even_odd = 0;
+ if (order % 2 != 0)
+ even_odd = 1;
+
+ /* Compute the characteristic value. */
+ status = gsl_sf_mathieu_a(order, qq, &aa);
+ if (status != GSL_SUCCESS)
+ {
+ return status;
+ }
+
+ /* Compute the series coefficients. */
+ status = gsl_sf_mathieu_a_coeff(order, qq, aa.val, coeff);
+ if (status != GSL_SUCCESS)
+ {
+ return status;
+ }
+
+ if (even_odd == 0)
+ {
+ for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
+ {
+ amax = GSL_MAX(amax, fabs(coeff[kk]));
+ if (fabs(coeff[kk])/amax < maxerr)
+ break;
+
+ j1c = gsl_sf_bessel_Jn(kk, u1);
+ if (kind == 1)
+ {
+ z2c = gsl_sf_bessel_Jn(kk, u2);
+ }
+ else /* kind = 2 */
+ {
+ z2c = gsl_sf_bessel_Yn(kk, u2);
+ }
+
+ fc = pow(-1.0, 0.5*order+kk)*coeff[kk];
+ fn += fc*j1c*z2c;
+ }
+
+ fn *= sqrt(pi/2.0)/coeff[0];
+ }
+ else
+ {
+ for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
+ {
+ amax = GSL_MAX(amax, fabs(coeff[kk]));
+ if (fabs(coeff[kk])/amax < maxerr)
+ break;
+
+ j1c = gsl_sf_bessel_Jn(kk, u1);
+ j1pc = gsl_sf_bessel_Jn(kk+1, u1);
+ if (kind == 1)
+ {
+ z2c = gsl_sf_bessel_Jn(kk, u2);
+ z2pc = gsl_sf_bessel_Jn(kk+1, u2);
+ }
+ else /* kind = 2 */
+ {
+ z2c = gsl_sf_bessel_Yn(kk, u2);
+ z2pc = gsl_sf_bessel_Yn(kk+1, u2);
+ }
+ fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
+ fn += fc*(j1c*z2pc + j1pc*z2c);
+ }
+
+ fn *= sqrt(pi/2.0)/coeff[0];
+ }
+
+ result->val = fn;
+ result->err = 2.0*GSL_DBL_EPSILON;
+ factor = fabs(fn);
+ if (factor > 1.0)
+ result->err *= factor;
+
+ return GSL_SUCCESS;
+}
+
+
+int gsl_sf_mathieu_Ms(int kind, int order, double qq, double zz,
+ gsl_sf_result *result)
+{
+ int even_odd, kk, mm, status;
+ double maxerr = 1e-14, amax, pi = M_PI, fn, factor;
+ double coeff[GSL_SF_MATHIEU_COEFF], fc;
+ double j1c, z2c, j1mc, z2mc, j1pc, z2pc;
+ double u1, u2;
+ gsl_sf_result aa;
+
+
+ /* Check for out of bounds parameters. */
+ if (qq <= 0.0)
+ {
+ GSL_ERROR("q must be greater than zero", GSL_EINVAL);
+ }
+ if (kind < 1 || kind > 2)
+ {
+ GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
+ }
+
+ mm = 0;
+ amax = 0.0;
+ fn = 0.0;
+ u1 = sqrt(qq)*exp(-1.0*zz);
+ u2 = sqrt(qq)*exp(zz);
+
+ even_odd = 0;
+ if (order % 2 != 0)
+ even_odd = 1;
+
+ /* Compute the characteristic value. */
+ status = gsl_sf_mathieu_b(order, qq, &aa);
+ if (status != GSL_SUCCESS)
+ {
+ return status;
+ }
+
+ /* Compute the series coefficients. */
+ status = gsl_sf_mathieu_b_coeff(order, qq, aa.val, coeff);
+ if (status != GSL_SUCCESS)
+ {
+ return status;
+ }
+
+ if (even_odd == 0)
+ {
+ for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
+ {
+ amax = GSL_MAX(amax, fabs(coeff[kk]));
+ if (fabs(coeff[kk])/amax < maxerr)
+ break;
+
+ j1mc = gsl_sf_bessel_Jn(kk, u1);
+ j1pc = gsl_sf_bessel_Jn(kk+2, u1);
+ if (kind == 1)
+ {
+ z2mc = gsl_sf_bessel_Jn(kk, u2);
+ z2pc = gsl_sf_bessel_Jn(kk+2, u2);
+ }
+ else /* kind = 2 */
+ {
+ z2mc = gsl_sf_bessel_Yn(kk, u2);
+ z2pc = gsl_sf_bessel_Yn(kk+2, u2);
+ }
+
+ fc = pow(-1.0, 0.5*order+kk+1)*coeff[kk];
+ fn += fc*(j1mc*z2pc - j1pc*z2mc);
+ }
+
+ fn *= sqrt(pi/2.0)/coeff[0];
+ }
+ else
+ {
+ for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
+ {
+ amax = GSL_MAX(amax, fabs(coeff[kk]));
+ if (fabs(coeff[kk])/amax < maxerr)
+ break;
+
+ j1c = gsl_sf_bessel_Jn(kk, u1);
+ j1pc = gsl_sf_bessel_Jn(kk+1, u1);
+ if (kind == 1)
+ {
+ z2c = gsl_sf_bessel_Jn(kk, u2);
+ z2pc = gsl_sf_bessel_Jn(kk+1, u2);
+ }
+ else /* kind = 2 */
+ {
+ z2c = gsl_sf_bessel_Yn(kk, u2);
+ z2pc = gsl_sf_bessel_Yn(kk+1, u2);
+ }
+
+ fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
+ fn += fc*(j1c*z2pc - j1pc*z2c);
+ }
+
+ fn *= sqrt(pi/2.0)/coeff[0];
+ }
+
+ result->val = fn;
+ result->err = 2.0*GSL_DBL_EPSILON;
+ factor = fabs(fn);
+ if (factor > 1.0)
+ result->err *= factor;
+
+ return GSL_SUCCESS;
+}
+
+
+int gsl_sf_mathieu_Mc_array(int kind, int nmin, int nmax, double qq,
+ double zz, gsl_sf_mathieu_workspace *work,
+ double result_array[])
+{
+ int even_odd, order, ii, kk, mm, status;
+ double maxerr = 1e-14, amax, pi = M_PI, fn;
+ double coeff[GSL_SF_MATHIEU_COEFF], fc;
+ double j1c, z2c, j1pc, z2pc;
+ double u1, u2;
+ double *aa = work->aa;
+
+
+ /* Initialize the result array to zeroes. */
+ for (ii=0; ii<nmax-nmin+1; ii++)
+ result_array[ii] = 0.0;
+
+ /* Check for out of bounds parameters. */
+ if (qq <= 0.0)
+ {
+ GSL_ERROR("q must be greater than zero", GSL_EINVAL);
+ }
+ if (kind < 1 || kind > 2)
+ {
+ GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
+ }
+
+ mm = 0;
+ amax = 0.0;
+ fn = 0.0;
+ u1 = sqrt(qq)*exp(-1.0*zz);
+ u2 = sqrt(qq)*exp(zz);
+
+ /* Compute all eigenvalues up to nmax. */
+ gsl_sf_mathieu_a_array(0, nmax, qq, work, aa);
+
+ for (ii=0, order=nmin; order<=nmax; ii++, order++)
+ {
+ even_odd = 0;
+ if (order % 2 != 0)
+ even_odd = 1;
+
+ /* Compute the series coefficients. */
+ status = gsl_sf_mathieu_a_coeff(order, qq, aa[order], coeff);
+ if (status != GSL_SUCCESS)
+ {
+ return status;
+ }
+
+ if (even_odd == 0)
+ {
+ for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
+ {
+ amax = GSL_MAX(amax, fabs(coeff[kk]));
+ if (fabs(coeff[kk])/amax < maxerr)
+ break;
+
+ j1c = gsl_sf_bessel_Jn(kk, u1);
+ if (kind == 1)
+ {
+ z2c = gsl_sf_bessel_Jn(kk, u2);
+ }
+ else /* kind = 2 */
+ {
+ z2c = gsl_sf_bessel_Yn(kk, u2);
+ }
+
+ fc = pow(-1.0, 0.5*order+kk)*coeff[kk];
+ fn += fc*j1c*z2c;
+ }
+
+ fn *= sqrt(pi/2.0)/coeff[0];
+ }
+ else
+ {
+ for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
+ {
+ amax = GSL_MAX(amax, fabs(coeff[kk]));
+ if (fabs(coeff[kk])/amax < maxerr)
+ break;
+
+ j1c = gsl_sf_bessel_Jn(kk, u1);
+ j1pc = gsl_sf_bessel_Jn(kk+1, u1);
+ if (kind == 1)
+ {
+ z2c = gsl_sf_bessel_Jn(kk, u2);
+ z2pc = gsl_sf_bessel_Jn(kk+1, u2);
+ }
+ else /* kind = 2 */
+ {
+ z2c = gsl_sf_bessel_Yn(kk, u2);
+ z2pc = gsl_sf_bessel_Yn(kk+1, u2);
+ }
+ fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
+ fn += fc*(j1c*z2pc + j1pc*z2c);
+ }
+
+ fn *= sqrt(pi/2.0)/coeff[0];
+ }
+
+ result_array[ii] = fn;
+ } /* order loop */
+
+ return GSL_SUCCESS;
+}
+
+
+int gsl_sf_mathieu_Ms_array(int kind, int nmin, int nmax, double qq,
+ double zz, gsl_sf_mathieu_workspace *work,
+ double result_array[])
+{
+ int even_odd, order, ii, kk, mm, status;
+ double maxerr = 1e-14, amax, pi = M_PI, fn;
+ double coeff[GSL_SF_MATHIEU_COEFF], fc;
+ double j1c, z2c, j1mc, z2mc, j1pc, z2pc;
+ double u1, u2;
+ double *bb = work->bb;
+
+
+ /* Initialize the result array to zeroes. */
+ for (ii=0; ii<nmax-nmin+1; ii++)
+ result_array[ii] = 0.0;
+
+ /* Check for out of bounds parameters. */
+ if (qq <= 0.0)
+ {
+ GSL_ERROR("q must be greater than zero", GSL_EINVAL);
+ }
+ if (kind < 1 || kind > 2)
+ {
+ GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
+ }
+
+ mm = 0;
+ amax = 0.0;
+ fn = 0.0;
+ u1 = sqrt(qq)*exp(-1.0*zz);
+ u2 = sqrt(qq)*exp(zz);
+
+ /* Compute all eigenvalues up to nmax. */
+ gsl_sf_mathieu_b_array(0, nmax, qq, work, bb);
+
+ for (ii=0, order=nmin; order<=nmax; ii++, order++)
+ {
+ even_odd = 0;
+ if (order % 2 != 0)
+ even_odd = 1;
+
+ /* Compute the series coefficients. */
+ status = gsl_sf_mathieu_b_coeff(order, qq, bb[order], coeff);
+ if (status != GSL_SUCCESS)
+ {
+ return status;
+ }
+
+ if (even_odd == 0)
+ {
+ for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
+ {
+ amax = GSL_MAX(amax, fabs(coeff[kk]));
+ if (fabs(coeff[kk])/amax < maxerr)
+ break;
+
+ j1mc = gsl_sf_bessel_Jn(kk, u1);
+ j1pc = gsl_sf_bessel_Jn(kk+2, u1);
+ if (kind == 1)
+ {
+ z2mc = gsl_sf_bessel_Jn(kk, u2);
+ z2pc = gsl_sf_bessel_Jn(kk+2, u2);
+ }
+ else /* kind = 2 */
+ {
+ z2mc = gsl_sf_bessel_Yn(kk, u2);
+ z2pc = gsl_sf_bessel_Yn(kk+2, u2);
+ }
+
+ fc = pow(-1.0, 0.5*order+kk+1)*coeff[kk];
+ fn += fc*(j1mc*z2pc - j1pc*z2mc);
+ }
+
+ fn *= sqrt(pi/2.0)/coeff[0];
+ }
+ else
+ {
+ for (kk=0; kk<GSL_SF_MATHIEU_COEFF; kk++)
+ {
+ amax = GSL_MAX(amax, fabs(coeff[kk]));
+ if (fabs(coeff[kk])/amax < maxerr)
+ break;
+
+ j1c = gsl_sf_bessel_Jn(kk, u1);
+ j1pc = gsl_sf_bessel_Jn(kk+1, u1);
+ if (kind == 1)
+ {
+ z2c = gsl_sf_bessel_Jn(kk, u2);
+ z2pc = gsl_sf_bessel_Jn(kk+1, u2);
+ }
+ else /* kind = 2 */
+ {
+ z2c = gsl_sf_bessel_Yn(kk, u2);
+ z2pc = gsl_sf_bessel_Yn(kk+1, u2);
+ }
+
+ fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
+ fn += fc*(j1c*z2pc - j1pc*z2c);
+ }
+
+ fn *= sqrt(pi/2.0)/coeff[0];
+ }
+
+ result_array[ii] = fn;
+ } /* order loop */
+
+ return GSL_SUCCESS;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-30 22:20:48 +0000