[a2db1ab] | 1 | from sas.sascalc.data_util.calcthread import CalcThread |
---|
| 2 | from sas.sascalc.dataloader.data_info import Data1D |
---|
| 3 | from scipy.fftpack import dct |
---|
[5e2267e] | 4 | from scipy.integrate import trapz, cumtrapz |
---|
[a2db1ab] | 5 | import numpy as np |
---|
| 6 | from time import sleep |
---|
| 7 | |
---|
[d03228e] | 8 | class FourierThread(CalcThread): |
---|
[c728295] | 9 | def __init__(self, raw_data, extrapolated_data, bg, updatefn=None, |
---|
| 10 | completefn=None): |
---|
[a2db1ab] | 11 | CalcThread.__init__(self, updatefn=updatefn, completefn=completefn) |
---|
| 12 | self.data = raw_data |
---|
| 13 | self.background = bg |
---|
| 14 | self.extrapolation = extrapolated_data |
---|
| 15 | |
---|
[37d53225] | 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 | |
---|
[a2db1ab] | 23 | def compute(self): |
---|
| 24 | qs = self.extrapolation.x |
---|
| 25 | iqs = self.extrapolation.y |
---|
| 26 | q = self.data.x |
---|
| 27 | background = self.background |
---|
| 28 | |
---|
[412c509] | 29 | xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs) |
---|
| 30 | |
---|
[a2db1ab] | 31 | self.ready(delay=0.0) |
---|
[412c509] | 32 | self.update(msg="Fourier transform in progress.") |
---|
[a2db1ab] | 33 | self.ready(delay=0.0) |
---|
[412c509] | 34 | |
---|
[37d53225] | 35 | if self.check_if_cancelled(): return |
---|
[a2db1ab] | 36 | try: |
---|
[37d53225] | 37 | # ----- 1D Correlation Function ----- |
---|
[457f735] | 38 | gamma1 = dct((iqs-background)*qs**2) |
---|
[9b90bf8] | 39 | Q = gamma1.max() |
---|
| 40 | gamma1 /= Q |
---|
[412c509] | 41 | |
---|
[37d53225] | 42 | if self.check_if_cancelled(): return |
---|
| 43 | |
---|
| 44 | # ----- 3D Correlation Function ----- |
---|
[412c509] | 45 | # gamma3(R) = 1/R int_{0}^{R} gamma1(x) dx |
---|
[2a54ba5] | 46 | # trapz uses the trapezium rule to calculate the integral |
---|
[a859f99] | 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 |
---|
[f7e6b30] | 49 | n = len(xs) |
---|
[a859f99] | 50 | gamma3 = cumtrapz(gamma1[:n], xs[:n])/xs[1:n] |
---|
[5e2267e] | 51 | gamma3 = np.hstack((1.0, gamma3)) # Gamma_3(0) is defined as 1 |
---|
[37d53225] | 52 | |
---|
| 53 | if self.check_if_cancelled(): return |
---|
| 54 | |
---|
| 55 | # ----- Interface Distribution function ----- |
---|
[a309667] | 56 | idf = dct(-qs**4 * (iqs-background)) |
---|
[37d53225] | 57 | |
---|
[a309667] | 58 | if self.check_if_cancelled(): return |
---|
[37d53225] | 59 | |
---|
[a309667] | 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 |
---|
[9b90bf8] | 64 | idf[0] = trapz(-qs**4 * (iqs-background), qs) |
---|
[a309667] | 65 | idf /= Q # Normalise using scattering invariant |
---|
[37d53225] | 66 | |
---|
[412c509] | 67 | except Exception as e: |
---|
| 68 | import logging |
---|
| 69 | logger = logging.getLogger(__name__) |
---|
| 70 | logger.error(e) |
---|
| 71 | |
---|
[a2db1ab] | 72 | self.update(msg="Fourier transform failed.") |
---|
[412c509] | 73 | self.complete(transforms=None) |
---|
[a2db1ab] | 74 | return |
---|
| 75 | if self.isquit(): |
---|
| 76 | return |
---|
| 77 | self.update(msg="Fourier transform completed.") |
---|
| 78 | |
---|
[457f735] | 79 | transform1 = Data1D(xs, gamma1) |
---|
[f7e6b30] | 80 | transform3 = Data1D(xs, gamma3) |
---|
[9b90bf8] | 81 | idf = Data1D(xs, idf) |
---|
[a2db1ab] | 82 | |
---|
[37d53225] | 83 | transforms = (transform1, transform3, idf) |
---|
[457f735] | 84 | |
---|
| 85 | self.complete(transforms=transforms) |
---|
[d03228e] | 86 | |
---|
| 87 | class 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 | |
---|
[626f833] | 111 | self.complete(transforms=None) |
---|