source: sasmodels/sasmodels/sesans.py @ d7af1c6

Last change on this file since d7af1c6 was d7af1c6, checked in by Adam Washington <adam.washington@…>, 6 years ago

First attempt at log spaced sesans

  • Property mode set to 100644
File size: 3.2 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 j0, 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        #import logging; logging.info("creating SESANS transform")
47        self.q = z
48        self._set_hankel(SElength, lam, zaccept, Rmax)
49
50    def apply(self, Iq):
51        # tye: (np.ndarray) -> np.ndarray
52        G0 = np.dot(self._H0, Iq)
53        G = np.dot(self._H.T, Iq)
54        P = G - G0
55        return P
56
57    def _set_hankel(self, SElength, lam, zaccept, Rmax):
58        # type: (np.ndarray, float, float) -> None
59        # Force float32 arrays, otherwise run into memory problems on some machines
60        SElength = np.asarray(SElength, dtype='float32')
61
62        # Rmax = #value in text box somewhere in FitPage?
63        q_max = 2*pi / (SElength[1] - SElength[0])
64        q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1])
65        # q = np.arange(q_min, q_max, q_min, dtype='float32')
66        q = np.exp(np.arange(np.log(q_min), np.log(q_max), np.log(2),
67                             dtype=np.float32))
68        q = np.hstack([0], q)
69
70        H0 = np.pi * (q[1:]**2 - q[-1]**2)
71
72        # repq = np.tile(q, (SElength.size, 1)).T
73        H = np.outer(q, SElength)
74        j1(H, out=H)
75        H *= q
76        H = H[1:] - H[:-1]
77        H *= 2 * np.pi / SElength
78
79        lam = np.asarray(lam, dtype=np.float32)
80        reptheta = np.outer(q, lam)
81        reptheta /= np.float32(2*np.pi)
82        np.arcsin(reptheta, out=reptheta)
83        # reptheta = np.arcsin(repq*replam/2*np.pi)
84        mask = reptheta > zaccept
85        H[mask] = 0
86
87        # H = np.zeros((q.size, SElength.size), dtype=np.float32)
88        # H0 = q * 0
89        assert(H.shape == (q.size, SElength.size))
90
91        self.q_calc = q
92        self._H, self._H0 = H, H0
Note: See TracBrowser for help on using the repository browser.