78 static void ikv_temme(
double v,
double x,
double *Iv,
double *Kv);
80 double iv(
double v,
double x)
104 if (
v != 2.0 *
floor(
v / 2.0)) {
148 double sum, term, prefactor, factor;
153 if (prefactor == INFINITY) {
163 factor = (
mu - (2 * k - 1) * (2 * k - 1)) / (8 *
x) / k;
173 return sum * prefactor;
196 #define N_UFACTORS 11
197 #define N_UFACTOR_TERMS 31
199 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
200 0, 0, 0, 0, 0, 0, 0, 1},
201 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
202 0, 0, 0, 0, -0.20833333333333334, 0.0, 0.125, 0.0},
203 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
204 0, 0.3342013888888889, 0.0, -0.40104166666666669, 0.0, 0.0703125, 0.0,
206 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
207 -1.0258125964506173, 0.0, 1.8464626736111112, 0.0,
208 -0.89121093750000002, 0.0, 0.0732421875, 0.0, 0.0, 0.0},
209 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
210 4.6695844234262474, 0.0, -11.207002616222995, 0.0, 8.78912353515625,
211 0.0, -2.3640869140624998, 0.0, 0.112152099609375, 0.0, 0.0, 0.0, 0.0},
212 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -28.212072558200244, 0.0,
213 84.636217674600744, 0.0, -91.818241543240035, 0.0, 42.534998745388457,
214 0.0, -7.3687943594796312, 0.0, 0.22710800170898438, 0.0, 0.0, 0.0,
216 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 212.5701300392171, 0.0,
217 -765.25246814118157, 0.0, 1059.9904525279999, 0.0,
218 -699.57962737613275, 0.0, 218.19051174421159, 0.0,
219 -26.491430486951554, 0.0, 0.57250142097473145, 0.0, 0.0, 0.0, 0.0,
221 {0, 0, 0, 0, 0, 0, 0, 0, 0, -1919.4576623184068, 0.0,
222 8061.7221817373083, 0.0, -13586.550006434136, 0.0, 11655.393336864536,
223 0.0, -5305.6469786134048, 0.0, 1200.9029132163525, 0.0,
224 -108.09091978839464, 0.0, 1.7277275025844574, 0.0, 0.0, 0.0, 0.0, 0.0,
226 {0, 0, 0, 0, 0, 0, 20204.291330966149, 0.0, -96980.598388637503, 0.0,
227 192547.0012325315, 0.0, -203400.17728041555, 0.0, 122200.46498301747,
228 0.0, -41192.654968897557, 0.0, 7109.5143024893641, 0.0,
229 -493.915304773088, 0.0, 6.074042001273483, 0.0, 0.0, 0.0, 0.0, 0.0,
231 {0, 0, 0, -242919.18790055133, 0.0, 1311763.6146629769, 0.0,
232 -2998015.9185381061, 0.0, 3763271.2976564039, 0.0,
233 -2813563.2265865342, 0.0, 1268365.2733216248, 0.0,
234 -331645.17248456361, 0.0, 45218.768981362737, 0.0,
235 -2499.8304818112092, 0.0, 24.380529699556064, 0.0, 0.0, 0.0, 0.0, 0.0,
237 {3284469.8530720375, 0.0, -19706819.11843222, 0.0, 50952602.492664628,
238 0.0, -74105148.211532637, 0.0, 66344512.274729028, 0.0,
239 -37567176.660763353, 0.0, 13288767.166421819, 0.0,
240 -2785618.1280864552, 0.0, 308186.40461266245, 0.0,
241 -13886.089753717039, 0.0, 110.01714026924674, 0.0, 0.0, 0.0, 0.0, 0.0,
242 0.0, 0.0, 0.0, 0.0, 0.0}
250 double *i_value,
double *k_value)
252 double i_prefactor, k_prefactor;
253 double t, t2, eta,
z;
254 double i_sum, k_sum, term, divisor;
287 for (k = 1; k <
n; k += 2) {
297 k_sum += (
n % 2 == 0) ? term : -term;
316 if (k_value !=
NULL) {
318 *k_value = k_prefactor * k_sum;
321 if (i_value !=
NULL) {
323 *i_value = i_prefactor * i_sum;
327 *i_value = (i_prefactor * i_sum
341 #define BOOST_ASSERT(a) assert(a)
343 #define BOOST_ASSERT(a)
354 double f,
h,
p,
q, coef, sum, sum1, tolerance;
377 gamma2 = (2 + gp + gm) *
c / 2;
380 p = (gp + 1) / (2 *
b);
381 q = (1 + gm) *
b / 2;
390 for (k = 1; k <
MAXITER; k++) {
391 f = (k *
f +
p +
q) / (k * k -
v *
v);
395 coef *=
x *
x / (4 * k);
398 if (
fabs(coef *
f) <
fabs(sum) * tolerance) {
430 tiny = 1 /
sqrt(DBL_MAX);
433 for (k = 1; k <
MAXITER; k++) {
465 static int CF2_ik(
double v,
double x,
double *Kv,
double *Kv1)
468 double S,
C,
Q,
D,
f,
a,
b,
q,
delta, tolerance, current, prev;
491 for (k = 2; k <
MAXITER; k++) {
500 q = (prev - (
b - 2) * current) /
a;
517 *Kv1 = *Kv * (0.5f +
v +
x + (
v *
v - 0.25f) *
f) /
x;
532 static void ikv_temme(
double v,
double x,
double *Iv_p,
double *Kv_p)
536 double u, Iv, Kv, Kv1, Ku, Ku1, fv;
537 double W, current, prev, next;
567 Iv = (
v == 0) ? 1 : 0;
576 if (reflect && (kind &
need_i)) {
577 double z = (u +
n % 2);
579 Iv =
sin((
double)
M_PI *
z) == 0 ? Iv : INFINITY;
580 if (Iv == INFINITY || Iv == -INFINITY) {
603 for (k = 1; k <=
n; k++) {
604 next = 2 * (u + k) * current /
x + prev;
611 double lim = (4 *
v *
v + 10) / (8 *
x);
616 if ((lim <
MACHEP * 10) && (
x > 100)) {
628 Iv =
W / (Kv * fv + Kv1);
636 double z = (u +
n % 2);