source: sasview/src/sas/sascalc/corfunc/transform_thread.py @ f7e6b30

ESS_GUIESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_opencl
Last change on this file since f7e6b30 was f7e6b30, checked in by awashington, 13 months ago

Remote artificial limits on corfunc range.

With the new navigation toolbars, they aren't needed.

  • Property mode set to 100644
File size: 3.7 KB
Line 
1from sas.sascalc.data_util.calcthread import CalcThread
2from sas.sascalc.dataloader.data_info import Data1D
3from scipy.fftpack import dct
4from scipy.integrate import trapz, cumtrapz
5import numpy as np
6from time import sleep
7
8class FourierThread(CalcThread):
9    def __init__(self, raw_data, extrapolated_data, bg, updatefn=None,
10        completefn=None):
11        CalcThread.__init__(self, updatefn=updatefn, completefn=completefn)
12        self.data = raw_data
13        self.background = bg
14        self.extrapolation = extrapolated_data
15
16    def check_if_cancelled(self):
17        if self.isquit():
18            self.update("Fourier transform cancelled.")
19            self.complete(transforms=None)
20            return True
21        return False
22
23    def compute(self):
24        qs = self.extrapolation.x
25        iqs = self.extrapolation.y
26        q = self.data.x
27        background = self.background
28
29        xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs)
30
31        self.ready(delay=0.0)
32        self.update(msg="Fourier transform in progress.")
33        self.ready(delay=0.0)
34
35        if self.check_if_cancelled(): return
36        try:
37            # ----- 1D Correlation Function -----
38            gamma1 = dct((iqs-background)*qs**2)
39            Q = gamma1.max()
40            gamma1 /= Q
41
42            if self.check_if_cancelled(): return
43
44            # ----- 3D Correlation Function -----
45            # gamma3(R) = 1/R int_{0}^{R} gamma1(x) dx
46            # trapz uses the trapezium rule to calculate the integral
47            # gamma3 = [trapz(gamma1[:n], xs[:n])/xs[n-1] for n in range(2, len(xs[mask]) + 1)]j
48            # gamma3.insert(0, 1.0) # Gamma_3(0) is defined as 1
49            n = len(xs)
50            gamma3 = cumtrapz(gamma1[:n], xs[:n])/xs[1:n]
51            gamma3 = np.hstack((1.0, gamma3)) # Gamma_3(0) is defined as 1
52
53            if self.check_if_cancelled(): return
54
55            # ----- Interface Distribution function -----
56            idf = dct(-qs**4 * (iqs-background))
57
58            if self.check_if_cancelled(): return
59
60            # Manually calculate IDF(0.0), since scipy DCT tends to give us a
61            # very large negative value.
62            # IDF(x) = int_0^inf q^4 * I(q) * cos(q*x) * dq
63            # => IDF(0) = int_0^inf q^4 * I(q) * dq
64            idf[0] = trapz(-qs**4 * (iqs-background), qs)
65            idf /= Q # Normalise using scattering invariant
66
67        except Exception as e:
68            import logging
69            logger = logging.getLogger(__name__)
70            logger.error(e)
71
72            self.update(msg="Fourier transform failed.")
73            self.complete(transforms=None)
74            return
75        if self.isquit():
76            return
77        self.update(msg="Fourier transform completed.")
78
79        transform1 = Data1D(xs, gamma1)
80        transform3 = Data1D(xs, gamma3)
81        idf = Data1D(xs, idf)
82
83        transforms = (transform1, transform3, idf)
84
85        self.complete(transforms=transforms)
86
87class HilbertThread(CalcThread):
88    def __init__(self, raw_data, extrapolated_data, bg, updatefn=None,
89        completefn=None):
90        CalcThread.__init__(self, updatefn=updatefn, completefn=completefn)
91        self.data = raw_data
92        self.background = bg
93        self.extrapolation = extrapolated_data
94
95    def compute(self):
96        qs = self.extrapolation.x
97        iqs = self.extrapolation.y
98        q = self.data.x
99        background = self.background
100
101        self.ready(delay=0.0)
102        self.update(msg="Starting Hilbert transform.")
103        self.ready(delay=0.0)
104        if self.isquit():
105            return
106
107        # TODO: Implement hilbert transform
108
109        self.update(msg="Hilbert transform completed.")
110
111        self.complete(transforms=None)
Note: See TracBrowser for help on using the repository browser.