Changeset e6f1410 in sasmodels for sasmodels/models/lib/sph_j1c.c
 Timestamp:
 Mar 6, 2016 8:16:29 PM (7 years ago)
 Branches:
 master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 8dca856
 Parents:
 ad90df9
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sasmodels/models/lib/sph_j1c.c
rc437dbb re6f1410 5 5 * Note that the values differ from sasview ~ 5e12 rather than 5e14, but 6 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.)))) 7 * using double precision that are the source. 12 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 3e8 at qr=1.1 ** 17 // double machine precision eps 4e16 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.2e7 24 // single 4 0.68 0.7 2.3e7 25 // single 5 1.18 1.2 7.5e8 26 // double 3 0.01 0.03 2.3e13 27 // double 4 0.06 0.1 3.1e14 28 // double 5 0.16 0.2 5.0e15 29 // 30 // ** Note: relative error on single precision starts increase on the direct 31 // method at qr=1.1, rising from 3e8 to 5e5 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 13 44 14 45 double sph_j1c(double q); 15 46 double sph_j1c(double q) 16 47 { 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 1e6) 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 < 18e10) { 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); 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.))))); 63 51 } else { 64 // NB: can use a cutoff of 0.1f if the goal is 5 digits rather than 765 ret = 3.0*(sin_q/q  cos_q)/q2; // good above 0.18d, 1.03f66 printf("%.15g %.15g\n", q, ret);52 double sin_q, cos_q; 53 SINCOS(q, sin_q, cos_q); 54 return 3.0*(sin_q/q  cos_q)/(q*q); 67 55 } 68 return ret;69 */70 56 }
Note: See TracChangeset
for help on using the changeset viewer.