source: sasmodels/sasmodels/models/lib/sph_j1c.c @ e6f1410

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since e6f1410 was e6f1410, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

examine numerical precision of J1c, sinc, j1c and log gamma

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