[c97724e] | 1 | """ |
---|
[d459d4e] | 2 | Conversion of scattering cross section from SANS (I(q), or rather, ds/dO) in absolute |
---|
| 3 | units (cm-1)into SESANS correlation function G using a Hankel transformation, then converting |
---|
| 4 | the SESANS correlation function into polarisation from the SESANS experiment |
---|
[c97724e] | 5 | |
---|
[d459d4e] | 6 | Everything is in units of metres except specified otherwise (NOT TRUE!!!) |
---|
| 7 | Everything is in conventional units (nm for spin echo length) |
---|
[c97724e] | 8 | |
---|
| 9 | Wim Bouwman (w.g.bouwman@tudelft.nl), June 2013 |
---|
| 10 | """ |
---|
| 11 | |
---|
| 12 | from __future__ import division |
---|
| 13 | |
---|
[7ae2b7f] | 14 | import numpy as np # type: ignore |
---|
| 15 | from numpy import pi, exp # type: ignore |
---|
[c97724e] | 16 | from scipy.special import jv as besselj |
---|
[02e70ff] | 17 | #import direct_model.DataMixin as model |
---|
| 18 | |
---|
[09ebe8c] | 19 | def make_q(q_max, Rmax): |
---|
[384d114] | 20 | r""" |
---|
[09ebe8c] | 21 | Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ |
---|
[d459d4e] | 22 | to $q_max$. This is the integration range of the Hankel transform; bigger range and |
---|
| 23 | more points makes a better numerical integration. |
---|
| 24 | Smaller q_min will increase reliable spin echo length range. |
---|
| 25 | Rmax is the "radius" of the largest expected object and can be set elsewhere. |
---|
| 26 | q_max is determined by the acceptance angle of the SESANS instrument. |
---|
[09ebe8c] | 27 | """ |
---|
[a154ad16] | 28 | from sas.sascalc.data_util.nxsunit import Converter |
---|
| 29 | |
---|
[c97724e] | 30 | q_min = dq = 0.1 * 2*pi / Rmax |
---|
[2684d45] | 31 | return np.arange(q_min, |
---|
| 32 | Converter(q_max[1])(q_max[0], |
---|
| 33 | units="1/A"), |
---|
| 34 | dq) |
---|
[02e70ff] | 35 | |
---|
[a06430c] | 36 | def make_all_q(data): |
---|
[d459d4e] | 37 | """ |
---|
| 38 | Return a $q$ vector suitable for calculating the total scattering cross section for |
---|
| 39 | calculating the effect of finite acceptance angles on Time of Flight SESANS instruments. |
---|
| 40 | If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed. |
---|
| 41 | If the instrument has a rectangular acceptance, 2 all_q vectors are needed. |
---|
| 42 | If the instrument has a circular acceptance, 1 all_q vector is needed |
---|
| 43 | |
---|
| 44 | """ |
---|
| 45 | if not data.has_no_finite_acceptance: |
---|
[02e70ff] | 46 | return [] |
---|
[d459d4e] | 47 | elif data.has_yz_acceptance(data): |
---|
[02e70ff] | 48 | # compute qx, qy |
---|
| 49 | Qx, Qy = np.meshgrid(qx, qy) |
---|
| 50 | return [Qx, Qy] |
---|
| 51 | else: |
---|
| 52 | # else only need q |
---|
[d459d4e] | 53 | # data.has_z_acceptance |
---|
[02e70ff] | 54 | return [q] |
---|
| 55 | |
---|
| 56 | def transform(data, q_calc, Iq_calc, qmono, Iq_mono): |
---|
[d459d4e] | 57 | """ |
---|
| 58 | Decides which transform type is to be used, based on the experiment data file contents (header) |
---|
| 59 | (2016-03-19: currently controlled from parameters script) |
---|
| 60 | nqmono is the number of q vectors to be used for the detector integration |
---|
| 61 | """ |
---|
[02e70ff] | 62 | nqmono = len(qmono) |
---|
| 63 | if nqmono == 0: |
---|
| 64 | result = call_hankel(data, q_calc, Iq_calc) |
---|
| 65 | elif nqmono == 1: |
---|
| 66 | q = qmono[0] |
---|
| 67 | result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) |
---|
| 68 | else: |
---|
| 69 | Qx, Qy = [qmono[0], qmono[1]] |
---|
| 70 | Qx = np.reshape(Qx, nqx, nqy) |
---|
| 71 | Qy = np.reshape(Qy, nqx, nqy) |
---|
| 72 | Iq_mono = np.reshape(Iq_mono, nqx, nqy) |
---|
| 73 | qx = Qx[0, :] |
---|
| 74 | qy = Qy[:, 0] |
---|
| 75 | result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) |
---|
| 76 | |
---|
| 77 | return result |
---|
| 78 | |
---|
| 79 | def call_hankel(data, q_calc, Iq_calc): |
---|
[2684d45] | 80 | return hankel((data.x, data.x_unit), |
---|
| 81 | (data.lam, data.lam_unit), |
---|
| 82 | (data.sample.thickness, |
---|
| 83 | data.sample.thickness_unit), |
---|
[02e70ff] | 84 | q_calc, Iq_calc) |
---|
| 85 | |
---|
| 86 | def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): |
---|
| 87 | return hankel(data.x, data.lam * 1e-9, |
---|
| 88 | data.sample.thickness / 10, |
---|
| 89 | q_calc, Iq_calc) |
---|
| 90 | |
---|
[d459d4e] | 91 | def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): |
---|
[251b40b] | 92 | return hankel(data.x, data.y, data.lam * 1e-9, |
---|
[02e70ff] | 93 | data.sample.thickness / 10, |
---|
| 94 | q_calc, Iq_calc) |
---|
| 95 | |
---|
| 96 | def TotalScatter(model, parameters): #Work in progress!! |
---|
| 97 | # 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) |
---|
| 98 | allq = np.linspace(0,4*pi/wavelength,1000) |
---|
[251b40b] | 99 | allIq = 1 |
---|
[02e70ff] | 100 | integral = allq*allIq |
---|
| 101 | |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still |
---|
| 105 | #============================================================================== |
---|
| 106 | # 2D Cosine Transform if "wavelength" is a vector |
---|
| 107 | #============================================================================== |
---|
| 108 | #allq is the q-space needed to create the total scattering cross-section |
---|
| 109 | |
---|
| 110 | Gprime = np.zeros_like(wavelength, 'd') |
---|
| 111 | s = np.zeros_like(wavelength, 'd') |
---|
| 112 | sd = np.zeros_like(wavelength, 'd') |
---|
| 113 | Gprime = np.zeros_like(wavelength, 'd') |
---|
| 114 | f = np.zeros_like(wavelength, 'd') |
---|
[251b40b] | 115 | for i, wavelength_i in enumerate(wavelength): |
---|
| 116 | z = magfield*wavelength_i |
---|
| 117 | allq=np.linspace() #for calculating the Q-range of the scattering power integral |
---|
| 118 | allIq=np.linspace() # This is the model applied to the allq q-space. Needs to refference the model somehow |
---|
| 119 | alldq = (allq[1]-allq[0])*1e10 |
---|
| 120 | sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) |
---|
| 121 | s[i]=1-exp(-sigma) |
---|
| 122 | for j, Iqy_j, qy_j in enumerate(qy): |
---|
| 123 | for k, Iqz_k, qz_k in enumerate(qz): |
---|
| 124 | Iq = np.sqrt(Iqy_j^2+Iqz_k^2) |
---|
| 125 | q = np.sqrt(qy_j^2 + qz_k^2) |
---|
| 126 | Gintegral = Iq*cos(z*Qz_k) |
---|
| 127 | Gprime[i] += Gintegral |
---|
| 128 | # sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] |
---|
| 129 | # s[i] += 1-exp(Totalscatter(modelname)*thickness) |
---|
| 130 | # For now, work with standard 2-phase scatter |
---|
| 131 | |
---|
| 132 | |
---|
| 133 | sd[i] += Iq |
---|
| 134 | f[i] = 1-s[i]+sd[i] |
---|
| 135 | P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] |
---|
[02e70ff] | 136 | |
---|
| 137 | |
---|
| 138 | |
---|
| 139 | |
---|
| 140 | def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): |
---|
| 141 | #============================================================================== |
---|
| 142 | # HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS |
---|
| 143 | #============================================================================== |
---|
| 144 | #acceptq is the q-space needed to create limited acceptance effect |
---|
| 145 | SElength= wavelength*magfield |
---|
| 146 | G = np.zeros_like(SElength, 'd') |
---|
| 147 | threshold=2*pi*theta/wavelength |
---|
[251b40b] | 148 | for i, SElength_i in enumerate(SElength): |
---|
| 149 | allq=np.linspace() #for calculating the Q-range of the scattering power integral |
---|
| 150 | allIq=np.linspace() # This is the model applied to the allq q-space. Needs to refference the model somehow |
---|
| 151 | alldq = (allq[1]-allq[0])*1e10 |
---|
| 152 | sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) |
---|
| 153 | s[i]=1-exp(-sigma) |
---|
| 154 | |
---|
| 155 | dq = (q[1]-q[0])*1e10 |
---|
| 156 | a = (x<threshold) |
---|
| 157 | acceptq = a*q |
---|
| 158 | acceptIq = a*Iq |
---|
| 159 | |
---|
| 160 | G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) |
---|
| 161 | |
---|
| 162 | # G[i]=np.sum(integral) |
---|
| 163 | |
---|
| 164 | G *= dq*1e10*2*pi |
---|
| 165 | |
---|
| 166 | P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) |
---|
[02e70ff] | 167 | |
---|
[c97724e] | 168 | def hankel(SElength, wavelength, thickness, q, Iq): |
---|
[384d114] | 169 | r""" |
---|
[c97724e] | 170 | Compute the expected SESANS polarization for a given SANS pattern. |
---|
| 171 | |
---|
[384d114] | 172 | Uses the hankel transform followed by the exponential. The values for *zz* |
---|
| 173 | (or spin echo length, or delta), wavelength and sample thickness should |
---|
| 174 | come from the dataset. $q$ should be chosen such that the oscillations |
---|
| 175 | in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). |
---|
[c97724e] | 176 | |
---|
[384d114] | 177 | *SElength* [A] is the set of $z$ points at which to compute the |
---|
| 178 | Hankel transform |
---|
[c97724e] | 179 | |
---|
| 180 | *wavelength* [m] is the wavelength of each individual point *zz* |
---|
| 181 | |
---|
| 182 | *thickness* [cm] is the sample thickness. |
---|
| 183 | |
---|
[384d114] | 184 | *q* [A$^{-1}$] is the set of $q$ points at which the model has been |
---|
| 185 | computed. These should be equally spaced. |
---|
[c97724e] | 186 | |
---|
[384d114] | 187 | *I* [cm$^{-1}$] is the value of the SANS model at *q* |
---|
[c97724e] | 188 | """ |
---|
[2684d45] | 189 | |
---|
| 190 | from sas.sascalc.data_util.nxsunit import Converter |
---|
| 191 | wavelength = Converter(wavelength[1])(wavelength[0],"A") |
---|
| 192 | thickness = Converter(thickness[1])(thickness[0],"A") |
---|
| 193 | Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters |
---|
| 194 | SElength = Converter(SElength[1])(SElength[0],"A") |
---|
| 195 | |
---|
[190fc2b] | 196 | G = np.zeros_like(SElength, 'd') |
---|
[02e70ff] | 197 | #============================================================================== |
---|
| 198 | # Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS |
---|
| 199 | #============================================================================== |
---|
[190fc2b] | 200 | for i, SElength_i in enumerate(SElength): |
---|
| 201 | integral = besselj(0, q*SElength_i)*Iq*q |
---|
| 202 | G[i] = np.sum(integral) |
---|
[c94577f] | 203 | G0 = np.sum(Iq*q) |
---|
[384d114] | 204 | |
---|
| 205 | # [m^-1] step size in q, needed for integration |
---|
[2684d45] | 206 | dq = (q[1]-q[0]) |
---|
[384d114] | 207 | |
---|
| 208 | # integration step, convert q into [m**-1] and 2 pi circle integration |
---|
[2684d45] | 209 | G *= dq*2*pi |
---|
| 210 | G0 = np.sum(Iq*q)*dq*2*np.pi |
---|
[384d114] | 211 | |
---|
[c94577f] | 212 | P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0)) |
---|
[c97724e] | 213 | |
---|
| 214 | return P |
---|