Changeset 237c9cf in sasmodels for explore

Timestamp:

Oct 4, 2017 1:20:48 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:

5772737

Parents:

3a220e6

Message:

flexible_cylinder: improve low q Sdebye calculation; adjust code so different branches use the same equation numbers

File:

• explore/precision.py (modified) (4 diffs)

explore/precision.py

@@ -326 +326 @@
     name="debye",
     mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4,
-    np_function=lambda x: 2*(np.exp(-x**2) + x**2 - 1)/x**4,
+    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4,
     ocl_function=make_ocl("""
     const double qsq = q*q;
-    if (qsq < -0.75*0.75) { // Pade approximation around 0
+    if (qsq < 1.0) { // Pade approximation
         const double x = qsq;
-        if (0) {
+        if (0) { // 0.36 single
             // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 4}]
             return (x*x/180. + 1.)/((1./30.*x + 1./3.)*x + 1);
-        } else if (0) {
+        } else if (0) { // 1.0 for single
             // padeapproximant[2*exp[-x^2] + x^2-1)/x^4, {x, 0, 6}]
-            // works for single precision, with qsq < 1
             const double A1=1./24., A2=1./84, A3=-1./3360;
             const double B1=3./8., B2=3./56., B3=1./336.;
             return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.);
-        } else if (0) {
+        } else if (1) { // 1.0 for single, 0.25 for double
             // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
             const double A1=1./15., A2=1./60, A3=0., A4=1./75600.;
@@ -346 +345 @@
             return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.)
                   /((((B4*x + B3)*x + B2)*x + B1)*x + 1.);
-
-        } else {
+        } else { // 1.0 for single, 0.5 for double
             // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
-            // works for double precision, with qsq < 9/16
             const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.;
             const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.;
@@ -355 +352 @@
                   /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.);
         }
-    } else if (qsq < 0.1) { // Taylor series around 0.
+    } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double
         const double x = qsq;
         const double C0 = +1.;
@@ -369 +366 @@
         return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
     } else {
-        return 2.*(exp(-qsq) + qsq - 1.)/(qsq*qsq);
+        return 2.*(expm1(-qsq) + qsq)/(qsq*qsq);
     }
     """, "sas_debye"),