Changeset eb2946f in sasmodels for sasmodels/models/lib/sas_J1.c

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

File:

• sasmodels/models/lib/sas_J1.c (modified) (6 diffs)

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