36 #ifndef VIGRA_BESSEL_HXX
37 #define VIGRA_BESSEL_HXX
39 #include "mathutil.hxx"
42 #include <boost/math/special_functions/bessel.hpp>
54 int msta1(REAL x,
int mp)
66 nn = int(n1-(n1-n0)/(1.0-f0/f1));
79 int msta2(REAL x,
int n,
int mp)
81 double a0,ejn,hmp,f0,f1,f,obj;
104 nn = int(n1-(n1-n0)/(1.0-f0/f1));
133 template <
class REAL>
134 void bessjyn(
int n, REAL x,
int &nm,
double *jn,
double *yn)
136 double t1,t2,f,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv;
137 double ec,bs,byk,p0,p1,q0,q1;
138 static double a[] = {
139 -0.7031250000000000e-1,
143 static double b[] = {
144 0.7324218750000000e-1,
147 -2.438052969955606e1};
148 static double a1[] = {
153 static double b1[] = {
157 2.724882731126854e1};
171 if (x <= 300.0 || n > (
int)(0.9*x))
187 f = 2.0*(k+1.0)/x*f1 - f2;
190 if ((k == 2*(
int)(k/2)) && (k != 0))
193 su += (-1)*((k & 2)-1)*f/(
double)k;
197 sv += (-1)*((k & 2)-1)*(
double)k*f/(k*k-1.0);
207 ec =
std::log(0.5*x) + M_EULER_GAMMA;
208 by0 = M_2_PI*(ec*jn[0]-4.0*su/s0);
210 by1 = M_2_PI*((ec-1.0)*jn[1]-jn[0]/x-4.0*sv/s0);
220 p0 += a[k]*std::pow(x,-2*k-2);
221 q0 += b[k]*std::pow(x,-2*k-3);
233 p1 += a1[k]*std::pow(x,-2*k-2);
234 q1 += b1[k]*std::pow(x,-2*k-3);
242 bjk = 2.0*(k-1.0)*bj1/x-bj0;
250 byk = 2.0*(k-1.0)*by1/x-by0;
277 throw std::domain_error(
"besselJ(n, x): x cannot be negative");
279 return n == 0 ? 1.0 : 0.0;
280 #if defined(HasBoostMath)
281 return boost::math::cyl_bessel_j((
double)n, x);
282 #elif defined(__GNUC__)
284 #elif defined(_MSC_VER)
287 int an =
abs(n), nr = n, s = an+2;
289 detail::bessjyn(an, x, nr, &t[0], &t[s]);
313 throw std::domain_error(
"besselY(n, x): x cannot be negative");
315 return -std::numeric_limits<double>::infinity();
316 #if defined(HasBoostMath)
317 return boost::math::cyl_neumann((
double)n, x);
318 #elif defined(__GNUC__)
320 #elif defined(_MSC_VER)
323 int an =
abs(n), nr = n, s = an+2;
325 detail::bessjyn(an, x, nr, &t[0], &t[s]);
326 if(n < 0.0 &&
odd(n))
337 #endif // VIGRA_BESSEL_HXX