Changeset c9d052d in sasmodels

Timestamp:

Jun 17, 2018 1:45:45 PM (7 years ago)

Author:

awashington

Branches:

master

Children:

298d2d4

Parents:

befe905 (diff), b032200 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'log_sesans' of github.com:rprospero/sasmodels into log_sesans

Files:

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

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)
@@ -56 +56 @@
     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

doc/guide/magnetism/magnetism.rst

@@ -39 +39 @@
 
 .. math::
-    -- &= ((1-u_i)(1-u_f))^{1/4} \\
-    -+ &= ((1-u_i)(u_f))^{1/4} \\
-    +- &= ((u_i)(1-u_f))^{1/4} \\
-    ++ &= ((u_i)(u_f))^{1/4}
+    -- &= (1-u_i)(1-u_f) \\
+    -+ &= (1-u_i)(u_f) \\
+    +- &= (u_i)(1-u_f) \\
+    ++ &= (u_i)(u_f)
 
 Ideally the experiment would measure the pure spin states independently and
@@ -104 +104 @@
 | 2015-05-02 Steve King
 | 2017-11-15 Paul Kienzle
+| 2018-06-02 Adam Washington

sasmodels/kernel_iq.c

@@ -83 +83 @@
   in_spin = clip(in_spin, 0.0, 1.0);
   out_spin = clip(out_spin, 0.0, 1.0);
-  // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that
+  // Previous version of this function took the square root of the weights,
+  // under the assumption that 
+  //
   //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...)
-  // which is likely to be the case for simple models.
-  weight[0] = sqrt((1.0-in_spin) * (1.0-out_spin)); // dd
-  weight[1] = sqrt((1.0-in_spin) * out_spin);       // du.real
-  weight[2] = sqrt(in_spin * (1.0-out_spin));       // ud.real
-  weight[3] = sqrt(in_spin * out_spin);             // uu
+  //
+  // However, since the weights are applied to the final intensity and
+  // are not interned inside the I(q) function, we want the full
+  // weight and not the square root.  Any function using
+  // set_spin_weights as part of calculating an amplitude will need to
+  // manually take that square root, but there is currently no such
+  // function.
+  weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd
+  weight[1] = (1.0-in_spin) * out_spin;       // du
+  weight[2] = in_spin * (1.0-out_spin);       // ud
+  weight[3] = in_spin * out_spin;             // uu
   weight[4] = weight[1]; // du.imag
   weight[5] = weight[2]; // ud.imag