1 | """ |
---|
2 | Conversion of scattering cross section from SANS in absolute |
---|
3 | units into SESANS using a Hankel transformation |
---|
4 | |
---|
5 | Everything is in units of metres except specified otherwise |
---|
6 | |
---|
7 | Wim Bouwman (w.g.bouwman@tudelft.nl), June 2013 |
---|
8 | """ |
---|
9 | |
---|
10 | from __future__ import division |
---|
11 | |
---|
12 | import numpy as np |
---|
13 | from numpy import pi, exp |
---|
14 | from scipy.special import jv as besselj |
---|
15 | #import direct_model.DataMixin as model |
---|
16 | |
---|
17 | def make_q(q_max, Rmax): |
---|
18 | r""" |
---|
19 | Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ |
---|
20 | to $q_max$. |
---|
21 | """ |
---|
22 | q_min = dq = 0.1 * 2*pi / Rmax |
---|
23 | return np.arange(q_min, q_max, dq) |
---|
24 | |
---|
25 | def make_all_q(data): |
---|
26 | if not data.needs_all_q: |
---|
27 | return [] |
---|
28 | elif needs_Iqxy(data): |
---|
29 | # compute qx, qy |
---|
30 | Qx, Qy = np.meshgrid(qx, qy) |
---|
31 | return [Qx, Qy] |
---|
32 | else: |
---|
33 | # else only need q |
---|
34 | return [q] |
---|
35 | |
---|
36 | def transform(data, q_calc, Iq_calc, qmono, Iq_mono): |
---|
37 | nqmono = len(qmono) |
---|
38 | if nqmono == 0: |
---|
39 | result = call_hankel(data, q_calc, Iq_calc) |
---|
40 | elif nqmono == 1: |
---|
41 | q = qmono[0] |
---|
42 | result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) |
---|
43 | else: |
---|
44 | Qx, Qy = [qmono[0], qmono[1]] |
---|
45 | Qx = np.reshape(Qx, nqx, nqy) |
---|
46 | Qy = np.reshape(Qy, nqx, nqy) |
---|
47 | Iq_mono = np.reshape(Iq_mono, nqx, nqy) |
---|
48 | qx = Qx[0, :] |
---|
49 | qy = Qy[:, 0] |
---|
50 | result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) |
---|
51 | |
---|
52 | return result |
---|
53 | |
---|
54 | def call_hankel(data, q_calc, Iq_calc): |
---|
55 | return hankel(data.x, data.lam * 1e-9, |
---|
56 | data.sample.thickness / 10, |
---|
57 | q_calc, Iq_calc) |
---|
58 | |
---|
59 | def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): |
---|
60 | return hankel(data.x, data.lam * 1e-9, |
---|
61 | data.sample.thickness / 10, |
---|
62 | q_calc, Iq_calc) |
---|
63 | |
---|
64 | def Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): |
---|
65 | return hankel(data.x, data.y, data.lam * 1e-9, |
---|
66 | data.sample.thickness / 10, |
---|
67 | q_calc, Iq_calc) |
---|
68 | |
---|
69 | def TotalScatter(model, parameters): #Work in progress!! |
---|
70 | # Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) |
---|
71 | allq = np.linspace(0,4*pi/wavelength,1000) |
---|
72 | allIq = 1 |
---|
73 | integral = allq*allIq |
---|
74 | |
---|
75 | |
---|
76 | |
---|
77 | def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still |
---|
78 | #============================================================================== |
---|
79 | # 2D Cosine Transform if "wavelength" is a vector |
---|
80 | #============================================================================== |
---|
81 | #allq is the q-space needed to create the total scattering cross-section |
---|
82 | |
---|
83 | Gprime = np.zeros_like(wavelength, 'd') |
---|
84 | s = np.zeros_like(wavelength, 'd') |
---|
85 | sd = np.zeros_like(wavelength, 'd') |
---|
86 | Gprime = np.zeros_like(wavelength, 'd') |
---|
87 | f = np.zeros_like(wavelength, 'd') |
---|
88 | for i, wavelength_i in enumerate(wavelength): |
---|
89 | z = magfield*wavelength_i |
---|
90 | allq=np.linspace() #for calculating the Q-range of the scattering power integral |
---|
91 | allIq=np.linspace() # This is the model applied to the allq q-space. Needs to refference the model somehow |
---|
92 | alldq = (allq[1]-allq[0])*1e10 |
---|
93 | sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) |
---|
94 | s[i]=1-exp(-sigma) |
---|
95 | for j, Iqy_j, qy_j in enumerate(qy): |
---|
96 | for k, Iqz_k, qz_k in enumerate(qz): |
---|
97 | Iq = np.sqrt(Iqy_j^2+Iqz_k^2) |
---|
98 | q = np.sqrt(qy_j^2 + qz_k^2) |
---|
99 | Gintegral = Iq*cos(z*Qz_k) |
---|
100 | Gprime[i] += Gintegral |
---|
101 | # sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] |
---|
102 | # s[i] += 1-exp(Totalscatter(modelname)*thickness) |
---|
103 | # For now, work with standard 2-phase scatter |
---|
104 | |
---|
105 | |
---|
106 | sd[i] += Iq |
---|
107 | f[i] = 1-s[i]+sd[i] |
---|
108 | P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] |
---|
109 | |
---|
110 | |
---|
111 | |
---|
112 | |
---|
113 | def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): |
---|
114 | #============================================================================== |
---|
115 | # HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS |
---|
116 | #============================================================================== |
---|
117 | #acceptq is the q-space needed to create limited acceptance effect |
---|
118 | SElength= wavelength*magfield |
---|
119 | G = np.zeros_like(SElength, 'd') |
---|
120 | threshold=2*pi*theta/wavelength |
---|
121 | for i, SElength_i in enumerate(SElength): |
---|
122 | allq=np.linspace() #for calculating the Q-range of the scattering power integral |
---|
123 | allIq=np.linspace() # This is the model applied to the allq q-space. Needs to refference the model somehow |
---|
124 | alldq = (allq[1]-allq[0])*1e10 |
---|
125 | sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) |
---|
126 | s[i]=1-exp(-sigma) |
---|
127 | |
---|
128 | dq = (q[1]-q[0])*1e10 |
---|
129 | a = (x<threshold) |
---|
130 | acceptq = a*q |
---|
131 | acceptIq = a*Iq |
---|
132 | |
---|
133 | G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) |
---|
134 | |
---|
135 | # G[i]=np.sum(integral) |
---|
136 | |
---|
137 | G *= dq*1e10*2*pi |
---|
138 | |
---|
139 | P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) |
---|
140 | |
---|
141 | def hankel(SElength, wavelength, thickness, q, Iq): |
---|
142 | r""" |
---|
143 | Compute the expected SESANS polarization for a given SANS pattern. |
---|
144 | |
---|
145 | Uses the hankel transform followed by the exponential. The values for *zz* |
---|
146 | (or spin echo length, or delta), wavelength and sample thickness should |
---|
147 | come from the dataset. $q$ should be chosen such that the oscillations |
---|
148 | in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). |
---|
149 | |
---|
150 | *SElength* [A] is the set of $z$ points at which to compute the |
---|
151 | Hankel transform |
---|
152 | |
---|
153 | *wavelength* [m] is the wavelength of each individual point *zz* |
---|
154 | |
---|
155 | *thickness* [cm] is the sample thickness. |
---|
156 | |
---|
157 | *q* [A$^{-1}$] is the set of $q$ points at which the model has been |
---|
158 | computed. These should be equally spaced. |
---|
159 | |
---|
160 | *I* [cm$^{-1}$] is the value of the SANS model at *q* |
---|
161 | """ |
---|
162 | G = np.zeros_like(SElength, 'd') |
---|
163 | #============================================================================== |
---|
164 | # Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS |
---|
165 | #============================================================================== |
---|
166 | for i, SElength_i in enumerate(SElength): |
---|
167 | integral = besselj(0, q*SElength_i)*Iq*q |
---|
168 | G[i] = np.sum(integral) |
---|
169 | |
---|
170 | # [m^-1] step size in q, needed for integration |
---|
171 | dq = (q[1]-q[0])*1e10 |
---|
172 | |
---|
173 | # integration step, convert q into [m**-1] and 2 pi circle integration |
---|
174 | G *= dq*1e10*2*pi |
---|
175 | |
---|
176 | P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) |
---|
177 | |
---|
178 | return P |
---|