Changes in / [5a483877:7d4b2ae] in sasmodels

Files:

• explore/J1c.py (deleted)
• explore/log_gamma.py (deleted)
• explore/rpa.py (deleted)
• explore/sinc.py (deleted)
• explore/sph_j1c.py (deleted)
• sasmodels/models/fuzzy_sphere.py (deleted)
• sasmodels/models/img/fuzzy_sphere.jpg (deleted)
• sasmodels/models/lib/lanczos_gamma.c (modified) (1 diff)
• sasmodels/models/lib/sph_j1c.c (modified) (1 diff)
• test/rpa.py (added)

sasmodels/models/lib/lanczos_gamma.c

@@ -2 +2 @@
 double lanczos_gamma(double q)
 {
-    // Lanczos approximation to the Log Gamma function.
+    // Lanczos approximation to the Gamma function.
 
     double x,y,tmp,ser;

sasmodels/models/lib/sph_j1c.c

@@ -5 +5 @@
 * Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but
 * in this case it is likely cancellation errors in the original expression
-* using double precision that are the source.
+* using double precision that are the source.  Single precision only
+* requires the first 3 terms.  Double precision requires the 4th term.
+* The fifth term is not needed, and is commented out.
+* Taylor expansion:
+*      1.0 + q2*(-3./30. + q2*(3./840.))+ q2*(-3./45360. + q2*(3./3991680.))))
 */
-
-// The choice of the number of terms in the series and the cutoff value for
-// switching between series and direct calculation depends on the numeric
-// precision.
-//
-// Point where direct calculation reaches machine precision:
-//
-//   single machine precision eps 3e-8 at qr=1.1 **
-//   double machine precision eps 4e-16 at qr=1.1
-//
-// Point where Taylor series reaches machine precision (eps), where taylor
-// series matches direct calculation (cross) and the error at that point:
-//
-//   prec   n eps  cross  error
-//   single 3 0.28  0.4   6.2e-7
-//   single 4 0.68  0.7   2.3e-7
-//   single 5 1.18  1.2   7.5e-8
-//   double 3 0.01  0.03  2.3e-13
-//   double 4 0.06  0.1   3.1e-14
-//   double 5 0.16  0.2   5.0e-15
-//
-// ** Note: relative error on single precision starts increase on the direct
-// method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3.  This should be
-// safe for the sans range, with objects of 100 nm supported to a q of 0.1
-// while maintaining 5 digits of precision.  For usans/sesans, the objects
-// are larger but the q is smaller, so again it should be fine.
-//
-// See explore/sph_j1c.py for code to explore these ranges.
-
-// Use 4th order series
-#if FLOAT_SIZE>4
-#define SPH_J1C_CUTOFF 0.1
-#else
-#define SPH_J1C_CUTOFF 0.7
-#endif
 
 double sph_j1c(double q);
 double sph_j1c(double q)
 {
-    if (q < SPH_J1C_CUTOFF) {
-        const double q2 = q*q;
-        return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.)))));
+    const double q2 = q*q;
+    double sin_q, cos_q;
+
+    SINCOS(q, sin_q, cos_q);
+
+#if FLOAT_SIZE>4
+// For double precision, choose a cutoff of 0.18, which is the lower limit
+// for the trig expression to 14 digits.  If only 12 digits are needed, then
+// only 4 terms of the Taylor expansion are needed.
+#define CUTOFF 0.18
+#else
+// For single precision, choose a cutoff halfway between the single precision
+// lower limit for the trig expression of 1.03 and the upper limit of 1.3
+// for the Taylor series expansion with 5 terms (or 9 if you count the zeros).
+// For 5 digits of precision, you can drop two terms of the taylor series
+// and choose a cutoff of 0.3.
+#define CUTOFF 1.15
+#endif
+
+    const double bessel = (q < CUTOFF)
+        ? (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))))
+        : 3.0*(sin_q/q - cos_q)/q2;
+    return bessel;
+
+/*
+    // Code to test various expressions
+    // tested using sympy.mpmath.mp with
+    //    mp.dps = 50 (which is good for q down to 1e-6)
+    //    def j1c(x): return 3*(mp.sin(x)/x - mp.cos(x))/x**2
+    double ret;
+    if (sizeof(q2) > 8) {
+
+        ret = 3.0*(sinl(q)/q - cosl(q))/q2;
+printf("%.15Lf %.15Lg\n", q, ret);
+    //} else if (q < 0.384038453352533) {
+    //} else if (q < 0.18) {
+    } else if (q < 18e-10) {
+        // NB: all are good to 5 digits below 0.18f
+        // good is defined as 14 digits for double and 7 for single
+        //ret = 1.0 + q2*q2*(3./840.) - q2*(3./30.); // good below 0.02d 0.34f
+        ret = square((1. + 3./5600.*q2*q2) - q2/20.); // good below 0.03d 0.08f
+        //ret = 1.0 - q2*(3./30.); // good below 0.02d 0.07f
+        //ret = 1.0 + q2*(-3./30. + q2*(3./840.)); // good below 0.02d 0.34f
+        //ret = 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))); // good below 0.1d 0.8f, 12 digits below 0.18d
+        //ret = 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))); // good below 0.18d 1.3f
+printf("%.15g %.15g\n", q, ret);
     } else {
-        double sin_q, cos_q;
-        SINCOS(q, sin_q, cos_q);
-        return 3.0*(sin_q/q - cos_q)/(q*q);
+        // NB: can use a cutoff of 0.1f if the goal is 5 digits rather than 7
+        ret = 3.0*(sin_q/q - cos_q)/q2; // good above 0.18d, 1.03f
+printf("%.15g %.15g\n", q, ret);
     }
+    return ret;
+*/
 }