/* 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 #include #include #include #include #include "mathieu.h" #define MAX(a, b) ((a) > (b) ? (a) : (b)) int gsl_sf_mathieu_mc_1(int order, double qq, double zz, gsl_sf_result *result) { int even_odd, kk, mm, status; double maxerr = 1e-14, amax, pi = acos(-1.0), fn; double coeff[NUM_MATHIEU_COEFF], aa, fc, fj, fjp; double j1c, j2c, j1pc, j2pc; double u1, u2; /* Check for out of bounds parameters. */ if (qq == 0.0) { return GSL_FAILURE; } mm = 0; amax = 0.0; fn = 0.0; /* rm3 = (0.0,0.0); */ /* rm1p = 0.0; */ /* rm3p = (0.0,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_c_charv(order, qq, &aa); if (status != GSL_SUCCESS) { return status; } /* Compute the series coefficients. */ status = gsl_sf_mathieu_c_coeff(order, qq, aa, coeff); if (status != GSL_SUCCESS) { return status; } if (even_odd == 0) { for (kk=0; kkval = fn; result->err = GSL_DBL_EPSILON*fabs(fn); return GSL_SUCCESS; } int gsl_sf_mathieu_ms_1(int order, double qq, double zz, gsl_sf_result *result) { int even_odd, kk, mm, status; double maxerr = 1e-14, amax, pi = acos(-1.0), fn; double coeff[NUM_MATHIEU_COEFF], aa, fc, fj, fjp, fjm; double j1c, j2c, j1mc, j2mc, j1pc, j2pc; double u1, u2; /* Check for out of bounds parameters. */ if (qq == 0.0) { return GSL_FAILURE; } mm = 0; amax = 0.0; fn = 0.0; /* rm3 = (0.0,0.0); */ /* rm1p = 0.0; */ /* rm3p = (0.0,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_s_charv(order, qq, &aa); if (status != GSL_SUCCESS) { return status; } /* Compute the series coefficients. */ status = gsl_sf_mathieu_s_coeff(order, qq, aa, coeff); if (status != GSL_SUCCESS) { return status; } if (even_odd == 0) { for (kk=0; kkval = fn; result->err = GSL_DBL_EPSILON*fabs(fn); return GSL_SUCCESS; }