Changeset f0afad2 in sasmodels for sasmodels/models/poly_gauss_coil.py

Timestamp:

Oct 18, 2016 10:54:57 PM (8 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:

b829b16

Parents:

6d96b66

Message:

poly_gauss_coil: improved accuracy

File:

• sasmodels/models/poly_gauss_coil.py (modified) (2 diffs)

sasmodels/models/poly_gauss_coil.py

@@ -59 +59 @@
 """
 
-from numpy import inf, exp, power
+import numpy as np
+from numpy import inf, expm1, power
 
 name = "poly_gauss_coil"
@@ -83 +84 @@
     # pylint: disable = missing-docstring
     u = polydispersity - 1.0
-    z = (q*rg)**2 / (1.0 + 2.0*u)
+    z = q**2 * (rg**2 / (1.0 + 2.0*u))
+
     # need to trap the case of the polydispersity being 1 (ie, monodisperse!)
     if polydispersity == 1.0:
-        inten = i_zero * 2.0 * (exp(-z) + z - 1.0)
+        result = 2.0 * (expm1(-z) + z)
+        index = q != 0.
+        result[index] /= z[index]**2
+        result[~index] = 1.0
     else:
-        inten = i_zero * 2.0 * (power(1.0 + u*z, -1.0/u) + z - 1.0) / (1.0 + u)
-    index = q != 0.
-    inten[~index] = i_zero
-    inten[index] /= z[index]**2
-    return inten
+        # Taylor series around z=0 of (2*(1+uz)^(-1/u) + z - 1) / (z^2(u+1))
+        p = [
+            #(-1 - 20*u - 155*u**2 - 580*u**3 - 1044*u**4 - 720*u**5) / 2520.,
+            #( 1 + 14*u + 71*u**2 + 154*u**3 + 120*u**4) / 360.,
+            #(-1 - 9*u - 26*u**2 - 24*u**3) / 60.,
+            ( 1 + 5*u + 6*u**2) / 12.,
+            (-1 - 2*u) / 3.,
+            ( 1 ),
+            ]
+        result = 2.0 * (power(1.0 + u*z, -1.0/u) + z - 1.0) / (1.0 + u)
+        index = z > 1e-4
+        result[index] /= z[index]**2
+        result[~index] = np.polyval(p, z[~index])
+    return i_zero * result
 Iq.vectorized = True  # Iq accepts an array of q values