|
Packit |
67cb25 |
/* specfunc/mathieu_radfunc.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2002 Lowell Johnson
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: L. Johnson */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_mathieu.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_mathieu_Mc_e(int kind, int order, double qq, double zz,
|
|
Packit |
67cb25 |
gsl_sf_result *result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int even_odd, kk, mm, status;
|
|
Packit |
67cb25 |
double maxerr = 1e-14, amax, pi = M_PI, fn, factor;
|
|
Packit |
67cb25 |
double coeff[GSL_SF_MATHIEU_COEFF], fc;
|
|
Packit |
67cb25 |
double j1c, z2c, j1pc, z2pc;
|
|
Packit |
67cb25 |
double u1, u2;
|
|
Packit |
67cb25 |
gsl_sf_result aa;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check for out of bounds parameters. */
|
|
Packit |
67cb25 |
if (qq <= 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("q must be greater than zero", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (kind < 1 || kind > 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
mm = 0;
|
|
Packit |
67cb25 |
amax = 0.0;
|
|
Packit |
67cb25 |
fn = 0.0;
|
|
Packit |
67cb25 |
u1 = sqrt(qq)*exp(-1.0*zz);
|
|
Packit |
67cb25 |
u2 = sqrt(qq)*exp(zz);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
even_odd = 0;
|
|
Packit |
67cb25 |
if (order % 2 != 0)
|
|
Packit |
67cb25 |
even_odd = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the characteristic value. */
|
|
Packit |
67cb25 |
status = gsl_sf_mathieu_a_e(order, qq, &aa);
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the series coefficients. */
|
|
Packit |
67cb25 |
status = gsl_sf_mathieu_a_coeff(order, qq, aa.val, coeff);
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (even_odd == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (kk=0; kk
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
amax = GSL_MAX(amax, fabs(coeff[kk]));
|
|
Packit |
67cb25 |
if (fabs(coeff[kk])/amax < maxerr)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
j1c = gsl_sf_bessel_Jn(kk, u1);
|
|
Packit |
67cb25 |
if (kind == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Jn(kk, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else /* kind = 2 */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Yn(kk, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fc = pow(-1.0, 0.5*order+kk)*coeff[kk];
|
|
Packit |
67cb25 |
fn += fc*j1c*z2c;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fn *= sqrt(pi/2.0)/coeff[0];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (kk=0; kk
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
amax = GSL_MAX(amax, fabs(coeff[kk]));
|
|
Packit |
67cb25 |
if (fabs(coeff[kk])/amax < maxerr)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
j1c = gsl_sf_bessel_Jn(kk, u1);
|
|
Packit |
67cb25 |
j1pc = gsl_sf_bessel_Jn(kk+1, u1);
|
|
Packit |
67cb25 |
if (kind == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Jn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Jn(kk+1, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else /* kind = 2 */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Yn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Yn(kk+1, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
|
|
Packit |
67cb25 |
fn += fc*(j1c*z2pc + j1pc*z2c);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fn *= sqrt(pi/2.0)/coeff[0];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = fn;
|
|
Packit |
67cb25 |
result->err = 2.0*GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
factor = fabs(fn);
|
|
Packit |
67cb25 |
if (factor > 1.0)
|
|
Packit |
67cb25 |
result->err *= factor;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_mathieu_Ms_e(int kind, int order, double qq, double zz,
|
|
Packit |
67cb25 |
gsl_sf_result *result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int even_odd, kk, mm, status;
|
|
Packit |
67cb25 |
double maxerr = 1e-14, amax, pi = M_PI, fn, factor;
|
|
Packit |
67cb25 |
double coeff[GSL_SF_MATHIEU_COEFF], fc;
|
|
Packit |
67cb25 |
double j1c, z2c, j1mc, z2mc, j1pc, z2pc;
|
|
Packit |
67cb25 |
double u1, u2;
|
|
Packit |
67cb25 |
gsl_sf_result aa;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check for out of bounds parameters. */
|
|
Packit |
67cb25 |
if (qq <= 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("q must be greater than zero", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (kind < 1 || kind > 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Handle the trivial cases where order = 0. */
|
|
Packit |
67cb25 |
if (order == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
mm = 0;
|
|
Packit |
67cb25 |
amax = 0.0;
|
|
Packit |
67cb25 |
fn = 0.0;
|
|
Packit |
67cb25 |
u1 = sqrt(qq)*exp(-1.0*zz);
|
|
Packit |
67cb25 |
u2 = sqrt(qq)*exp(zz);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
even_odd = 0;
|
|
Packit |
67cb25 |
if (order % 2 != 0)
|
|
Packit |
67cb25 |
even_odd = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the characteristic value. */
|
|
Packit |
67cb25 |
status = gsl_sf_mathieu_b_e(order, qq, &aa);
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the series coefficients. */
|
|
Packit |
67cb25 |
status = gsl_sf_mathieu_b_coeff(order, qq, aa.val, coeff);
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (even_odd == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (kk=0; kk
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
amax = GSL_MAX(amax, fabs(coeff[kk]));
|
|
Packit |
67cb25 |
if (fabs(coeff[kk])/amax < maxerr)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
j1mc = gsl_sf_bessel_Jn(kk, u1);
|
|
Packit |
67cb25 |
j1pc = gsl_sf_bessel_Jn(kk+2, u1);
|
|
Packit |
67cb25 |
if (kind == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2mc = gsl_sf_bessel_Jn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Jn(kk+2, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else /* kind = 2 */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2mc = gsl_sf_bessel_Yn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Yn(kk+2, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fc = pow(-1.0, 0.5*order+kk+1)*coeff[kk];
|
|
Packit |
67cb25 |
fn += fc*(j1mc*z2pc - j1pc*z2mc);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fn *= sqrt(pi/2.0)/coeff[0];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (kk=0; kk
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
amax = GSL_MAX(amax, fabs(coeff[kk]));
|
|
Packit |
67cb25 |
if (fabs(coeff[kk])/amax < maxerr)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
j1c = gsl_sf_bessel_Jn(kk, u1);
|
|
Packit |
67cb25 |
j1pc = gsl_sf_bessel_Jn(kk+1, u1);
|
|
Packit |
67cb25 |
if (kind == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Jn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Jn(kk+1, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else /* kind = 2 */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Yn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Yn(kk+1, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
|
|
Packit |
67cb25 |
fn += fc*(j1c*z2pc - j1pc*z2c);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fn *= sqrt(pi/2.0)/coeff[0];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = fn;
|
|
Packit |
67cb25 |
result->err = 2.0*GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
factor = fabs(fn);
|
|
Packit |
67cb25 |
if (factor > 1.0)
|
|
Packit |
67cb25 |
result->err *= factor;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_mathieu_Mc_array(int kind, int nmin, int nmax, double qq,
|
|
Packit |
67cb25 |
double zz, gsl_sf_mathieu_workspace *work,
|
|
Packit |
67cb25 |
double result_array[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int even_odd, order, ii, kk, mm, status;
|
|
Packit |
67cb25 |
double maxerr = 1e-14, amax, pi = M_PI, fn;
|
|
Packit |
67cb25 |
double coeff[GSL_SF_MATHIEU_COEFF], fc;
|
|
Packit |
67cb25 |
double j1c, z2c, j1pc, z2pc;
|
|
Packit |
67cb25 |
double u1, u2;
|
|
Packit |
67cb25 |
double *aa = work->aa;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialize the result array to zeroes. */
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
result_array[ii] = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check for out of bounds parameters. */
|
|
Packit |
67cb25 |
if (qq <= 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("q must be greater than zero", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (kind < 1 || kind > 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
mm = 0;
|
|
Packit |
67cb25 |
amax = 0.0;
|
|
Packit |
67cb25 |
u1 = sqrt(qq)*exp(-1.0*zz);
|
|
Packit |
67cb25 |
u2 = sqrt(qq)*exp(zz);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute all eigenvalues up to nmax. */
|
|
Packit |
67cb25 |
gsl_sf_mathieu_a_array(0, nmax, qq, work, aa);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii=0, order=nmin; order<=nmax; ii++, order++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
fn = 0.0;
|
|
Packit |
67cb25 |
even_odd = 0;
|
|
Packit |
67cb25 |
if (order % 2 != 0)
|
|
Packit |
67cb25 |
even_odd = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the series coefficients. */
|
|
Packit |
67cb25 |
status = gsl_sf_mathieu_a_coeff(order, qq, aa[order], coeff);
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (even_odd == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (kk=0; kk
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
amax = GSL_MAX(amax, fabs(coeff[kk]));
|
|
Packit |
67cb25 |
if (fabs(coeff[kk])/amax < maxerr)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
j1c = gsl_sf_bessel_Jn(kk, u1);
|
|
Packit |
67cb25 |
if (kind == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Jn(kk, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else /* kind = 2 */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Yn(kk, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fc = pow(-1.0, 0.5*order+kk)*coeff[kk];
|
|
Packit |
67cb25 |
fn += fc*j1c*z2c;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fn *= sqrt(pi/2.0)/coeff[0];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (kk=0; kk
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
amax = GSL_MAX(amax, fabs(coeff[kk]));
|
|
Packit |
67cb25 |
if (fabs(coeff[kk])/amax < maxerr)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
j1c = gsl_sf_bessel_Jn(kk, u1);
|
|
Packit |
67cb25 |
j1pc = gsl_sf_bessel_Jn(kk+1, u1);
|
|
Packit |
67cb25 |
if (kind == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Jn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Jn(kk+1, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else /* kind = 2 */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Yn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Yn(kk+1, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
|
|
Packit |
67cb25 |
fn += fc*(j1c*z2pc + j1pc*z2c);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fn *= sqrt(pi/2.0)/coeff[0];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result_array[ii] = fn;
|
|
Packit |
67cb25 |
} /* order loop */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_mathieu_Ms_array(int kind, int nmin, int nmax, double qq,
|
|
Packit |
67cb25 |
double zz, gsl_sf_mathieu_workspace *work,
|
|
Packit |
67cb25 |
double result_array[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int even_odd, order, ii, kk, mm, status;
|
|
Packit |
67cb25 |
double maxerr = 1e-14, amax, pi = M_PI, fn;
|
|
Packit |
67cb25 |
double coeff[GSL_SF_MATHIEU_COEFF], fc;
|
|
Packit |
67cb25 |
double j1c, z2c, j1mc, z2mc, j1pc, z2pc;
|
|
Packit |
67cb25 |
double u1, u2;
|
|
Packit |
67cb25 |
double *bb = work->bb;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialize the result array to zeroes. */
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
result_array[ii] = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check for out of bounds parameters. */
|
|
Packit |
67cb25 |
if (qq <= 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("q must be greater than zero", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (kind < 1 || kind > 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("kind must be 1 or 2", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
mm = 0;
|
|
Packit |
67cb25 |
amax = 0.0;
|
|
Packit |
67cb25 |
u1 = sqrt(qq)*exp(-1.0*zz);
|
|
Packit |
67cb25 |
u2 = sqrt(qq)*exp(zz);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute all eigenvalues up to nmax. */
|
|
Packit |
67cb25 |
gsl_sf_mathieu_b_array(0, nmax, qq, work, bb);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii=0, order=nmin; order<=nmax; ii++, order++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
fn = 0.0;
|
|
Packit |
67cb25 |
even_odd = 0;
|
|
Packit |
67cb25 |
if (order % 2 != 0)
|
|
Packit |
67cb25 |
even_odd = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Handle the trivial cases where order = 0. */
|
|
Packit |
67cb25 |
if (order == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result_array[ii] = 0.0;
|
|
Packit |
67cb25 |
continue;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the series coefficients. */
|
|
Packit |
67cb25 |
status = gsl_sf_mathieu_b_coeff(order, qq, bb[order], coeff);
|
|
Packit |
67cb25 |
if (status != GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (even_odd == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (kk=0; kk
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
amax = GSL_MAX(amax, fabs(coeff[kk]));
|
|
Packit |
67cb25 |
if (fabs(coeff[kk])/amax < maxerr)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
j1mc = gsl_sf_bessel_Jn(kk, u1);
|
|
Packit |
67cb25 |
j1pc = gsl_sf_bessel_Jn(kk+2, u1);
|
|
Packit |
67cb25 |
if (kind == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2mc = gsl_sf_bessel_Jn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Jn(kk+2, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else /* kind = 2 */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2mc = gsl_sf_bessel_Yn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Yn(kk+2, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fc = pow(-1.0, 0.5*order+kk+1)*coeff[kk];
|
|
Packit |
67cb25 |
fn += fc*(j1mc*z2pc - j1pc*z2mc);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fn *= sqrt(pi/2.0)/coeff[0];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (kk=0; kk
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
amax = GSL_MAX(amax, fabs(coeff[kk]));
|
|
Packit |
67cb25 |
if (fabs(coeff[kk])/amax < maxerr)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
j1c = gsl_sf_bessel_Jn(kk, u1);
|
|
Packit |
67cb25 |
j1pc = gsl_sf_bessel_Jn(kk+1, u1);
|
|
Packit |
67cb25 |
if (kind == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Jn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Jn(kk+1, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else /* kind = 2 */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
z2c = gsl_sf_bessel_Yn(kk, u2);
|
|
Packit |
67cb25 |
z2pc = gsl_sf_bessel_Yn(kk+1, u2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fc = pow(-1.0, 0.5*(order-1)+kk)*coeff[kk];
|
|
Packit |
67cb25 |
fn += fc*(j1c*z2pc - j1pc*z2c);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
fn *= sqrt(pi/2.0)/coeff[0];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result_array[ii] = fn;
|
|
Packit |
67cb25 |
} /* order loop */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "eval.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_mathieu_Mc(int kind, int order, double qq, double zz)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_mathieu_Mc_e(kind, order, qq, zz, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_mathieu_Ms(int kind, int order, double qq, double zz)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_mathieu_Ms_e(kind, order, qq, zz, &result));
|
|
Packit |
67cb25 |
}
|