Go to the documentation of this file.
85 #define STRUVE_MAXITER 10000
87 #define SUM_TINY 1e-100
88 #define GOOD_EPS 1e-12
89 #define ACCEPTABLE_EPS 1e-7
90 #define ACCEPTABLE_ATOL 1e-300
92 #define MIN(a, b) ((a) < (b) ? (a) : (b))
100 static double struve_hl(
double v,
double x,
int is_h);
114 double value[4], err[4], tmp;
120 tmp = (
n % 2 == 0) ? -1 : 1;
140 if (
n == -
v - 0.5 &&
n > 0) {
142 return (
n % 2 == 0 ? 1 : -1) *
bessel_j(
n + 0.5,
z);
145 return iv(
n + 0.5,
z);
150 if (
z >= 0.7*
v + 12) {
179 if (err[1] < err[
n])
n = 1;
180 if (err[2] < err[
n])
n = 2;
210 double term, sum, maxterm, scaleexp, tmp;
211 double2 cterm, csum, cdiv,
z2, c2v, ctmp;
221 if (tmp < -600 || tmp > 600) {
244 cdiv =
dd_mul(cdiv, ctmp);
248 cterm =
dd_div(cterm, cdiv);
250 csum =
dd_add(csum, cterm);
255 if (
fabs(term) > maxterm) {
256 maxterm =
fabs(term);
263 *err =
fabs(term) +
fabs(maxterm) * 1
e-22;
266 sum *=
exp(scaleexp);
267 *err *=
exp(scaleexp);
270 if (sum == 0 && term == 0 &&
v < 0 && !is_h) {
287 double term, cterm, sum, maxterm;
303 cterm *=
z/2 / (
n + 1);
306 term = cterm *
iv(
n +
v + 0.5,
z) / (
n + 0.5);
307 cterm *= -
z/2 / (
n + 1);
310 if (
fabs(term) > maxterm) {
311 maxterm =
fabs(term);
318 *err =
fabs(term) +
fabs(maxterm) * 1
e-16;
321 *err += 1
e-300 *
fabs(cterm);
334 double term, sum, maxterm;
371 for (
n = 0;
n < maxiter; ++
n) {
372 term *= sgn * (1 + 2*
n) * (1 + 2*
n - 2*
v) / (
z*
z);
374 if (
fabs(term) > maxterm) {
375 maxterm =
fabs(term);
394 *err =
fabs(term) +
fabs(maxterm) * 1
e-16;
402 return cbesy_wrap_real(
v,
x);
407 return cbesj_wrap_real(
v,
x);
Namespace containing all symbols from the Eigen library.
static double bessel_y(double v, double x)
static double struve_hl(double v, double x, int is_h)
double gammasgn(double x)
Array< double, 1, 3 > e(1./3., 0.5, 2.)
static double2 dd_mul(const double2 a, const double2 b)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
static double bessel_j(double v, double x)
const EIGEN_DEVICE_FUNC LogReturnType log() const
static double2 dd_div(const double2 a, const double2 b)
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
double struve_bessel_series(double v, double z, int is_h, double *err)
double struve_asymp_large_z(double v, double z, int is_h, double *err)
static double dd_to_double(const double2 a)
static double2 dd_add(const double2 a, const double2 b)
double struve_h(double v, double z)
double struve_power_series(double v, double x, int is_h, double *err)
double struve_l(double v, double z)
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
static double2 dd_create_d(double hi)
Array< int, Dynamic, 1 > v
double iv(double v, double x)
Jet< T, N > sqrt(const Jet< T, N > &f)
gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:05:01