Changeset c437dbb in sasmodels for sasmodels/models/lib/sph_j1c.c
- Timestamp:
- Mar 1, 2016 1:40:08 PM (8 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- ce0b154
- Parents:
- cf85329
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/lib/sph_j1c.c
ra3e78c3 rc437dbb 10 10 * Taylor expansion: 11 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 12 */ 16 13 … … 23 20 SINCOS(q, sin_q, cos_q); 24 21 25 const double bessel = (q < 0.384038453352533) 26 ? (1.0 + q2*(-3./30. + q2*(3./840.))) 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.))))) 27 38 : 3.0*(sin_q/q - cos_q)/q2; 28 29 39 return bessel; 30 40 31 41 /* 32 42 // 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.)))); 43 // tested using sympy.mpmath.mp with 44 // mp.dps = 50 (which is good for q down to 1e-6) 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 < 18e-10) { 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); 41 63 } else { 42 return 3.0*(sin_q/q - cos_q)/q2; 64 // NB: can use a cutoff of 0.1f if the goal is 5 digits rather than 7 65 ret = 3.0*(sin_q/q - cos_q)/q2; // good above 0.18d, 1.03f 66 printf("%.15g %.15g\n", q, ret); 43 67 } 68 return ret; 44 69 */ 45 70 }
Note: See TracChangeset
for help on using the changeset viewer.