Changeset 49d1d42f in sasmodels

Timestamp:

Mar 11, 2015 8:58:05 PM (10 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:

aa4946b

Parents:

63b32bb

Message:

working demo of 1D resolution calculator (still incomplete and incorrect)

File:

• sasmodels/resolution.py (modified) (6 diffs)

sasmodels/resolution.py

@@ -4 +4 @@
 
 SLIT_SMEAR_POINTS = 500
+
+class MatrixResolution:
+    def apply(self, Iq):
+        return np.dot(Iq, self.resolution_matrix)
+
+class Pinhole1D(MatrixResolution):
+    def __init__(self, q, q_width):
+        self.q, self.q_width = q, q_width
+        self.q_calc = pinhole_extend_q(q,q_width)
+        self.resolution_matrix = \
+            pinhole_resolution(self.q_calc, self.q, self.q_width)
+
+class Slit1D(MatrixResolution):
+    def __init__(self, q, qx_width, qy_width):
+        if np.isscalar(qx_width):
+            qx_width = qx_width*np.ones(len(q))
+        if np.isscalar(qy_width):
+            qy_width = qy_width*np.ones(len(q))
+        self.q, self.qx_width, self.qy_width = q, qx_width, qy_width
+        self.q_calc = slit_extend_q(q, qx_width, qy_width)
+        self.resolution_matrix = \
+            slit_resolution(self.q_calc, self.q, self.qx_width, self.qy_width)
 
 def pinhole_resolution(q_calc, q, q_width):
@@ -19 +41 @@
     G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:] )
     weights = G[1:] - G[:-1]
-    weights /= sum(weights, axis=1)
+    weights /= np.sum(weights, axis=1)
     return weights
+
+
+def pinhole_extend_q(q, q_width):
+    return q
+
 
 def slit_resolution(q_calc, q, qx_width, qy_width):
@@ -27 +54 @@
     qy_min, qy_max = 0.0, edges[-1]
 
+    # Make q_calc into a row vector, and q, qx_width, qy_width into columns
+    # Make weights match [ q_calc X q ]
     weights = np.zeros((len(q),len(q_calc)),'d')
+    q_calc = q_calc[None,:]
+    q, qx_width, qy_width, edges = [
+        v[:,None] for v in (q, qx_width, qy_width, edges)]
+
     # Loop for width (height is analytical).
     # Condition: height >>> width, otherwise, below is not accurate enough.
@@ -37 +70 @@
     # In general, this weight numerical iteration for width>0 might be a rough
     # approximation, but it must be good enough when height >>> width.
-    E_sq = edges**2[:,None]
+    E_sq = edges**2
     y_points = SLIT_SMEAR_POINTS if np.any(qy_width>0) else 1
     qy_step = 0 if y_points == 1 else qy_width/(y_points-1)
@@ -44 +77 @@
         qx_low = qy
         qx_high = sqrt(qx_low**2 + qx_width**2)
-        in_x = (q_calc[:,None]>=qx_low[None,:])*(q_calc[:,None]<=qx_high[None,:])
-        qy_sq = qy**2[None,:]
-        weights += (sqrt(E_sq[1:]-qy_sq) - sqrt(E_sq[:-1]-qy_sq))*in_x
-    weights /= sum(weights, axis=1)
+        in_x = (q_calc>=qx_low)*(q_calc<=qx_high)
+        qy_sq = qy**2
+        weights += (sqrt(E_sq[1:]-qy_sq) - sqrt(qy_sq - E_sq[:-1]))*in_x
+    weights /= np.sum(weights, axis=1)
     return weights
+
+
+def slit_extend_q(q, qx_width, qy_width):
+    return q
+
 
 def bin_edges(x):
@@ -59 +97 @@
         ])
     return edges
+
+
+############################################################################
+# usage demo
+############################################################################
+
+def _eval_demo_1d(resolution, title):
+    from sasmodels import core
+    from sasmodels.models import cylinder
+    ## or alternatively:
+    # cylinder = core.load_model_definition('cylinder')
+    model = core.load_model(cylinder)
+
+    kernel = core.make_kernel(model, [resolution.q_calc])
+    Iq_calc = core.call_kernel(kernel, {'length':210, 'radius':500})
+    Iq = resolution.apply(Iq_calc)
+
+    import matplotlib.pyplot as plt
+    plt.loglog(resolution.q_calc, Iq_calc, label='unsmeared')
+    plt.loglog(resolution.q, Iq, label='smeared', hold=True)
+    plt.legend()
+    plt.title(title)
+    plt.xlabel("Q (1/Ang)")
+    plt.ylabel("I(Q) (1/cm)")
+
+def demo_pinhole_1d():
+    q = np.logspace(-3,-1,400)
+    dq = 0.1*q
+    resolution = Pinhole1D(q, dq)
+    _eval_demo_1d(resolution, title="10% dQ/Q Pinhole Resolution")
+
+def demo_slit_1d():
+    q = np.logspace(-3,-1,400)
+    qx_width = 0.005
+    qy_width = 0.0
+    resolution = Slit1D(q, qx_width, qy_width)
+    _eval_demo_1d(resolution, title="0.005 Qx Slit Resolution")
+
+def demo():
+    import matplotlib.pyplot as plt
+    plt.subplot(121)
+    demo_pinhole_1d()
+    plt.subplot(122)
+    demo_slit_1d()
+    plt.show()
+
+
+if __name__ == "__main__":
+    demo()
+
+