core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change
on this file since a3e78c3 was
a3e78c3,
checked in by Paul Kienzle <pkienzle@…>, 8 years ago
|
simplify calculation for sph_j1c
|
-
Property mode set to
100644
|
File size:
1.5 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. 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 | * Expression returned from Herbie (herbie.uwpise.org/demo): |
---|
13 | * const double t = ((1. + 3.*q2*q2/5600.) - q2/20.); |
---|
14 | * return t*t; |
---|
15 | */ |
---|
16 | |
---|
17 | double sph_j1c(double q); |
---|
18 | double sph_j1c(double q) |
---|
19 | { |
---|
20 | const double q2 = q*q; |
---|
21 | double sin_q, cos_q; |
---|
22 | |
---|
23 | SINCOS(q, sin_q, cos_q); |
---|
24 | |
---|
25 | const double bessel = (q < 0.384038453352533) |
---|
26 | ? (1.0 + q2*(-3./30. + q2*(3./840.))) |
---|
27 | : 3.0*(sin_q/q - cos_q)/q2; |
---|
28 | |
---|
29 | return bessel; |
---|
30 | |
---|
31 | /* |
---|
32 | // Code to test various expressions |
---|
33 | if (sizeof(q2) > 4) { |
---|
34 | return 3.0*(sin_q/q - cos_q)/q2; |
---|
35 | } else if (q < 0.384038453352533) { |
---|
36 | //const double t = ((1. + 3.*q2*q2/5600.) - q2/20.); return t*t; |
---|
37 | return 1.0 + q2*q2*(3./840.) - q2*(3./30.); |
---|
38 | //return 1.0 + q2*(-3./30. + q2*(3./840.)); |
---|
39 | //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))); |
---|
40 | //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))); |
---|
41 | } else { |
---|
42 | return 3.0*(sin_q/q - cos_q)/q2; |
---|
43 | } |
---|
44 | */ |
---|
45 | } |
---|
Note: See
TracBrowser
for help on using the repository browser.