36#ifndef VIGRA_BESSEL_HXX
37#define VIGRA_BESSEL_HXX
39#include "mathutil.hxx"
42#include <boost/math/special_functions/bessel.hpp>
54int msta1(REAL x,
int mp)
61 f0 = 0.5*std::log10(6.28*n0) - n0*std::log10(1.36*a0/n0)-mp;
63 f1 = 0.5*std::log10(6.28*n1) - n1*std::log10(1.36*a0/n1)-mp;
66 nn = int(n1-(n1-n0)/(1.0-f0/f1));
67 f = 0.5*std::log10(6.28*nn) - nn*std::log10(1.36*a0/nn)-mp;
79int msta2(REAL x,
int n,
int mp)
81 double a0,ejn,hmp,f0,f1,f,obj;
86 ejn = 0.5*std::log10(6.28*n) - n*std::log10(1.36*a0/n);
99 f0 = 0.5*std::log10(6.28*n0) - n0*std::log10(1.36*a0/n0)-obj;
101 f1 = 0.5*std::log10(6.28*n1) - n1*std::log10(1.36*a0/n1)-obj;
104 nn = int(n1-(n1-n0)/(1.0-f0/f1));
105 f = 0.5*std::log10(6.28*nn) - nn*std::log10(1.36*a0/nn)-obj;
134void 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;
139 -0.7031250000000000e-1,
144 0.7324218750000000e-1,
147 -2.438052969955606e1};
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);
223 cu = std::sqrt(M_2_PI/x);
224 bj0 = cu*(p0*std::cos(t1)-q0*std::sin(t1));
225 by0 = cu*(p0*std::sin(t1)+q0*std::cos(t1));
233 p1 += a1[k]*std::pow(x,-2*k-2);
234 q1 += b1[k]*std::pow(x,-2*k-3);
236 bj1 = cu*(p1*std::cos(t2)-q1*std::sin(t2));
237 by1 = cu*(p1*std::sin(t2)+q1*std::cos(t2));
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))
double besselJ(int n, double x)
Bessel function of the first kind.
Definition: bessel.hxx:274
double besselY(int n, double x)
Bessel function of the second kind.
Definition: bessel.hxx:310
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
bool odd(int t)
Check if an integer is odd.