|
Packit |
67cb25 |
/* specfunc/mathieu_charv.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2002, 2009 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 <stdio.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_eigen.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_mathieu.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* prototypes */
|
|
Packit |
67cb25 |
static double solve_cubic(double c2, double c1, double c0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double ceer(int order, double qq, double aa, int nterms)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double term, term1;
|
|
Packit |
67cb25 |
int ii, n1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (order == 0)
|
|
Packit |
67cb25 |
term = 0.0;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
term = 2.0*qq*qq/aa;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (order != 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
n1 = order/2 - 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
term = qq*qq/(aa - 4.0*(ii+1)*(ii+1) - term);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
term += order*order;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
term1 = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
term1 = qq*qq/
|
|
Packit |
67cb25 |
(aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (order == 0)
|
|
Packit |
67cb25 |
term1 *= 2.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (term + term1 - aa);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double ceor(int order, double qq, double aa, int nterms)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double term, term1;
|
|
Packit |
67cb25 |
int ii, n1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
term = qq;
|
|
Packit |
67cb25 |
n1 = (int)((float)order/2.0 - 0.5);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
term = qq*qq/(aa - (2.0*ii + 1.0)*(2.0*ii + 1.0) - term);
|
|
Packit |
67cb25 |
term += order*order;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
term1 = 0.0;
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
term1 = qq*qq/
|
|
Packit |
67cb25 |
(aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (term + term1 - aa);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double seer(int order, double qq, double aa, int nterms)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double term, term1;
|
|
Packit |
67cb25 |
int ii, n1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
term = 0.0;
|
|
Packit |
67cb25 |
n1 = order/2 - 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
term = qq*qq/(aa - 4.0*(ii + 1)*(ii + 1) - term);
|
|
Packit |
67cb25 |
term += order*order;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
term1 = 0.0;
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
term1 = qq*qq/
|
|
Packit |
67cb25 |
(aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (term + term1 - aa);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double seor(int order, double qq, double aa, int nterms)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double term, term1;
|
|
Packit |
67cb25 |
int ii, n1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
term = -1.0*qq;
|
|
Packit |
67cb25 |
n1 = (int)((float)order/2.0 - 0.5);
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
term = qq*qq/(aa - (2.0*ii + 1.0)*(2.0*ii + 1.0) - term);
|
|
Packit |
67cb25 |
term += order*order;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
term1 = 0.0;
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
term1 = qq*qq/
|
|
Packit |
67cb25 |
(aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (term + term1 - aa);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*----------------------------------------------------------------------------
|
|
Packit |
67cb25 |
* Asymptotic and approximation routines for the characteristic value.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Adapted from F.A. Alhargan's paper,
|
|
Packit |
67cb25 |
* "Algorithms for the Computation of All Mathieu Functions of Integer
|
|
Packit |
67cb25 |
* Orders," ACM Transactions on Mathematical Software, Vol. 26, No. 3,
|
|
Packit |
67cb25 |
* September 2000, pp. 390-407.
|
|
Packit |
67cb25 |
*--------------------------------------------------------------------------*/
|
|
Packit |
67cb25 |
static double asymptotic(int order, double qq)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double asymp;
|
|
Packit |
67cb25 |
double nn, n2, n4, n6;
|
|
Packit |
67cb25 |
double hh, ah, ah2, ah3, ah4, ah5;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Set up temporary variables to simplify the readability. */
|
|
Packit |
67cb25 |
nn = 2*order + 1;
|
|
Packit |
67cb25 |
n2 = nn*nn;
|
|
Packit |
67cb25 |
n4 = n2*n2;
|
|
Packit |
67cb25 |
n6 = n4*n2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
hh = 2*sqrt(qq);
|
|
Packit |
67cb25 |
ah = 16*hh;
|
|
Packit |
67cb25 |
ah2 = ah*ah;
|
|
Packit |
67cb25 |
ah3 = ah2*ah;
|
|
Packit |
67cb25 |
ah4 = ah3*ah;
|
|
Packit |
67cb25 |
ah5 = ah4*ah;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Equation 38, p. 397 of Alhargan's paper. */
|
|
Packit |
67cb25 |
asymp = -2*qq + nn*hh - 0.125*(n2 + 1);
|
|
Packit |
67cb25 |
asymp -= 0.25*nn*( n2 + 3)/ah;
|
|
Packit |
67cb25 |
asymp -= 0.25* ( 5*n4 + 34*n2 + 9)/ah2;
|
|
Packit |
67cb25 |
asymp -= 0.25*nn*( 33*n4 + 410*n2 + 405)/ah3;
|
|
Packit |
67cb25 |
asymp -= ( 63*n6 + 1260*n4 + 2943*n2 + 486)/ah4;
|
|
Packit |
67cb25 |
asymp -= nn*(527*n6 + 15617*n4 + 69001*n2 + 41607)/ah5;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return asymp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
|
|
Packit |
67cb25 |
static double solve_cubic(double c2, double c1, double c0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double qq, rr, ww, ss, tt;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
qq = (3*c1 - c2*c2)/9;
|
|
Packit |
67cb25 |
rr = (9*c2*c1 - 27*c0 - 2*c2*c2*c2)/54;
|
|
Packit |
67cb25 |
ww = qq*qq*qq + rr*rr;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (ww >= 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double t1 = rr + sqrt(ww);
|
|
Packit |
67cb25 |
ss = fabs(t1)/t1*pow(fabs(t1), 1/3.);
|
|
Packit |
67cb25 |
t1 = rr - sqrt(ww);
|
|
Packit |
67cb25 |
tt = fabs(t1)/t1*pow(fabs(t1), 1/3.);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double theta = acos(rr/sqrt(-qq*qq*qq));
|
|
Packit |
67cb25 |
ss = 2*sqrt(-qq)*cos((theta + 4*M_PI)/3.);
|
|
Packit |
67cb25 |
tt = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (ss + tt - c2/3);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute an initial approximation for the characteristic value. */
|
|
Packit |
67cb25 |
static double approx_c(int order, double qq)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double approx;
|
|
Packit |
67cb25 |
double c0, c1, c2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (order < 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_VAL("Undefined order for Mathieu function", GSL_EINVAL, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
switch (order)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
case 0:
|
|
Packit |
67cb25 |
if (qq <= 4)
|
|
Packit |
67cb25 |
return (2 - sqrt(4 + 2*qq*qq)); /* Eqn. 31 */
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return asymptotic(order, qq);
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
case 1:
|
|
Packit |
67cb25 |
if (qq <= 4)
|
|
Packit |
67cb25 |
return (5 + 0.5*(qq - sqrt(5*qq*qq - 16*qq + 64))); /* Eqn. 32 */
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return asymptotic(order, qq);
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
case 2:
|
|
Packit |
67cb25 |
if (qq <= 3)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
c2 = -8.0; /* Eqn. 33 */
|
|
Packit |
67cb25 |
c1 = -48 - 3*qq*qq;
|
|
Packit |
67cb25 |
c0 = 20*qq*qq;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return asymptotic(order, qq);
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
case 3:
|
|
Packit |
67cb25 |
if (qq <= 6.25)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
c2 = -qq - 8; /* Eqn. 34 */
|
|
Packit |
67cb25 |
c1 = 16*qq - 128 - 2*qq*qq;
|
|
Packit |
67cb25 |
c0 = qq*qq*(qq + 8);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return asymptotic(order, qq);
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
default:
|
|
Packit |
67cb25 |
if (order < 70)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (1.7*order > 2*sqrt(qq))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Eqn. 30 */
|
|
Packit |
67cb25 |
double n2 = (double)(order*order);
|
|
Packit |
67cb25 |
double n22 = (double)((n2 - 1)*(n2 - 1));
|
|
Packit |
67cb25 |
double q2 = qq*qq;
|
|
Packit |
67cb25 |
double q4 = q2*q2;
|
|
Packit |
67cb25 |
approx = n2 + 0.5*q2/(n2 - 1);
|
|
Packit |
67cb25 |
approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4));
|
|
Packit |
67cb25 |
approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/
|
|
Packit |
67cb25 |
(64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9));
|
|
Packit |
67cb25 |
if (1.4*order < 2*sqrt(qq))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
approx += asymptotic(order, qq);
|
|
Packit |
67cb25 |
approx *= 0.5;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
approx = asymptotic(order, qq);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return approx;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return order*order;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
|
|
Packit |
67cb25 |
approx = solve_cubic(c2, c1, c0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if ( approx < 0 && sqrt(qq) > 0.1*order )
|
|
Packit |
67cb25 |
return asymptotic(order-1, qq);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return (order*order + fabs(approx));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double approx_s(int order, double qq)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double approx;
|
|
Packit |
67cb25 |
double c0, c1, c2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (order < 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR_VAL("Undefined order for Mathieu function", GSL_EINVAL, 0.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
switch (order)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
case 1:
|
|
Packit |
67cb25 |
if (qq <= 4)
|
|
Packit |
67cb25 |
return (5 - 0.5*(qq + sqrt(5*qq*qq + 16*qq + 64))); /* Eqn. 35 */
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return asymptotic(order-1, qq);
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
case 2:
|
|
Packit |
67cb25 |
if (qq <= 5)
|
|
Packit |
67cb25 |
return (10 - sqrt(36 + qq*qq)); /* Eqn. 36 */
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return asymptotic(order-1, qq);
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
case 3:
|
|
Packit |
67cb25 |
if (qq <= 6.25)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
c2 = qq - 8; /* Eqn. 37 */
|
|
Packit |
67cb25 |
c1 = -128 - 16*qq - 2*qq*qq;
|
|
Packit |
67cb25 |
c0 = qq*qq*(8 - qq);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return asymptotic(order-1, qq);
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
default:
|
|
Packit |
67cb25 |
if (order < 70)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (1.7*order > 2*sqrt(qq))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Eqn. 30 */
|
|
Packit |
67cb25 |
double n2 = (double)(order*order);
|
|
Packit |
67cb25 |
double n22 = (double)((n2 - 1)*(n2 - 1));
|
|
Packit |
67cb25 |
double q2 = qq*qq;
|
|
Packit |
67cb25 |
double q4 = q2*q2;
|
|
Packit |
67cb25 |
approx = n2 + 0.5*q2/(n2 - 1);
|
|
Packit |
67cb25 |
approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4));
|
|
Packit |
67cb25 |
approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/
|
|
Packit |
67cb25 |
(64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9));
|
|
Packit |
67cb25 |
if (1.4*order < 2*sqrt(qq))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
approx += asymptotic(order-1, qq);
|
|
Packit |
67cb25 |
approx *= 0.5;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
approx = asymptotic(order-1, qq);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return approx;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return order*order;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
|
|
Packit |
67cb25 |
approx = solve_cubic(c2, c1, c0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if ( approx < 0 && sqrt(qq) > 0.1*order )
|
|
Packit |
67cb25 |
return asymptotic(order-1, qq);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return (order*order + fabs(approx));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_mathieu_a_e(int order, double qq, gsl_sf_result *result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int even_odd, nterms = 50, ii, counter = 0, maxcount = 1000;
|
|
Packit |
67cb25 |
int dir = 0; /* step direction for new search */
|
|
Packit |
67cb25 |
double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;
|
|
Packit |
67cb25 |
double aa_approx; /* current approximation for solution */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
even_odd = 0;
|
|
Packit |
67cb25 |
if (order % 2 != 0)
|
|
Packit |
67cb25 |
even_odd = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* If the argument is 0, then the coefficient is simply the square of
|
|
Packit |
67cb25 |
the order. */
|
|
Packit |
67cb25 |
if (qq == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->val = order*order;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Use symmetry characteristics of the functions to handle cases with
|
|
Packit |
67cb25 |
negative order and/or argument q. See Abramowitz & Stegun, 20.8.3. */
|
|
Packit |
67cb25 |
if (order < 0)
|
|
Packit |
67cb25 |
order *= -1;
|
|
Packit |
67cb25 |
if (qq < 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (even_odd == 0)
|
|
Packit |
67cb25 |
return gsl_sf_mathieu_a_e(order, -qq, result);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return gsl_sf_mathieu_b_e(order, -qq, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute an initial approximation for the characteristic value. */
|
|
Packit |
67cb25 |
aa_approx = approx_c(order, qq);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Save the original approximation for later comparison. */
|
|
Packit |
67cb25 |
aa_orig = aa = aa_approx;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Loop as long as the final value is not near the approximate value
|
|
Packit |
67cb25 |
(with a max limit to avoid potential infinite loop). */
|
|
Packit |
67cb25 |
while (counter < maxcount)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
a1 = aa + 0.001;
|
|
Packit |
67cb25 |
ii = 0;
|
|
Packit |
67cb25 |
if (even_odd == 0)
|
|
Packit |
67cb25 |
fa1 = ceer(order, qq, a1, nterms);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
fa1 = ceor(order, qq, a1, nterms);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (;;)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (even_odd == 0)
|
|
Packit |
67cb25 |
fa = ceer(order, qq, aa, nterms);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
fa = ceor(order, qq, aa, nterms);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
a2 = a1;
|
|
Packit |
67cb25 |
a1 = aa;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (fa == fa1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
aa -= (aa - a2)/(fa - fa1)*fa;
|
|
Packit |
67cb25 |
dela = fabs(aa - a2);
|
|
Packit |
67cb25 |
if (dela < GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (ii > 40)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->err = dela;
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
fa1 = fa;
|
|
Packit |
67cb25 |
ii++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* If the solution found is not near the original approximation,
|
|
Packit |
67cb25 |
tweak the approximate value, and try again. */
|
|
Packit |
67cb25 |
if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)) ||
|
|
Packit |
67cb25 |
(order > 10 && fabs(aa - aa_orig) > 1.5*order))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
counter++;
|
|
Packit |
67cb25 |
if (counter == maxcount)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->err = fabs(aa - aa_orig);
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (aa > aa_orig)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (dir == 1)
|
|
Packit |
67cb25 |
da /= 2;
|
|
Packit |
67cb25 |
dir = -1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (dir == -1)
|
|
Packit |
67cb25 |
da /= 2;
|
|
Packit |
67cb25 |
dir = 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
aa_approx += dir*da*counter;
|
|
Packit |
67cb25 |
aa = aa_approx;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
continue;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = aa;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* If we went through the maximum number of retries and still didn't
|
|
Packit |
67cb25 |
find the solution, let us know. */
|
|
Packit |
67cb25 |
if (counter == maxcount)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("Wrong characteristic Mathieu value", GSL_EFAILED);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_mathieu_b_e(int order, double qq, gsl_sf_result *result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int even_odd, nterms = 50, ii, counter = 0, maxcount = 1000;
|
|
Packit |
67cb25 |
int dir = 0; /* step direction for new search */
|
|
Packit |
67cb25 |
double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;
|
|
Packit |
67cb25 |
double aa_approx; /* current approximation for solution */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
even_odd = 0;
|
|
Packit |
67cb25 |
if (order % 2 != 0)
|
|
Packit |
67cb25 |
even_odd = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The order cannot be 0. */
|
|
Packit |
67cb25 |
if (order == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("Characteristic value undefined for order 0", GSL_EFAILED);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* If the argument is 0, then the coefficient is simply the square of
|
|
Packit |
67cb25 |
the order. */
|
|
Packit |
67cb25 |
if (qq == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->val = order*order;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Use symmetry characteristics of the functions to handle cases with
|
|
Packit |
67cb25 |
negative order and/or argument q. See Abramowitz & Stegun, 20.8.3. */
|
|
Packit |
67cb25 |
if (order < 0)
|
|
Packit |
67cb25 |
order *= -1;
|
|
Packit |
67cb25 |
if (qq < 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (even_odd == 0)
|
|
Packit |
67cb25 |
return gsl_sf_mathieu_b_e(order, -qq, result);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return gsl_sf_mathieu_a_e(order, -qq, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute an initial approximation for the characteristic value. */
|
|
Packit |
67cb25 |
aa_approx = approx_s(order, qq);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Save the original approximation for later comparison. */
|
|
Packit |
67cb25 |
aa_orig = aa = aa_approx;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Loop as long as the final value is not near the approximate value
|
|
Packit |
67cb25 |
(with a max limit to avoid potential infinite loop). */
|
|
Packit |
67cb25 |
while (counter < maxcount)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
a1 = aa + 0.001;
|
|
Packit |
67cb25 |
ii = 0;
|
|
Packit |
67cb25 |
if (even_odd == 0)
|
|
Packit |
67cb25 |
fa1 = seer(order, qq, a1, nterms);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
fa1 = seor(order, qq, a1, nterms);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (;;)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (even_odd == 0)
|
|
Packit |
67cb25 |
fa = seer(order, qq, aa, nterms);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
fa = seor(order, qq, aa, nterms);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
a2 = a1;
|
|
Packit |
67cb25 |
a1 = aa;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (fa == fa1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
aa -= (aa - a2)/(fa - fa1)*fa;
|
|
Packit |
67cb25 |
dela = fabs(aa - a2);
|
|
Packit |
67cb25 |
if (dela < 1e-18)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (ii > 40)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->err = dela;
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
fa1 = fa;
|
|
Packit |
67cb25 |
ii++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* If the solution found is not near the original approximation,
|
|
Packit |
67cb25 |
tweak the approximate value, and try again. */
|
|
Packit |
67cb25 |
if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)) ||
|
|
Packit |
67cb25 |
(order > 10 && fabs(aa - aa_orig) > 1.5*order))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
counter++;
|
|
Packit |
67cb25 |
if (counter == maxcount)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->err = fabs(aa - aa_orig);
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (aa > aa_orig)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (dir == 1)
|
|
Packit |
67cb25 |
da /= 2;
|
|
Packit |
67cb25 |
dir = -1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (dir == -1)
|
|
Packit |
67cb25 |
da /= 2;
|
|
Packit |
67cb25 |
dir = 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
aa_approx += dir*da*counter;
|
|
Packit |
67cb25 |
aa = aa_approx;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
continue;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = aa;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* If we went through the maximum number of retries and still didn't
|
|
Packit |
67cb25 |
find the solution, let us know. */
|
|
Packit |
67cb25 |
if (counter == maxcount)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("Wrong characteristic Mathieu value", GSL_EFAILED);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Eigenvalue solutions for characteristic values below. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* figi.c converted from EISPACK Fortran FIGI.F.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* given a nonsymmetric tridiagonal matrix such that the products
|
|
Packit |
67cb25 |
* of corresponding pairs of off-diagonal elements are all
|
|
Packit |
67cb25 |
* non-negative, this subroutine reduces it to a symmetric
|
|
Packit |
67cb25 |
* tridiagonal matrix with the same eigenvalues. if, further,
|
|
Packit |
67cb25 |
* a zero product only occurs when both factors are zero,
|
|
Packit |
67cb25 |
* the reduced matrix is similar to the original matrix.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* on input
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* n is the order of the matrix.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* t contains the input matrix. its subdiagonal is
|
|
Packit |
67cb25 |
* stored in the last n-1 positions of the first column,
|
|
Packit |
67cb25 |
* its diagonal in the n positions of the second column,
|
|
Packit |
67cb25 |
* and its superdiagonal in the first n-1 positions of
|
|
Packit |
67cb25 |
* the third column. t(1,1) and t(n,3) are arbitrary.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* on output
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* t is unaltered.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* d contains the diagonal elements of the symmetric matrix.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* e contains the subdiagonal elements of the symmetric
|
|
Packit |
67cb25 |
* matrix in its last n-1 positions. e(1) is not set.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* e2 contains the squares of the corresponding elements of e.
|
|
Packit |
67cb25 |
* e2 may coincide with e if the squares are not needed.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* ierr is set to
|
|
Packit |
67cb25 |
* zero for normal return,
|
|
Packit |
67cb25 |
* n+i if t(i,1)*t(i-1,3) is negative,
|
|
Packit |
67cb25 |
* -(3*n+i) if t(i,1)*t(i-1,3) is zero with one factor
|
|
Packit |
67cb25 |
* non-zero. in this case, the eigenvectors of
|
|
Packit |
67cb25 |
* the symmetric matrix are not simply related
|
|
Packit |
67cb25 |
* to those of t and should not be sought.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* questions and comments should be directed to burton s. garbow,
|
|
Packit |
67cb25 |
* mathematics and computer science div, argonne national laboratory
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* this version dated august 1983.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static int figi(int nn, double *tt, double *dd, double *ee,
|
|
Packit |
67cb25 |
double *e2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int ii;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (ii != 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
e2[ii] = tt[3*ii]*tt[3*(ii-1)+2];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (e2[ii] < 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* set error -- product of some pair of off-diagonal
|
|
Packit |
67cb25 |
elements is negative */
|
|
Packit |
67cb25 |
return (nn + ii);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (e2[ii] == 0.0 && (tt[3*ii] != 0.0 || tt[3*(ii-1)+2] != 0.0))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* set error -- product of some pair of off-diagonal
|
|
Packit |
67cb25 |
elements is zero with one member non-zero */
|
|
Packit |
67cb25 |
return (-1*(3*nn + ii));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ee[ii] = sqrt(e2[ii]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
dd[ii] = tt[3*ii+1];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_mathieu_a_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
unsigned int even_order = work->even_order, odd_order = work->odd_order,
|
|
Packit |
67cb25 |
extra_values = work->extra_values, ii, jj;
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
double *tt = work->tt, *dd = work->dd, *ee = work->ee, *e2 = work->e2,
|
|
Packit |
67cb25 |
*zz = work->zz, *aa = work->aa;
|
|
Packit |
67cb25 |
gsl_matrix_view mat, evec;
|
|
Packit |
67cb25 |
gsl_vector_view eval;
|
|
Packit |
67cb25 |
gsl_eigen_symmv_workspace *wmat = work->wmat;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (order_max > work->size || order_max <= order_min || order_min < 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("invalid range [order_min,order_max]", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Convert the nonsymmetric tridiagonal matrix to a symmetric tridiagonal
|
|
Packit |
67cb25 |
form. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tt[0] = 0.0;
|
|
Packit |
67cb25 |
tt[1] = 0.0;
|
|
Packit |
67cb25 |
tt[2] = qq;
|
|
Packit |
67cb25 |
for (ii=1; ii
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
tt[3*ii] = qq;
|
|
Packit |
67cb25 |
tt[3*ii+1] = 4*ii*ii;
|
|
Packit |
67cb25 |
tt[3*ii+2] = qq;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
tt[3*even_order-3] = qq;
|
|
Packit |
67cb25 |
tt[3*even_order-2] = 4*(even_order - 1)*(even_order - 1);
|
|
Packit |
67cb25 |
tt[3*even_order-1] = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tt[3] *= 2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = figi((signed int)even_order, tt, dd, ee, e2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("Internal error in tridiagonal Mathieu matrix", GSL_EFAILED);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Fill the period \pi matrix. */
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
zz[ii] = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
zz[0] = dd[0];
|
|
Packit |
67cb25 |
zz[1] = ee[1];
|
|
Packit |
67cb25 |
for (ii=1; ii
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
zz[ii*even_order+ii-1] = ee[ii];
|
|
Packit |
67cb25 |
zz[ii*even_order+ii] = dd[ii];
|
|
Packit |
67cb25 |
zz[ii*even_order+ii+1] = ee[ii+1];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
zz[even_order*(even_order-1)+even_order-2] = ee[even_order-1];
|
|
Packit |
67cb25 |
zz[even_order*even_order-1] = dd[even_order-1];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute (and sort) the eigenvalues of the matrix. */
|
|
Packit |
67cb25 |
mat = gsl_matrix_view_array(zz, even_order, even_order);
|
|
Packit |
67cb25 |
eval = gsl_vector_subvector(work->eval, 0, even_order);
|
|
Packit |
67cb25 |
evec = gsl_matrix_submatrix(work->evec, 0, 0, even_order, even_order);
|
|
Packit |
67cb25 |
gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
|
|
Packit |
67cb25 |
gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
aa[2*ii] = gsl_vector_get(&eval.vector, ii);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Fill the period 2\pi matrix. */
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
zz[ii] = 0.0;
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
for (jj=0; jj
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (ii == jj)
|
|
Packit |
67cb25 |
zz[ii*odd_order+jj] = (2*ii + 1)*(2*ii + 1);
|
|
Packit |
67cb25 |
else if (ii == jj + 1 || ii + 1 == jj)
|
|
Packit |
67cb25 |
zz[ii*odd_order+jj] = qq;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
zz[0] += qq;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute (and sort) the eigenvalues of the matrix. */
|
|
Packit |
67cb25 |
mat = gsl_matrix_view_array(zz, odd_order, odd_order);
|
|
Packit |
67cb25 |
eval = gsl_vector_subvector(work->eval, 0, odd_order);
|
|
Packit |
67cb25 |
evec = gsl_matrix_submatrix(work->evec, 0, 0, odd_order, odd_order);
|
|
Packit |
67cb25 |
gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
|
|
Packit |
67cb25 |
gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
aa[2*ii+1] = gsl_vector_get(&eval.vector, ii);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii = order_min ; ii <= order_max ; ii++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result_array[ii - order_min] = aa[ii];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_mathieu_b_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
unsigned int even_order = work->even_order-1, odd_order = work->odd_order,
|
|
Packit |
67cb25 |
extra_values = work->extra_values, ii, jj;
|
|
Packit |
67cb25 |
double *zz = work->zz, *bb = work->bb;
|
|
Packit |
67cb25 |
gsl_matrix_view mat, evec;
|
|
Packit |
67cb25 |
gsl_vector_view eval;
|
|
Packit |
67cb25 |
gsl_eigen_symmv_workspace *wmat = work->wmat;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (order_max > work->size || order_max <= order_min || order_min < 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("invalid range [order_min,order_max]", GSL_EINVAL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Fill the period \pi matrix. */
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
zz[ii] = 0.0;
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
for (jj=0; jj
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (ii == jj)
|
|
Packit |
67cb25 |
zz[ii*even_order+jj] = 4*(ii + 1)*(ii + 1);
|
|
Packit |
67cb25 |
else if (ii == jj + 1 || ii + 1 == jj)
|
|
Packit |
67cb25 |
zz[ii*even_order+jj] = qq;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute (and sort) the eigenvalues of the matrix. */
|
|
Packit |
67cb25 |
mat = gsl_matrix_view_array(zz, even_order, even_order);
|
|
Packit |
67cb25 |
eval = gsl_vector_subvector(work->eval, 0, even_order);
|
|
Packit |
67cb25 |
evec = gsl_matrix_submatrix(work->evec, 0, 0, even_order, even_order);
|
|
Packit |
67cb25 |
gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
|
|
Packit |
67cb25 |
gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
bb[0] = 0.0;
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
bb[2*(ii+1)] = gsl_vector_get(&eval.vector, ii);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Fill the period 2\pi matrix. */
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
zz[ii] = 0.0;
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
for (jj=0; jj
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (ii == jj)
|
|
Packit |
67cb25 |
zz[ii*odd_order+jj] = (2*ii + 1)*(2*ii + 1);
|
|
Packit |
67cb25 |
else if (ii == jj + 1 || ii + 1 == jj)
|
|
Packit |
67cb25 |
zz[ii*odd_order+jj] = qq;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
zz[0] -= qq;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute (and sort) the eigenvalues of the matrix. */
|
|
Packit |
67cb25 |
mat = gsl_matrix_view_array(zz, odd_order, odd_order);
|
|
Packit |
67cb25 |
eval = gsl_vector_subvector(work->eval, 0, odd_order);
|
|
Packit |
67cb25 |
evec = gsl_matrix_submatrix(work->evec, 0, 0, odd_order, odd_order);
|
|
Packit |
67cb25 |
gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
|
|
Packit |
67cb25 |
gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii=0; ii
|
|
Packit |
67cb25 |
bb[2*ii+1] = gsl_vector_get(&eval.vector, ii);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (ii = order_min ; ii <= order_max ; ii++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result_array[ii - order_min] = bb[ii];
|
|
Packit |
67cb25 |
}
|
|
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_a(int order, double qq)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_mathieu_a_e(order, qq, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_mathieu_b(int order, double qq)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_mathieu_b_e(order, qq, &result));
|
|
Packit |
67cb25 |
}
|