Changeset 237b7392 in sasmodels
 Timestamp:
 Feb 28, 2017 1:20:33 PM (2 years ago)
 Parents:
 23054a1 (diff), 7fa0e08 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.  gitauthor:
 Jeff Krzywon <krzywon@…> (02/28/17 13:20:33)
 gitcommitter:
 GitHub <noreply@…> (02/28/17 13:20:33)
 Files:

 1 deleted
 1 edited
Legend:
 Unmodified
 Added
 Removed

sasmodels/sesans.py
r94d13f1 r7fa0e08 14 14 import numpy as np # type: ignore 15 15 from numpy import pi, exp # type: ignore 16 from scipy.special import j0 17 18 class SesansTransform(object): 19 """ 20 SpinEcho SANS transform calculator. Similar to a resolution function, 21 the SesansTransform object takes I(q) for the set of *q_calc* values and 22 produces a transformed dataset 23 24 *SElength* (A) is the set of spinecho lengths in the measured data. 25 26 *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin 27 echo encoding dimension (for ToF: Q of min(R) and max(lam)). 28 29 *Rmax* (A) is the maximum size sensitivity; larger radius requires more 30 computation time. 31 """ 32 #: SElength from the data in the original data units; not used by transform 33 #: but the GUI uses it, so make sure that it is present. 34 q = None # type: np.ndarray 35 36 #: q values to calculate when computing transform 37 q_calc = None # type: np.ndarray 38 39 # transform arrays 40 _H = None # type: np.ndarray 41 _H0 = None # type: np.ndarray 42 43 def __init__(self, z, SElength, zaccept, Rmax): 44 # type: (np.ndarray, float, float) > None 45 #import logging; logging.info("creating SESANS transform") 46 self.q = z 47 self._set_hankel(SElength, zaccept, Rmax) 48 49 def apply(self, Iq): 50 # tye: (np.ndarray) > np.ndarray 51 G0 = np.dot(self._H0, Iq) 52 G = np.dot(self._H.T, Iq) 53 P = G  G0 54 return P 55 56 def _set_hankel(self, SElength, zaccept, Rmax): 57 # type: (np.ndarray, float, float) > None 58 # Force float32 arrays, otherwise run into memory problems on some machines 59 SElength = np.asarray(SElength, dtype='float32') 60 61 #Rmax = #value in text box somewhere in FitPage? 62 q_max = 2*pi / (SElength[1]  SElength[0]) 63 q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[1]) 64 q = np.arange(q_min, q_max, q_min, dtype='float32') 65 dq = q_min 66 67 H0 = np.float32(dq/(2*pi)) * q 68 69 repq = np.tile(q, (SElength.size, 1)).T 70 repSE = np.tile(SElength, (q.size, 1)) 71 H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq 72 73 self.q_calc = q 74 self._H, self._H0 = H, H0 16 from scipy.special import jv as besselj 17 #import direct_model.DataMixin as model 18 19 def make_q(q_max, Rmax): 20 r""" 21 Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ 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. 27 """ 28 from sas.sascalc.data_util.nxsunit import Converter 29 30 q_min = dq = 0.1 * 2*pi / Rmax 31 return np.arange(q_min, 32 Converter(q_max[1])(q_max[0], 33 units="1/A"), 34 dq) 35 36 def make_all_q(data): 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: 46 return [] 47 elif data.has_yz_acceptance(data): 48 # compute qx, qy 49 Qx, Qy = np.meshgrid(qx, qy) 50 return [Qx, Qy] 51 else: 52 # else only need q 53 # data.has_z_acceptance 54 return [q] 55 56 def transform(data, q_calc, Iq_calc, qmono, Iq_mono): 57 """ 58 Decides which transform type is to be used, based on the experiment data file contents (header) 59 (20160319: currently controlled from parameters script) 60 nqmono is the number of q vectors to be used for the detector integration 61 """ 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): 80 return hankel((data.x, data.x_unit), 81 (data.lam, data.lam_unit), 82 (data.sample.thickness, 83 data.sample.thickness_unit), 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 * 1e9, 88 data.sample.thickness / 10, 89 q_calc, Iq_calc) 90 91 def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): 92 return hankel(data.x, data.y, data.lam * 1e9, 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) 99 allIq = 1 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 qspace needed to create the total scattering crosssection 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') 115 for i, wavelength_i in enumerate(wavelength): 116 z = magfield*wavelength_i 117 allq=np.linspace() #for calculating the Qrange of the scattering power integral 118 allIq=np.linspace() # This is the model applied to the allq qspace. 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]=1exp(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] += 1exp(Totalscatter(modelname)*thickness) 130 # For now, work with standard 2phase scatter 131 132 133 sd[i] += Iq 134 f[i] = 1s[i]+sd[i] 135 P[i] = (1sd[i]/f[i])+1/f[i]*Gprime[i] 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 qspace needed to create limited acceptance effect 145 SElength= wavelength*magfield 146 G = np.zeros_like(SElength, 'd') 147 threshold=2*pi*theta/wavelength 148 for i, SElength_i in enumerate(SElength): 149 allq=np.linspace() #for calculating the Qrange of the scattering power integral 150 allIq=np.linspace() # This is the model applied to the allq qspace. 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]=1exp(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)*(GG[0])) 167 168 def hankel(SElength, wavelength, thickness, q, Iq): 169 r""" 170 Compute the expected SESANS polarization for a given SANS pattern. 171 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}$). 176 177 *SElength* [A] is the set of $z$ points at which to compute the 178 Hankel transform 179 180 *wavelength* [m] is the wavelength of each individual point *zz* 181 182 *thickness* [cm] is the sample thickness. 183 184 *q* [A$^{1}$] is the set of $q$ points at which the model has been 185 computed. These should be equally spaced. 186 187 *I* [cm$^{1}$] is the value of the SANS model at *q* 188 """ 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 196 G = np.zeros_like(SElength, 'd') 197 #============================================================================== 198 # Hankel Transform method if "wavelength" is a scalar; monochromatic SESANS 199 #============================================================================== 200 for i, SElength_i in enumerate(SElength): 201 integral = besselj(0, q*SElength_i)*Iq*q 202 G[i] = np.sum(integral) 203 G0 = np.sum(Iq*q) 204 205 # [m^1] step size in q, needed for integration 206 dq = (q[1]q[0]) 207 208 # integration step, convert q into [m**1] and 2 pi circle integration 209 G *= dq*2*pi 210 G0 = np.sum(Iq*q)*dq*2*np.pi 211 212 P = exp(thickness*wavelength**2/(4*pi**2)*(GG0)) 213 214 return P
Note: See TracChangeset
for help on using the changeset viewer.