Changeset 3f8584a2 in sasmodels for sasmodels/models/lib/sas_J0.c

Timestamp:

Jul 27, 2016 3:53:07 AM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

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:

c047acf

Parents:

56547a8

Message:

clean separation between float/double bessel functions

File:

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

sasmodels/models/lib/sas_J0.c

@@ -53 +53 @@
  * except YP, YQ which are designed for absolute error. */
 
+#if FLOAT_SIZE>4
+//Cephes double precission
 
  constant double PPJ0[8] = {
@@ -144 +146 @@
   };
 
- constant double MOJ0[8] = {
+double j0(double x);
+double j0(double x)
+{
+    double w, z, p, q, xn;
+
+    //const double TWOOPI = 6.36619772367581343075535E-1;
+    const double SQ2OPI = 7.9788456080286535587989E-1;
+    const double PIO4 = 7.85398163397448309616E-1;
+
+    const double DR1 = 5.78318596294678452118E0;
+    const double DR2 = 3.04712623436620863991E1;
+
+
+    if( x < 0 )
+        x = -x;
+
+    if( x <= 5.0 ) {
+        z = x * x;
+        if( x < 1.0e-5 )
+            return( 1.0 - z/4.0 );
+
+        p = (z - DR1) * (z - DR2);
+        p = p * polevl( z, RPJ0, 3)/p1evl( z, RQJ0, 8 );
+        return( p );
+    }
+
+    w = 5.0/x;
+    q = 25.0/(x*x);
+    p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 );
+    q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 );
+    xn = x - PIO4;
+
+    double sn, cn;
+    SINCOS(xn, sn, cn);
+    p = p * cn - w * q * sn;
+
+    return( p * SQ2OPI / sqrt(x) );
+}
+#else
+
+ constant float MOJ0[8] = {
         -6.838999669318810E-002,
         1.864949361379502E-001,
@@ -155 +197 @@
   };
 
- constant double PHJ0[8] = {
+ constant float PHJ0[8] = {
         3.242077816988247E+001,
         -3.630592630518434E+001,
@@ -166 +208 @@
   };
 
- constant double JPJ0[8] = {
+ constant float JPJ0[8] = {
         -6.068350350393235E-008,
         6.388945720783375E-006,
@@ -177 +219 @@
  };
 
-double sas_J0(double x);
-double sas_J0(double x)
+//Cephes single precission
+float j0f(float x);
+float j0f(float x)
 {
-
-//Cephes single precission
-#if FLOAT_SIZE>4
-    double w, z, p, q, xn;
-
-    //const double TWOOPI = 6.36619772367581343075535E-1;
-    const double SQ2OPI = 7.9788456080286535587989E-1;
-    const double PIO4 = 7.85398163397448309616E-1;
-
-    const double DR1 = 5.78318596294678452118E0;
-    const double DR2 = 3.04712623436620863991E1;
-
-
-    if( x < 0 )
-            x = -x;
-
-    if( x <= 5.0 ) {
-            z = x * x;
-            if( x < 1.0e-5 )
-                    return( 1.0 - z/4.0 );
-
-            p = (z - DR1) * (z - DR2);
-            p = p * polevl( z, RPJ0, 3)/p1evl( z, RQJ0, 8 );
-            return( p );
-        }
-
-    w = 5.0/x;
-    q = 25.0/(x*x);
-    p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 );
-    q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 );
-    xn = x - PIO4;
-
-    double sn, cn;
-    SINCOS(xn, sn, cn);
-    p = p * cn - w * q * sn;
-
-    return( p * SQ2OPI / sqrt(x) );
-//Cephes single precission
-#else
-    double xx, w, z, p, q, xn;
+    float xx, w, z, p, q, xn;
 
     //const double YZ1 =  0.43221455686510834878;
     //const double YZ2 = 22.401876406482861405;
     //const double YZ3 = 64.130620282338755553;
-    const double DR1 =  5.78318596294678452118;
-    const double PIO4F = 0.7853981633974483096;
+    const float DR1 =  5.78318596294678452118;
+    const float PIO4F = 0.7853981633974483096;
 
     if( x < 0 )
-            xx = -x;
+        xx = -x;
     else
-            xx = x;
+        xx = x;
 
     if( x <= 2.0 ) {
-            z = xx * xx;
-            if( x < 1.0e-3 )
-                    return( 1.0 - 0.25*z );
-
-            p = (z-DR1) * polevl( z, JPJ0, 4);
-            return( p );
-        }
+        z = xx * xx;
+        if( x < 1.0e-3 )
+            return( 1.0 - 0.25*z );
+
+        p = (z-DR1) * polevl( z, JPJ0, 4);
+        return( p );
+    }
 
     q = 1.0/x;
@@ -249 +253 @@
     p = p * cos(xn + xx);
     return(p);
+}
 #endif
 
-}
-
-
+#if FLOAT_SIZE>4
+#define sas_J0 j0
+#else
+#define sas_J0 j0f
+#endif