123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349 |
- /* specfunc/mathieu_coeff.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 3 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 "gsl__config.h"
- #include <stdlib.h>
- #include <math.h>
- #include "gsl_sf_mathieu.h"
- /*****************************************************************************
- * backward_recurse
- *
- * Purpose:
- ****************************************************************************/
- static void backward_recurse_c(double aa, double qq, double xx, double *ff,
- double *gx, int even_odd, int ni)
- {
- int ii, nn;
- double g1;
- g1 = *gx;
- ff[ni] = xx;
- if (even_odd == 0)
- {
- for (ii=0; ii<ni; ii++)
- {
- nn = GSL_SF_MATHIEU_COEFF - ii - 1;
- ff[ni-ii-1] = -1.0/((4*nn*nn - aa)/qq + ff[ni-ii]);
- }
- if (ni == GSL_SF_MATHIEU_COEFF - 1)
- ff[0] *= 2.0;
- }
- else
- {
- for (ii=0; ii<ni; ii++)
- {
- nn = GSL_SF_MATHIEU_COEFF - ii - 1;
- ff[ni-ii-1] = -1.0/(((2*nn + 1)*(2*nn + 1) - aa)/qq + ff[ni-ii]);
- }
- }
- *gx = ff[0] - g1;
- }
- static void backward_recurse_s(double aa, double qq, double xx, double *ff,
- double *gx, int even_odd, int ni)
- {
- int ii, nn;
- double g1;
- g1 = *gx;
- ff[ni] = xx;
- if (even_odd == 0)
- {
- for (ii=0; ii<ni; ii++)
- {
- nn = GSL_SF_MATHIEU_COEFF - ii - 1;
- ff[ni-ii-1] = -1.0/((4*(nn + 1)*(nn + 1) - aa)/qq + ff[ni-ii]);
- }
- }
- else
- {
- for (ii=0; ii<ni; ii++)
- {
- nn = GSL_SF_MATHIEU_COEFF - ii - 1;
- ff[ni-ii-1] = -1.0/(((2*nn + 1)*(2*nn + 1) - aa)/qq + ff[ni-ii]);
- }
- }
- *gx = ff[0] - g1;
- }
- int gsl_sf_mathieu_a_coeff(int order, double qq, double aa, double coeff[])
- {
- int ni, nn, ii, even_odd;
- double eps, g1, g2, x1, x2, e1, e2, de, xh, sum, ratio,
- ff[GSL_SF_MATHIEU_COEFF];
- eps = 1e-14;
- coeff[0] = 1.0;
-
- even_odd = 0;
- if (order % 2 != 0)
- even_odd = 1;
- /* If the coefficient array is not large enough to hold all necessary
- coefficients, error out. */
- if (order > GSL_SF_MATHIEU_COEFF)
- return GSL_FAILURE;
-
- /* Handle the trivial case where q = 0. */
- if (qq == 0.0)
- {
- for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
- coeff[ii] = 0.0;
- coeff[order/2] = 1.0;
-
- return GSL_SUCCESS;
- }
-
- if (order < 5)
- {
- nn = 0;
- sum = 0.0;
- if (even_odd == 0)
- ratio = aa/qq;
- else
- ratio = (aa - 1 - qq)/qq;
- }
- else
- {
- if (even_odd == 0)
- {
- coeff[1] = aa/qq;
- coeff[2] = (aa - 4)/qq*coeff[1] - 2;
- sum = coeff[0] + coeff[1] + coeff[2];
- for (ii=3; ii<order/2+1; ii++)
- {
- coeff[ii] = (aa - 4*(ii - 1)*(ii - 1))/qq*coeff[ii-1] -
- coeff[ii-2];
- sum += coeff[ii];
- }
- }
- else
- {
- coeff[1] = (aa - 1)/qq - 1;
- sum = coeff[0] + coeff[1];
- for (ii=2; ii<order/2+1; ii++)
- {
- coeff[ii] = (aa - (2*ii - 1)*(2*ii - 1))/qq*coeff[ii-1] -
- coeff[ii-2];
- sum += coeff[ii];
- }
- }
- nn = ii - 1;
- ratio = coeff[nn]/coeff[nn-1];
- }
-
- ni = GSL_SF_MATHIEU_COEFF - nn - 1;
- /* Compute first two points to start root-finding. */
- if (even_odd == 0)
- x1 = -qq/(4.0*GSL_SF_MATHIEU_COEFF*GSL_SF_MATHIEU_COEFF);
- else
- x1 = -qq/((2.0*GSL_SF_MATHIEU_COEFF + 1.0)*(2.0*GSL_SF_MATHIEU_COEFF + 1.0));
- g1 = ratio;
- backward_recurse_c(aa, qq, x1, ff, &g1, even_odd, ni);
- x2 = g1;
- g2 = ratio;
- backward_recurse_c(aa, qq, x2, ff, &g2, even_odd, ni);
- /* Find the root. */
- while (1)
- {
- /* Compute the relative error. */
- e1 = g1 - x1;
- e2 = g2 - x2;
- de = e1 - e2;
- /* If we are close enough to the root, break... */
- if (fabs(de) < eps)
- break;
- /* Otherwise, determine the next guess and try again. */
- xh = (e1*x2 - e2*x1)/de;
- x1 = x2;
- g1 = g2;
- x2 = xh;
- g2 = ratio;
- backward_recurse_c(aa, qq, x2, ff, &g2, even_odd, ni);
- }
- /* Compute the rest of the coefficients. */
- sum += coeff[nn];
- for (ii=nn+1; ii<GSL_SF_MATHIEU_COEFF; ii++)
- {
- coeff[ii] = ff[ii-nn-1]*coeff[ii-1];
- sum += coeff[ii];
- /* If the coefficients are getting really small, set the remainder
- to zero. */
- if (fabs(coeff[ii]) < 1e-20)
- {
- for (; ii<GSL_SF_MATHIEU_COEFF;)
- coeff[ii++] = 0.0;
- }
- }
-
- /* Normalize the coefficients. */
- for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
- coeff[ii] /= sum;
- return GSL_SUCCESS;
- }
- int gsl_sf_mathieu_b_coeff(int order, double qq, double aa, double coeff[])
- {
- int ni, nn, ii, even_odd;
- double eps, g1, g2, x1, x2, e1, e2, de, xh, sum, ratio,
- ff[GSL_SF_MATHIEU_COEFF];
- eps = 1e-10;
- coeff[0] = 1.0;
-
- even_odd = 0;
- if (order % 2 != 0)
- even_odd = 1;
- /* If the coefficient array is not large enough to hold all necessary
- coefficients, error out. */
- if (order > GSL_SF_MATHIEU_COEFF)
- return GSL_FAILURE;
-
- /* Handle the trivial case where q = 0. */
- if (qq == 0.0)
- {
- for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
- coeff[ii] = 0.0;
- coeff[(order-1)/2] = 1.0;
-
- return GSL_SUCCESS;
- }
-
- if (order < 5)
- {
- nn = 0;
- sum = 0.0;
- if (even_odd == 0)
- ratio = (aa - 4)/qq;
- else
- ratio = (aa - 1 - qq)/qq;
- }
- else
- {
- if (even_odd == 0)
- {
- coeff[1] = (aa - 4)/qq;
- sum = 2*coeff[0] + 4*coeff[1];
- for (ii=2; ii<order/2; ii++)
- {
- coeff[ii] = (aa - 4*ii*ii)/qq*coeff[ii-1] - coeff[ii-2];
- sum += 2*(ii + 1)*coeff[ii];
- }
- }
- else
- {
- coeff[1] = (aa - 1)/qq + 1;
- sum = coeff[0] + 3*coeff[1];
- for (ii=2; ii<order/2+1; ii++)
- {
- coeff[ii] = (aa - (2*ii - 1)*(2*ii - 1))/qq*coeff[ii-1] -
- coeff[ii-2];
- sum += (2*(ii + 1) - 1)*coeff[ii];
- }
- }
- nn = ii - 1;
- ratio = coeff[nn]/coeff[nn-1];
- }
-
- ni = GSL_SF_MATHIEU_COEFF - nn - 1;
- /* Compute first two points to start root-finding. */
- if (even_odd == 0)
- x1 = -qq/(4.0*(GSL_SF_MATHIEU_COEFF + 1.0)*(GSL_SF_MATHIEU_COEFF + 1.0));
- else
- x1 = -qq/((2.0*GSL_SF_MATHIEU_COEFF + 1.0)*(2.0*GSL_SF_MATHIEU_COEFF + 1.0));
- g1 = ratio;
- backward_recurse_s(aa, qq, x1, ff, &g1, even_odd, ni);
- x2 = g1;
- g2 = ratio;
- backward_recurse_s(aa, qq, x2, ff, &g2, even_odd, ni);
- /* Find the root. */
- while (1)
- {
- /* Compute the relative error. */
- e1 = g1 - x1;
- e2 = g2 - x2;
- de = e1 - e2;
- /* If we are close enough to the root, break... */
- if (fabs(de) < eps)
- break;
- /* Otherwise, determine the next guess and try again. */
- xh = (e1*x2 - e2*x1)/de;
- x1 = x2;
- g1 = g2;
- x2 = xh;
- g2 = ratio;
- backward_recurse_s(aa, qq, x2, ff, &g2, even_odd, ni);
- }
- /* Compute the rest of the coefficients. */
- sum += 2*(nn + 1)*coeff[nn];
- for (ii=nn+1; ii<GSL_SF_MATHIEU_COEFF; ii++)
- {
- coeff[ii] = ff[ii-nn-1]*coeff[ii-1];
- sum += 2*(ii + 1)*coeff[ii];
- /* If the coefficients are getting really small, set the remainder
- to zero. */
- if (fabs(coeff[ii]) < 1e-20)
- {
- for (; ii<GSL_SF_MATHIEU_COEFF;)
- coeff[ii++] = 0.0;
- }
- }
-
- /* Normalize the coefficients. */
- for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
- coeff[ii] /= sum;
- return GSL_SUCCESS;
- }
|