Changeset a1c32c2 in sasmodels for explore

Timestamp:

Dec 12, 2017 1:08:07 PM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

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

Children:

c11d09f

Parents:

cfa28d3

Message:

Use gauss-legendre integration for cross-checking against P(r)

Location:

explore

Files:

• asymint.py (modified) (1 diff)
• realspace.py (modified) (5 diffs)

explore/asymint.py

@@ -223 +223 @@
     DA, DB, DC = 2300, 21, 58
     SLDA, SLDB, SLDC = "5", "-0.3", "11.5"
+    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3
+    #SLD_SOLVENT,CONTRAST = 0, 4
     if 1: # C shortest
         B, C = C, B

explore/realspace.py

@@ -7 +7 @@
 from numpy import pi, radians, sin, cos, sqrt
 from numpy.random import poisson, uniform
+from numpy.polynomial.legendre import leggauss
 from scipy.integrate import simps
 from scipy.special import j1 as J1
@@ -362 +363 @@
 
 def cylinder_Iq(q, radius, length):
-    height = length/2
+    z, w = leggauss(76)
+    cos_alpha = (z+1)/2
+    sin_alpha = sqrt(1.0 - cos_alpha**2)
     Iq = np.empty_like(q)
-    theta = np.linspace(0.0, pi/2, 500)
     for k, qk in enumerate(q):
-        qab, qc = qk*sin(theta), qk*cos(theta)
-        Iq[k] = simps((j0(qc*height)*sas_2J1x_x(qab*radius))**2 * sin(theta))
+        qab, qc = qk*sin_alpha, qk*cos_alpha
+        Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2)
+        Iq[k] = np.sum(w*Fq**2)
     Iq = Iq/Iq[0]
     return Iq
@@ -376 +379 @@
 
 def csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core):
+    z, w = leggauss(76)
+
     sld_solvent = 0
     overlapping = False
@@ -382 +387 @@
     tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc
 
-    Iq = np.zeros_like(q)
-    for cos_alpha in np.linspace(0.0, 1.0, 100):
+    outer_sum = np.zeros_like(q)
+    for cos_alpha, outer_w in zip((z+1)/2, w):
         sin_alpha = sqrt(1.0-cos_alpha*cos_alpha)
         qc = q*cos_alpha
         siC = c*j0(c*qc/2)
         siCt = tC*j0(tC*qc/2)
-        for beta in np.linspace(0.0, pi/2, 100):
+        inner_sum = np.zeros_like(q)
+        for beta, inner_w in zip((z + 1)*pi/4, w):
             qa, qb = q*sin_alpha*sin(beta), q*sin_alpha*cos(beta)
             siA = a*j0(a*qa/2)
@@ -395 +401 @@
             siBt = tB*j0(tB*qb/2)
             if overlapping:
-                Iq += (dr0*siA*siB*siC
-                       + drA*(siAt-siA)*siB*siC
-                       + drB*siAt*(siBt-siB)*siC
-                       + drC*siAt*siBt*(siCt-siC))**2
+                Fq = (dr0*siA*siB*siC
+                      + drA*(siAt-siA)*siB*siC
+                      + drB*siAt*(siBt-siB)*siC
+                      + drC*siAt*siBt*(siCt-siC))
             else:
-                Iq += (dr0*siA*siB*siC
-                       + drA*(siAt-siA)*siB*siC
-                       + drB*siA*(siBt-siB)*siC
-                       + drC*siA*siB*(siCt-siC))**2
+                Fq = (dr0*siA*siB*siC
+                      + drA*(siAt-siA)*siB*siC
+                      + drB*siA*(siBt-siB)*siC
+                      + drC*siA*siB*(siCt-siC))
+            inner_sum += inner_w * Fq**2
+        outer_sum += outer_w * inner_sum
+    Iq = outer_sum / 4  # = outer*um*zm*8.0/(4.0*M_PI)
     return Iq/Iq[0]