[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 |
---|
[22194b1] | 47 | mask = xs <= 200.0 # Only calculate gamma3 up to x=200 (as this is all that's plotted) |
---|
[a859f99] | 48 | # gamma3 = [trapz(gamma1[:n], xs[:n])/xs[n-1] for n in range(2, len(xs[mask]) + 1)]j |
---|
| 49 | # gamma3.insert(0, 1.0) # Gamma_3(0) is defined as 1 |
---|
| 50 | n = len(xs[mask]) |
---|
| 51 | gamma3 = cumtrapz(gamma1[:n], xs[:n])/xs[1:n] |
---|
[5e2267e] | 52 | gamma3 = np.hstack((1.0, gamma3)) # Gamma_3(0) is defined as 1 |
---|
[37d53225] | 53 | |
---|
| 54 | if self.check_if_cancelled(): return |
---|
| 55 | |
---|
| 56 | # ----- Interface Distribution function ----- |
---|
[a309667] | 57 | idf = dct(-qs**4 * (iqs-background)) |
---|
[37d53225] | 58 | |
---|
[a309667] | 59 | if self.check_if_cancelled(): return |
---|
[37d53225] | 60 | |
---|
[a309667] | 61 | # Manually calculate IDF(0.0), since scipy DCT tends to give us a |
---|
| 62 | # very large negative value. |
---|
| 63 | # IDF(x) = int_0^inf q^4 * I(q) * cos(q*x) * dq |
---|
| 64 | # => IDF(0) = int_0^inf q^4 * I(q) * dq |
---|
[9b90bf8] | 65 | idf[0] = trapz(-qs**4 * (iqs-background), qs) |
---|
[a309667] | 66 | idf /= Q # Normalise using scattering invariant |
---|
[37d53225] | 67 | |
---|
[412c509] | 68 | except Exception as e: |
---|
| 69 | import logging |
---|
| 70 | logger = logging.getLogger(__name__) |
---|
| 71 | logger.error(e) |
---|
| 72 | |
---|
[a2db1ab] | 73 | self.update(msg="Fourier transform failed.") |
---|
[412c509] | 74 | self.complete(transforms=None) |
---|
[a2db1ab] | 75 | return |
---|
| 76 | if self.isquit(): |
---|
| 77 | return |
---|
| 78 | self.update(msg="Fourier transform completed.") |
---|
| 79 | |
---|
[457f735] | 80 | transform1 = Data1D(xs, gamma1) |
---|
[22194b1] | 81 | transform3 = Data1D(xs[xs <= 200], gamma3) |
---|
[9b90bf8] | 82 | idf = Data1D(xs, idf) |
---|
[a2db1ab] | 83 | |
---|
[37d53225] | 84 | transforms = (transform1, transform3, idf) |
---|
[457f735] | 85 | |
---|
| 86 | self.complete(transforms=transforms) |
---|
[d03228e] | 87 | |
---|
| 88 | class HilbertThread(CalcThread): |
---|
| 89 | def __init__(self, raw_data, extrapolated_data, bg, updatefn=None, |
---|
| 90 | completefn=None): |
---|
| 91 | CalcThread.__init__(self, updatefn=updatefn, completefn=completefn) |
---|
| 92 | self.data = raw_data |
---|
| 93 | self.background = bg |
---|
| 94 | self.extrapolation = extrapolated_data |
---|
| 95 | |
---|
| 96 | def compute(self): |
---|
| 97 | qs = self.extrapolation.x |
---|
| 98 | iqs = self.extrapolation.y |
---|
| 99 | q = self.data.x |
---|
| 100 | background = self.background |
---|
| 101 | |
---|
| 102 | self.ready(delay=0.0) |
---|
| 103 | self.update(msg="Starting Hilbert transform.") |
---|
| 104 | self.ready(delay=0.0) |
---|
| 105 | if self.isquit(): |
---|
| 106 | return |
---|
| 107 | |
---|
| 108 | # TODO: Implement hilbert transform |
---|
| 109 | |
---|
| 110 | self.update(msg="Hilbert transform completed.") |
---|
| 111 | |
---|
[626f833] | 112 | self.complete(transforms=None) |
---|