Changeset c462169 in sasmodels

Timestamp:

Apr 2, 2018 10:48:10 AM (6 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:

29afc50

Parents:

b3703f5

git-author:

Paul Kienzle <pkienzle@…> (04/02/18 10:15:58)

git-committer:

Paul Kienzle <pkienzle@…> (04/02/18 10:48:10)

Message:

show u-sub integration as used in models in explore/asymint.py

File:

• explore/asymint.py (modified) (7 diffs)

explore/asymint.py

@@ -30 +30 @@
 import numpy as np
 import mpmath as mp
-from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10
+from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10, arccos
 from numpy.polynomial.legendre import leggauss
 from scipy.integrate import dblquad, simps, romb, romberg
@@ -41 +41 @@
 shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1]
 Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2]
+
+DTYPE = 'd'
 
 class MPenv:
@@ -72 +74 @@
     sas_sinx_x = staticmethod(sp.sas_sinx_x)
     pi = np.pi
-    mpf = staticmethod(float)
+    #mpf = staticmethod(float)
+    mpf = staticmethod(lambda x: np.array(x, DTYPE))
 
 SLD = 3
@@ -220 +223 @@
     #A, B, C = 4450, 14000, 47
     #A, B, C = 445, 140, 47  # integer for the sake of mpf
-    A, B, C = 6800, 114, 1380
-    DA, DB, DC = 2300, 21, 58
+    A, B, C = 114, 1380, 6800
+    DA, DB, DC = 21, 58, 2300
     SLDA, SLDB, SLDC = "5", "-0.3", "11.5"
+    ## default parameters from sasmodels
+    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2
+    ## swap A-B-C to C-B-A
+    #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA
     #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3
     #SLD_SOLVENT,CONTRAST = 0, 4
@@ -233 +240 @@
         DA, DC = DC, DA
         SLDA, SLDC = SLDC, SLDA
+    #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
     NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
     NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv)
@@ -348 +356 @@
     return n**2, Iq
 
+def gauss_quad_usub(q, n=150, dtype=DTYPE):
+    """
+    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
+
+    Use *u = sin theta* substitution, and restrict integration over a single
+    quadrant for shapes that are mirror symmetric about AB, AC and BC planes.
+
+    Note that this doesn't work for fcc/bcc paracrystals, which instead step
+    over the entire 4 pi surface uniformly in theta-phi.
+    """
+    z, w = leggauss(n)
+    cos_theta = 0.5 * (z + 1)
+    theta = arccos(cos_theta)
+    phi = pi/2*(0.5 * (z + 1))
+    Atheta, Aphi = np.meshgrid(theta, phi)
+    Aw = w[None, :] * w[:, None]
+    q, Atheta, Aphi, Aw = [np.asarray(v, dtype=dtype) for v in (q, Atheta, Aphi, Aw)]
+    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
+    Iq = np.sum(Zq*Aw)*0.25
+    return n**2, Iq
+
 def gridded_2d(q, n=300):
     """
@@ -395 +424 @@
     print("gauss-1025", *gauss_quad_2d(Q, n=1025))
     print("gauss-2049", *gauss_quad_2d(Q, n=2049))
+    print("gauss-20 usub", *gauss_quad_usub(Q, n=20))
+    print("gauss-76 usub", *gauss_quad_usub(Q, n=76))
+    print("gauss-150 usub", *gauss_quad_usub(Q, n=150))
     #gridded_2d(Q, n=2**8+1)
     gridded_2d(Q, n=2**10+1)