[9c461c7] | 1 | /** |
---|
[4f2478e] | 2 | * Spherical Bessel function 3*j1(x)/x |
---|
[9c461c7] | 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 |
---|
[e6f1410] | 7 | * using double precision that are the source. |
---|
[9c461c7] | 8 | */ |
---|
[ba32cdd] | 9 | double sph_j1c(double q); |
---|
[9c461c7] | 10 | |
---|
[e6f1410] | 11 | // The choice of the number of terms in the series and the cutoff value for |
---|
| 12 | // switching between series and direct calculation depends on the numeric |
---|
| 13 | // precision. |
---|
| 14 | // |
---|
| 15 | // Point where direct calculation reaches machine precision: |
---|
| 16 | // |
---|
| 17 | // single machine precision eps 3e-8 at qr=1.1 ** |
---|
| 18 | // double machine precision eps 4e-16 at qr=1.1 |
---|
| 19 | // |
---|
| 20 | // Point where Taylor series reaches machine precision (eps), where taylor |
---|
| 21 | // series matches direct calculation (cross) and the error at that point: |
---|
| 22 | // |
---|
| 23 | // prec n eps cross error |
---|
| 24 | // single 3 0.28 0.4 6.2e-7 |
---|
| 25 | // single 4 0.68 0.7 2.3e-7 |
---|
| 26 | // single 5 1.18 1.2 7.5e-8 |
---|
| 27 | // double 3 0.01 0.03 2.3e-13 |
---|
| 28 | // double 4 0.06 0.1 3.1e-14 |
---|
| 29 | // double 5 0.16 0.2 5.0e-15 |
---|
| 30 | // |
---|
| 31 | // ** Note: relative error on single precision starts increase on the direct |
---|
| 32 | // method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3. This should be |
---|
| 33 | // safe for the sans range, with objects of 100 nm supported to a q of 0.1 |
---|
| 34 | // while maintaining 5 digits of precision. For usans/sesans, the objects |
---|
| 35 | // are larger but the q is smaller, so again it should be fine. |
---|
| 36 | // |
---|
| 37 | // See explore/sph_j1c.py for code to explore these ranges. |
---|
[9c461c7] | 38 | |
---|
[e6f1410] | 39 | // Use 4th order series |
---|
[c437dbb] | 40 | #if FLOAT_SIZE>4 |
---|
[e6f1410] | 41 | #define SPH_J1C_CUTOFF 0.1 |
---|
[c437dbb] | 42 | #else |
---|
[e6f1410] | 43 | #define SPH_J1C_CUTOFF 0.7 |
---|
[c437dbb] | 44 | #endif |
---|
[9c461c7] | 45 | |
---|
[e6f1410] | 46 | double sph_j1c(double q) |
---|
| 47 | { |
---|
| 48 | if (q < SPH_J1C_CUTOFF) { |
---|
| 49 | const double q2 = q*q; |
---|
| 50 | return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.))))); |
---|
[4f2478e] | 51 | } else { |
---|
[e6f1410] | 52 | double sin_q, cos_q; |
---|
| 53 | SINCOS(q, sin_q, cos_q); |
---|
| 54 | return 3.0*(sin_q/q - cos_q)/(q*q); |
---|
[4f2478e] | 55 | } |
---|
[9c461c7] | 56 | } |
---|