1 | /** |
---|
2 | * Spherical Bessel function 3*j1(x)/x |
---|
3 | * |
---|
4 | * Used for low q to avoid cancellation error. |
---|
5 | * Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but |
---|
6 | * in this case it is likely cancellation errors in the original expression |
---|
7 | * using double precision that are the source. Single precision only |
---|
8 | * requires the first 3 terms. Double precision requires the 4th term. |
---|
9 | * The fifth term is not needed, and is commented out. |
---|
10 | * Taylor expansion: |
---|
11 | * 1.0 + q2*(-3./30. + q2*(3./840.))+ q2*(-3./45360. + q2*(3./3991680.)))) |
---|
12 | */ |
---|
13 | |
---|
14 | double sph_j1c(double q); |
---|
15 | double sph_j1c(double q) |
---|
16 | { |
---|
17 | const double q2 = q*q; |
---|
18 | double sin_q, cos_q; |
---|
19 | |
---|
20 | SINCOS(q, sin_q, cos_q); |
---|
21 | |
---|
22 | #if FLOAT_SIZE>4 |
---|
23 | // For double precision, choose a cutoff of 0.18, which is the lower limit |
---|
24 | // for the trig expression to 14 digits. If only 12 digits are needed, then |
---|
25 | // only 4 terms of the Taylor expansion are needed. |
---|
26 | #define CUTOFF 0.18 |
---|
27 | #else |
---|
28 | // For single precision, choose a cutoff halfway between the single precision |
---|
29 | // lower limit for the trig expression of 1.03 and the upper limit of 1.3 |
---|
30 | // for the Taylor series expansion with 5 terms (or 9 if you count the zeros). |
---|
31 | // For 5 digits of precision, you can drop two terms of the taylor series |
---|
32 | // and choose a cutoff of 0.3. |
---|
33 | #define CUTOFF 1.15 |
---|
34 | #endif |
---|
35 | |
---|
36 | const double bessel = (q < CUTOFF) |
---|
37 | ? (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.))))) |
---|
38 | : 3.0*(sin_q/q - cos_q)/q2; |
---|
39 | return bessel; |
---|
40 | |
---|
41 | /* |
---|
42 | // Code to test various expressions |
---|
43 | // tested using sympy.mpmath.mp with |
---|
44 | // mp.dps = 50 (which is good for q down to 1e-6) |
---|
45 | // def j1c(x): return 3*(mp.sin(x)/x - mp.cos(x))/x**2 |
---|
46 | double ret; |
---|
47 | if (sizeof(q2) > 8) { |
---|
48 | |
---|
49 | ret = 3.0*(sinl(q)/q - cosl(q))/q2; |
---|
50 | printf("%.15Lf %.15Lg\n", q, ret); |
---|
51 | //} else if (q < 0.384038453352533) { |
---|
52 | //} else if (q < 0.18) { |
---|
53 | } else if (q < 18e-10) { |
---|
54 | // NB: all are good to 5 digits below 0.18f |
---|
55 | // good is defined as 14 digits for double and 7 for single |
---|
56 | //ret = 1.0 + q2*q2*(3./840.) - q2*(3./30.); // good below 0.02d 0.34f |
---|
57 | ret = square((1. + 3./5600.*q2*q2) - q2/20.); // good below 0.03d 0.08f |
---|
58 | //ret = 1.0 - q2*(3./30.); // good below 0.02d 0.07f |
---|
59 | //ret = 1.0 + q2*(-3./30. + q2*(3./840.)); // good below 0.02d 0.34f |
---|
60 | //ret = 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))); // good below 0.1d 0.8f, 12 digits below 0.18d |
---|
61 | //ret = 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))); // good below 0.18d 1.3f |
---|
62 | printf("%.15g %.15g\n", q, ret); |
---|
63 | } else { |
---|
64 | // NB: can use a cutoff of 0.1f if the goal is 5 digits rather than 7 |
---|
65 | ret = 3.0*(sin_q/q - cos_q)/q2; // good above 0.18d, 1.03f |
---|
66 | printf("%.15g %.15g\n", q, ret); |
---|
67 | } |
---|
68 | return ret; |
---|
69 | */ |
---|
70 | } |
---|