core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change
on this file since 13ed84c 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
|
Rev | Line | |
---|
[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 |
---|
| 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. |
---|
[4f2478e] | 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; |
---|
[9c461c7] | 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 | |
---|
[4f2478e] | 25 | const double bessel = (q < 0.384038453352533) |
---|
[a3e78c3] | 26 | ? (1.0 + q2*(-3./30. + q2*(3./840.))) |
---|
[9c461c7] | 27 | : 3.0*(sin_q/q - cos_q)/q2; |
---|
| 28 | |
---|
| 29 | return bessel; |
---|
[4f2478e] | 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 | */ |
---|
[9c461c7] | 45 | } |
---|
Note: See
TracBrowser
for help on using the repository browser.