1 | from sas.sascalc.data_util.calcthread import CalcThread |
---|
2 | from sas.sascalc.dataloader.data_info import Data1D |
---|
3 | from scipy.fftpack import dct |
---|
4 | from scipy.integrate import trapz |
---|
5 | import numpy as np |
---|
6 | from time import sleep |
---|
7 | |
---|
8 | class FourierThread(CalcThread): |
---|
9 | def __init__(self, raw_data, extrapolated_data, bg, extrap_fn=None, |
---|
10 | updatefn=None, 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 | self.extrap_fn = extrap_fn |
---|
16 | |
---|
17 | def check_if_cancelled(self): |
---|
18 | if self.isquit(): |
---|
19 | self.update("Fourier transform cancelled.") |
---|
20 | self.complete(transforms=None) |
---|
21 | return True |
---|
22 | return False |
---|
23 | |
---|
24 | def compute(self): |
---|
25 | qs = self.extrapolation.x |
---|
26 | iqs = self.extrapolation.y |
---|
27 | q = self.data.x |
---|
28 | background = self.background |
---|
29 | |
---|
30 | xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs) |
---|
31 | |
---|
32 | self.ready(delay=0.0) |
---|
33 | self.update(msg="Fourier transform in progress.") |
---|
34 | self.ready(delay=0.0) |
---|
35 | |
---|
36 | if self.check_if_cancelled(): return |
---|
37 | try: |
---|
38 | # ----- 1D Correlation Function ----- |
---|
39 | gamma1 = dct((iqs-background)*qs**2) |
---|
40 | gamma1 = gamma1 / gamma1.max() |
---|
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 | mask = xs <= 200.0 # Only calculate gamma3 up to x=200 (as this is all that's plotted) |
---|
48 | gamma3 = [trapz(gamma1[:n], xs[:n])/xs[n-1] for n in range(2, len(xs[mask]) + 1)] |
---|
49 | gamma3.insert(0, 1.0) # Gamma_3(0) is defined as 1 |
---|
50 | gamma3 = np.array(gamma3) |
---|
51 | |
---|
52 | if self.check_if_cancelled(): return |
---|
53 | |
---|
54 | # ----- Interface Distribution function ----- |
---|
55 | dmax = 200.0 # Max real space value to calculate IDF up to |
---|
56 | dstep = 0.5 |
---|
57 | qmax = 1.0 # Max q value to integrate up to when calculating IDF |
---|
58 | |
---|
59 | # Units of x axis depend on qmax (for some reason?). This scales |
---|
60 | # the xgamma array appropriately, since qmax was set to 0.6 in |
---|
61 | # the original fortran code. |
---|
62 | x_scale = qmax / 0.6 |
---|
63 | |
---|
64 | xgamma = np.arange(0, dmax/x_scale, step=dstep/x_scale) |
---|
65 | idf = np.zeros(len(xgamma)) |
---|
66 | |
---|
67 | # nth moment = integral(q^n * I(q), q=0, q=inf) |
---|
68 | moment = np.zeros(5) |
---|
69 | for n in range(5): |
---|
70 | integrand = qs**n * (iqs-background) |
---|
71 | moment[n] = trapz(integrand[qs < qmax], qs[qs < qmax]) |
---|
72 | if self.check_if_cancelled(): return |
---|
73 | |
---|
74 | # idf(x) = -integral(q^4 * I(q)*cos(qx), q=0, q=inf) / 2nd moment |
---|
75 | # => idf(0) = -integral(q^4 * I(q), 0, inf) / (2nd moment) |
---|
76 | # = -(4th moment)/(2nd moment) |
---|
77 | idf[0] = -moment[4] / moment[2] |
---|
78 | for i in range(1, len(xgamma)): |
---|
79 | d = xgamma[i] |
---|
80 | |
---|
81 | integrand = -qs**4 * (iqs-background) * np.cos(d*qs) |
---|
82 | idf[i] = trapz(integrand[qs < qmax], qs[qs < qmax]) |
---|
83 | idf[i] /= moment[2] |
---|
84 | if self.check_if_cancelled(): return |
---|
85 | |
---|
86 | xgamma *= x_scale |
---|
87 | |
---|
88 | except Exception as e: |
---|
89 | import logging |
---|
90 | logger = logging.getLogger(__name__) |
---|
91 | logger.error(e) |
---|
92 | |
---|
93 | self.update(msg="Fourier transform failed.") |
---|
94 | self.complete(transforms=None) |
---|
95 | return |
---|
96 | if self.isquit(): |
---|
97 | return |
---|
98 | self.update(msg="Fourier transform completed.") |
---|
99 | |
---|
100 | transform1 = Data1D(xs, gamma1) |
---|
101 | transform3 = Data1D(xs[xs <= 200], gamma3) |
---|
102 | idf = Data1D(xgamma, idf) |
---|
103 | |
---|
104 | transforms = (transform1, transform3, idf) |
---|
105 | |
---|
106 | self.complete(transforms=transforms) |
---|
107 | |
---|
108 | class HilbertThread(CalcThread): |
---|
109 | def __init__(self, raw_data, extrapolated_data, bg, updatefn=None, |
---|
110 | completefn=None): |
---|
111 | CalcThread.__init__(self, updatefn=updatefn, completefn=completefn) |
---|
112 | self.data = raw_data |
---|
113 | self.background = bg |
---|
114 | self.extrapolation = extrapolated_data |
---|
115 | |
---|
116 | def compute(self): |
---|
117 | qs = self.extrapolation.x |
---|
118 | iqs = self.extrapolation.y |
---|
119 | q = self.data.x |
---|
120 | background = self.background |
---|
121 | |
---|
122 | self.ready(delay=0.0) |
---|
123 | self.update(msg="Starting Hilbert transform.") |
---|
124 | self.ready(delay=0.0) |
---|
125 | if self.isquit(): |
---|
126 | return |
---|
127 | |
---|
128 | # TODO: Implement hilbert transform |
---|
129 | |
---|
130 | self.update(msg="Hilbert transform completed.") |
---|
131 | |
---|
132 | self.complete(transforms=None) |
---|