summaryrefslogtreecommitdiff
path: root/gsl-1.9/specfunc/mathieu_angfunc.c
diff options
context:
space:
mode:
Diffstat (limited to 'gsl-1.9/specfunc/mathieu_angfunc.c')
-rw-r--r--gsl-1.9/specfunc/mathieu_angfunc.c343
1 files changed, 343 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/mathieu_angfunc.c b/gsl-1.9/specfunc/mathieu_angfunc.c
new file mode 100644
index 0000000..e572c22
--- /dev/null
+++ b/gsl-1.9/specfunc/mathieu_angfunc.c
@@ -0,0 +1,343 @@
+/* specfunc/mathieu_angfunc.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 <stdio.h>
+#include <math.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_sf_mathieu.h>
+
+
+int gsl_sf_mathieu_ce(int order, double qq, double zz, gsl_sf_result *result)
+{
+ int even_odd, ii, status;
+ double coeff[GSL_SF_MATHIEU_COEFF], norm, fn, factor;
+ gsl_sf_result aa;
+
+
+ norm = 0.0;
+ even_odd = 0;
+ if (order % 2 != 0)
+ even_odd = 1;
+
+ /* Handle the trivial case where q = 0. */
+ if (qq == 0.0)
+ {
+ norm = 1.0;
+ if (order == 0)
+ norm = sqrt(2.0);
+
+ fn = cos(order*zz)/norm;
+
+ result->val = fn;
+ result->err = 2.0*GSL_DBL_EPSILON;
+ factor = fabs(fn);
+ if (factor > 1.0)
+ result->err *= factor;
+
+ return GSL_SUCCESS;
+ }
+
+ /* Use symmetry characteristics of the functions to handle cases with
+ negative order. */
+ if (order < 0)
+ order *= -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)
+ {
+ fn = 0.0;
+ norm = coeff[0]*coeff[0];
+ for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
+ {
+ fn += coeff[ii]*cos(2.0*ii*zz);
+ norm += coeff[ii]*coeff[ii];
+ }
+ }
+ else
+ {
+ fn = 0.0;
+ for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
+ {
+ fn += coeff[ii]*cos((2.0*ii + 1.0)*zz);
+ norm += coeff[ii]*coeff[ii];
+ }
+ }
+
+ norm = sqrt(norm);
+ fn /= norm;
+
+ 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_se(int order, double qq, double zz, gsl_sf_result *result)
+{
+ int even_odd, ii, status;
+ double coeff[GSL_SF_MATHIEU_COEFF], norm, fn, factor;
+ gsl_sf_result aa;
+
+
+ norm = 0.0;
+ even_odd = 0;
+ if (order % 2 != 0)
+ even_odd = 1;
+
+ /* Handle the trivial cases where order = 0 and/or q = 0. */
+ if (order == 0)
+ {
+ result->val = 0.0;
+ result->err = 0.0;
+ return GSL_SUCCESS;
+ }
+
+ if (qq == 0.0)
+ {
+ norm = 1.0;
+ fn = sin(order*zz);
+
+ result->val = fn;
+ result->err = 2.0*GSL_DBL_EPSILON;
+ factor = fabs(fn);
+ if (factor > 1.0)
+ result->err *= factor;
+
+ return GSL_SUCCESS;
+ }
+
+ /* Use symmetry characteristics of the functions to handle cases with
+ negative order. */
+ if (order < 0)
+ order *= -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)
+ {
+ fn = 0.0;
+ for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
+ {
+ norm += coeff[ii]*coeff[ii];
+ fn += coeff[ii]*sin(2.0*(ii + 1)*zz);
+ }
+ }
+ else
+ {
+ fn = 0.0;
+ for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
+ {
+ norm += coeff[ii]*coeff[ii];
+ fn += coeff[ii]*sin((2.0*ii + 1)*zz);
+ }
+ }
+ norm = sqrt(norm);
+ fn /= norm;
+
+ 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_ce_array(int nmin, int nmax, double qq, double zz,
+ gsl_sf_mathieu_workspace *work,
+ double result_array[])
+{
+ int even_odd, order, ii, jj, status;
+ double coeff[GSL_SF_MATHIEU_COEFF], *aa = work->aa, norm;
+
+
+ /* Initialize the result array to zeroes. */
+ for (ii=0; ii<nmax-nmin+1; ii++)
+ result_array[ii] = 0.0;
+
+ /* Ensure that the workspace is large enough to accomodate. */
+ if (work->size < (unsigned int)nmax)
+ {
+ GSL_ERROR("Work space not large enough", GSL_EINVAL);
+ }
+
+ if (nmin < 0 || nmax < nmin)
+ {
+ GSL_ERROR("domain error", GSL_EDOM);
+ }
+
+ /* Compute all of the eigenvalues up to nmax. */
+ gsl_sf_mathieu_a_array(0, nmax, qq, work, aa);
+
+ for (ii=0, order=nmin; order<=nmax; ii++, order++)
+ {
+ norm = 0.0;
+ even_odd = 0;
+ if (order % 2 != 0)
+ even_odd = 1;
+
+ /* Handle the trivial case where q = 0. */
+ if (qq == 0.0)
+ {
+ norm = 1.0;
+ if (order == 0)
+ norm = sqrt(2.0);
+
+ result_array[ii] = cos(order*zz)/norm;
+
+ continue;
+ }
+
+ /* 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)
+ {
+ norm = coeff[0]*coeff[0];
+ for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
+ {
+ result_array[ii] += coeff[jj]*cos(2.0*jj*zz);
+ norm += coeff[jj]*coeff[jj];
+ }
+ }
+ else
+ {
+ for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
+ {
+ result_array[ii] += coeff[jj]*cos((2.0*jj + 1.0)*zz);
+ norm += coeff[jj]*coeff[jj];
+ }
+ }
+
+ norm = sqrt(norm);
+ result_array[ii] /= norm;
+ }
+
+ return GSL_SUCCESS;
+}
+
+
+int gsl_sf_mathieu_se_array(int nmin, int nmax, double qq, double zz,
+ gsl_sf_mathieu_workspace *work,
+ double result_array[])
+{
+ int even_odd, order, ii, jj, status;
+ double coeff[GSL_SF_MATHIEU_COEFF], *bb = work->bb, norm;
+
+
+ /* Initialize the result array to zeroes. */
+ for (ii=0; ii<nmax-nmin+1; ii++)
+ result_array[ii] = 0.0;
+
+ /* Ensure that the workspace is large enough to accomodate. */
+ if (work->size < (unsigned int)nmax)
+ {
+ GSL_ERROR("Work space not large enough", GSL_EINVAL);
+ }
+
+ if (nmin < 0 || nmax < nmin)
+ {
+ GSL_ERROR("domain error", GSL_EDOM);
+ }
+
+ /* Compute all of the eigenvalues up to nmax. */
+ gsl_sf_mathieu_b_array(0, nmax, qq, work, bb);
+
+ for (ii=0, order=nmin; order<=nmax; ii++, order++)
+ {
+ norm = 0.0;
+ even_odd = 0;
+ if (order % 2 != 0)
+ even_odd = 1;
+
+ /* Handle the trivial case where q = 0. */
+ if (qq == 0.0)
+ {
+ norm = 1.0;
+ result_array[ii] = sin(order*zz);
+
+ continue;
+ }
+
+ /* 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 (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
+ {
+ result_array[ii] += coeff[jj]*sin(2.0*(jj + 1)*zz);
+ norm += coeff[jj]*coeff[jj];
+ }
+ }
+ else
+ {
+ for (jj=0; jj<GSL_SF_MATHIEU_COEFF; jj++)
+ {
+ result_array[ii] += coeff[jj]*sin((2.0*jj + 1.0)*zz);
+ norm += coeff[jj]*coeff[jj];
+ }
+ }
+
+ norm = sqrt(norm);
+ result_array[ii] /= norm;
+ }
+
+ return GSL_SUCCESS;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-01 07:57:33 +0000