source: sasmodels/sasmodels/sesans.py @ 0c7b8d8

Last change on this file since 0c7b8d8 was 0c7b8d8, checked in by awashington, 5 years ago

Use 64 bit floats in SesansTransform?

The new log spacing creates small matrices, so there's no longer a
reason not to use full double precision values.

  • Property mode set to 100644
File size: 3.3 KB
Line 
1"""
2Conversion of scattering cross section from SANS (I(q), or rather, ds/dO) in absolute
3units (cm-1)into SESANS correlation function G using a Hankel transformation, then converting
4the SESANS correlation function into polarisation from the SESANS experiment
5
6Everything is in units of metres except specified otherwise (NOT TRUE!!!)
7Everything is in conventional units (nm for spin echo length)
8
9Wim Bouwman (w.g.bouwman@tudelft.nl), June 2013
10"""
11
12from __future__ import division
13
14import numpy as np  # type: ignore
15from numpy import pi  # type: ignore
16from scipy.special import j1
17
18
19class SesansTransform(object):
20    """
21    Spin-Echo SANS transform calculator.  Similar to a resolution function,
22    the SesansTransform object takes I(q) for the set of *q_calc* values and
23    produces a transformed dataset
24
25    *SElength* (A) is the set of spin-echo lengths in the measured data.
26
27    *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin
28    echo encoding dimension (for ToF: Q of min(R) and max(lam)).
29
30    *Rmax* (A) is the maximum size sensitivity; larger radius requires more
31    computation time.
32    """
33    #: SElength from the data in the original data units; not used by transform
34    #: but the GUI uses it, so make sure that it is present.
35    q = None  # type: np.ndarray
36
37    #: q values to calculate when computing transform
38    q_calc = None  # type: np.ndarray
39
40    # transform arrays
41    _H = None   # type: np.ndarray
42    _H0 = None  # type: np.ndarray
43
44    def __init__(self, z, SElength, lam, zaccept, Rmax):
45        # type: (np.ndarray, float, float) -> None
46        self.q = z
47        self._set_hankel(SElength, lam, zaccept, Rmax)
48
49    def apply(self, Iq):
50        # type: (np.ndarray) -> np.ndarray
51        """
52        Apply the SESANS transform to the computed I(q).
53        """
54        G0 = np.dot(self._H0, Iq)
55        G = np.dot(self._H.T, Iq)
56        P = G - G0
57        return P
58
59    def _set_hankel(self, SElength, lam, zaccept, Rmax):
60        # type: (np.ndarray, float, float) -> None
61        SElength = np.asarray(SElength, dtype='float64')
62
63        # Rmax = #value in text box somewhere in FitPage?
64        q_max = 2*pi / (SElength[1] - SElength[0])
65        q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1])
66        # q = np.arange(q_min, q_max, q_min, dtype='float32')
67        # q = np.exp(np.arange(np.log(q_min), np.log(q_max), np.log(2),
68        #                      dtype=np.float32))
69        q = np.exp(np.linspace(np.log(q_min), np.log(q_max), 10*SElength.size,
70                               dtype=np.float64))
71        q = np.hstack([[0], q])
72
73        H0 = (q[1:]**2 - q[:-1]**2) / (2 * np.pi) / 2
74
75        # repq = np.tile(q, (SElength.size, 1)).T
76        H = np.outer(q, SElength)
77        j1(H, out=H)
78        H *= q.reshape((-1, 1))
79        H = H[1:] - H[:-1]
80        H /= 2 * np.pi * SElength
81
82        lam = np.asarray(lam, dtype=np.float64)
83        reptheta = np.outer(q[1:], lam)
84        reptheta /= np.float64(2*np.pi)
85        np.arcsin(reptheta, out=reptheta)
86        # reptheta = np.arcsin(repq*replam/2*np.pi)
87        mask = reptheta > zaccept
88        # H[mask] = 0
89
90        # H = np.zeros((q.size, SElength.size), dtype=np.float32)
91        # H0 = q * 0
92        assert(H.shape == (q.size-1, SElength.size))
93
94        self.q_calc = q[1:]
95        self._H, self._H0 = H, H0
Note: See TracBrowser for help on using the repository browser.