Changeset eb2946f in sasmodels

Timestamp:

May 18, 2017 5:10:29 PM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

1b85b55

Parents:

452b168

Message:

improve precision of sasmodels special functions

Files:

• explore/precision.py (added)
• explore/precision/J1.py (deleted)
• explore/precision/J1c.py (deleted)
• explore/precision/cyl.py (deleted)
• explore/precision/lgamma.py (deleted)
• explore/precision/log_gamma.py (deleted)
• explore/precision/sinc.py (deleted)
• explore/precision/sph_j1c.py (deleted)
• sasmodels/models/lib/sas_3j1x_x.c (modified) (1 diff)
• sasmodels/models/lib/sas_J0.c (modified) (2 diffs)
• sasmodels/models/lib/sas_J1.c (modified) (6 diffs)
• sasmodels/models/lib/sas_erf.c (modified) (2 diffs)

sasmodels/models/lib/sas_3j1x_x.c

@@ -46 +46 @@
 double sas_3j1x_x(double q)
 {
-    if (q < SPH_J1C_CUTOFF) {
+    // 2017-05-18 PAK - support negative q
+    if (fabs(q) < SPH_J1C_CUTOFF) {
         const double q2 = q*q;
         return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.)))));

sasmodels/models/lib/sas_J0.c

@@ -236 +236 @@
         xx = x;
 
-    if( x <= 2.0 ) {
+    // 2017-05-18 PAK - support negative x
+    if( xx <= 2.0 ) {
         z = xx * xx;
-        if( x < 1.0e-3 )
+        if( xx < 1.0e-3 )
             return( 1.0 - 0.25*z );
 
@@ -245 +246 @@
     }
 
-    q = 1.0/x;
+    q = 1.0/xx;
     w = sqrt(q);
 

sasmodels/models/lib/sas_J1.c

@@ -109 +109 @@
 {
 
-    double w, z, p, q, xn;
+    double w, z, p, q, xn, abs_x, sign_x;
 
     const double Z1 = 1.46819706421238932572E1;
@@ -116 +116 @@
     const double SQ2OPI = 0.79788456080286535588;
 
-    w = x;
-    if( x < 0 )
-        w = -x;
-
-    if( w <= 5.0 ) {
-        z = x * x;
+    // 2017-05-18 PAK - mathematica and mpmath use J1(-x) = -J1(x)
+    if (x < 0) {
+        abs_x = -x;
+        sign_x = -1.0;
+    } else {
+        abs_x = x;
+        sign_x = 1.0;
+    }
+
+    if( abs_x <= 5.0 ) {
+        z = abs_x * abs_x;
         w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 );
-        w = w * x * (z - Z1) * (z - Z2);
-        return( w );
+        w = w * abs_x * (z - Z1) * (z - Z2);
+        return( sign_x * w );
     }
 
-    w = 5.0/x;
+    w = 5.0/abs_x;
     z = w * w;
-
     p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 );
     q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 );
-
-    xn = x - THPIO4;
+    xn = abs_x - THPIO4;
 
     double sn, cn;
@@ -139 +142 @@
     p = p * cn - w * q * sn;
 
-    return( p * SQ2OPI / sqrt(x) );
+    return( sign_x * p * SQ2OPI / sqrt(abs_x) );
 }
 
@@ -179 +182 @@
     };
 
-float cephes_j1f(float x)
+float cephes_j1f(float xx)
 {
 
-    float xx, w, z, p, q, xn;
+    float x, w, z, p, q, xn;
 
     const float Z1 = 1.46819706421238932572E1;
@@ -188 +191 @@
 
 
-    xx = x;
-    if( xx < 0 )
-        xx = -x;
-
-    if( xx <= 2.0 ) {
-        z = xx * xx;
-        p = (z-Z1) * xx * polevl( z, JPJ1, 4 );
-        return( p );
+    // 2017-05-18 PAK - mathematica and mpmath use J1(-x) = -J1(x)
+    x = xx;
+    if( x < 0 )
+        x = -xx;
+
+    if( x <= 2.0 ) {
+        z = x * x;
+        p = (z-Z1) * x * polevl( z, JPJ1, 4 );
+        return( xx < 0. ? -p : p );
     }
 
@@ -204 +208 @@
     w = q*q;
     xn = q * polevl( w, PH1J1, 7) - THPIO4F;
-    p = p * cos(xn + xx);
-
-    return(p);
+    p = p * cos(xn + x);
+
+    return( xx < 0. ? -p : p );
 }
 #endif

sasmodels/models/lib/sas_erf.c

@@ -165 +165 @@
         // the erf function instead and z < 1.0.
         //return (1.0 - cephes_erf(a));
-        z = x * x;
-        y = x * polevl(z, TD, 4) / p1evl(z, UD, 5);
-
-        return y;
+        // 2017-05-18 PAK - use erf(a) rather than erf(|a|)
+        z = a * a;
+        y = a * polevl(z, TD, 4) / p1evl(z, UD, 5);
+
+        return 1.0 - y;
     }
 
@@ -279 +280 @@
         //is explicit here for the case < 1.0
         //return (1.0 - sas_erf(a));
-        z = x * x;
-        y = x * polevl( z, TF, 6 );
-
-        return y;
+        // 2017-05-18 PAK - use erf(a) rather than erf(|a|)
+        z = a * a;
+        y = a * polevl( z, TF, 6 );
+
+        return 1.0 - y;
     }