Changes in / [598a354:298d2d4] in sasmodels

Files:

• doc/guide/sesans/sans_to_sesans.rst (modified) (1 diff)
• sasmodels/sesans.py (modified) (3 diffs)

doc/guide/sesans/sans_to_sesans.rst

@@ -31 +31 @@
 
 in which :math:`t` is the thickness of the sample and :math:`\lambda` is the wavelength of the neutrons.
+
+Log Spaced SESANS
+-----------------
+
+For computational efficiency, the integral in the Hankel transform is
+converted into a Reimann sum
+
+
+.. math:: G(\delta) \approx
+          2 \pi
+          \sum_{Q=q_{min}}^{q_{max}} J_0(Q \delta)
+          \frac{d \Sigma}{d \Omega} (Q)
+          Q \Delta Q \!
+
+However, this model approximates more than is strictly necessary.
+Specifically, it is approximating the entire integral, when it is only
+the scattering function that cannot be handled analytically.  A better
+approximation might be
+
+.. math:: G(\delta) \approx
+          \sum_{n=0} 2 \pi \frac{d \Sigma}{d \Omega} (q_n)
+          \int_{q_{n-1}}^{q_n} J_0(Q \delta) Q dQ
+          =
+          \sum_{n=0} \frac{2 \pi}{\delta} \frac{d \Sigma}{d \Omega} (q_n)
+          (q_n J_1(q_n \delta) - q_{n-1}J_1(q_{n-1} \delta))\!,
+
+Assume that vectors :math:`q_n` and :math:`I_n` represent the q points
+and corresponding intensity data, respectively.  Further assume that
+:math:`\delta_m` and :math:`G_m` are the spin echo lengths and
+corresponding Hankel transform value.
+
+.. math:: G_m = H_{nm} I_n
+
+where
+
+.. math:: H_{nm} = \frac{2 \pi}{\delta_m}
+          (q_n J_1(q_n \delta_m) - q_{n-1} J_1(q_{n-1} \delta_m))
+
+Also not that, for the limit as :math:`\delta_m` approaches zero,
+
+.. math:: G(0)
+          =
+          \sum_{n=0} \pi \frac{d \Sigma}{d \Omega} (q_n) (q_n^2 - q_{n-1}^2)

sasmodels/sesans.py

@@ -14 +14 @@
 import numpy as np  # type: ignore
 from numpy import pi  # type: ignore
-from scipy.special import j0
+from scipy.special import j1
+
 
 class SesansTransform(object):
@@ -38 +39 @@
 
     # transform arrays
-    _H = None  # type: np.ndarray
-    _H0 = None # type: np.ndarray
+    _H = None   # type: np.ndarray
+    _H0 = None  # type: np.ndarray
 
     def __init__(self, z, SElength, lam, zaccept, Rmax):
         # type: (np.ndarray, float, float) -> None
-        #import logging; logging.info("creating SESANS transform")
         self.q = z
         self._set_hankel(SElength, lam, zaccept, Rmax)
@@ -59 +59 @@
     def _set_hankel(self, SElength, lam, zaccept, Rmax):
         # type: (np.ndarray, float, float) -> None
-        # Force float32 arrays, otherwise run into memory problems on some machines
+        # Force float32 arrays, otherwise run into memory problems on
+        # some machines
         SElength = np.asarray(SElength, dtype='float32')
 
-        #Rmax = #value in text box somewhere in FitPage?
+        # Rmax = #value in text box somewhere in FitPage?
         q_max = 2*pi / (SElength[1] - SElength[0])
         q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1])
-        q = np.arange(q_min, q_max, q_min, dtype='float32')
-        dq = q_min
+        # q = np.arange(q_min, q_max, q_min, dtype='float32')
+        # q = np.exp(np.arange(np.log(q_min), np.log(q_max), np.log(2),
+        #                      dtype=np.float32))
+        q = np.exp(np.linspace(np.log(q_min), np.log(q_max), 10*SElength.size,
+                               dtype=np.float32))
+        q = np.hstack([[0], q])
 
-        H0 = np.float32(dq/(2*pi)) * q
+        H0 = np.pi * (q[1:]**2 - q[:-1]**2)
 
-        repq = np.tile(q, (SElength.size, 1)).T
-        repSE = np.tile(SElength, (q.size, 1))
-        H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq
+        # repq = np.tile(q, (SElength.size, 1)).T
+        H = np.outer(q, SElength)
+        j1(H, out=H)
+        H *= q.reshape((-1, 1))
+        H = H[1:] - H[:-1]
+        H *= 2 * np.pi / SElength
 
-        replam = np.tile(lam, (q.size, 1))
-        reptheta = np.arcsin(repq*replam/2*np.pi)
+        lam = np.asarray(lam, dtype=np.float32)
+        reptheta = np.outer(q[1:], lam)
+        reptheta /= np.float32(2*np.pi)
+        np.arcsin(reptheta, out=reptheta)
+        # reptheta = np.arcsin(repq*replam/2*np.pi)
         mask = reptheta > zaccept
-        H[mask] = 0
+        # H[mask] = 0
 
-        self.q_calc = q
+        # H = np.zeros((q.size, SElength.size), dtype=np.float32)
+        # H0 = q * 0
+        assert(H.shape == (q.size-1, SElength.size))
+
+        self.q_calc = q[1:]
         self._H, self._H0 = H, H0