Changeset 80ba1a2 in sasview
 Timestamp:
 Nov 12, 2015 6:11:10 PM (9 years ago)
 Branches:
 master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release4.1.1, release4.1.2, release4.2.2, release_4.0.1, ticket1009, ticket1094headless, ticket12422dresolution, ticket1243, ticket1249, ticket885, unittestsaveload
 Children:
 d430ee8, 662d8d87
 Parents:
 bc873053
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

src/sas/models/resolution.py
rbc873053 r80ba1a2 10 10 11 11 MINIMUM_RESOLUTION = 1e8 12 13 14 # When extrapolating to q, what is the minimum positive q relative to q_min 15 # that we wish to calculate? 16 MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION = 0.01 12 17 13 18 class Resolution(object): … … 148 153 149 154 150 def slit_resolution(q_calc, q, width, height ):155 def slit_resolution(q_calc, q, width, height, n_height=30): 151 156 r""" 152 157 Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given 153 $q_v$ = *width* and $q_h$ = *height*. 154 155 *width* and *height* are scalars since current instruments use the 156 same slit settings for all measurement points. 158 $q_\perp$ = *width* and $q_\parallel$ = *height*. *n_height* is 159 is the number of steps to use in the integration over $q_\parallel$ 160 when both $q_\perp$ and $q_\parallel$ are nonzero. 161 162 Each $q$ can have an independent width and height value even though 163 current instruments use the same slit setting for all measured points. 157 164 158 165 If slit height is large relative to width, use: … … 160 167 .. math:: 161 168 162 I_s(q_ o) = \frac{1}{\Delta q_v}163 \int_0^{\Delta q_ v} I(\sqrt{q_o^2 + u^2} du169 I_s(q_i) = \frac{1}{\Delta q_\perp} 170 \int_0^{\Delta q_\perp} I(\sqrt{q_i^2 + q_\perp^2} dq_\perp 164 171 165 172 If slit width is large relative to height, use: … … 167 174 .. math:: 168 175 169 I_s(q_o) = \frac{1}{2 \Delta q_v} 170 \int_{\Delta q_v}^{\Delta q_v} I(u) du 176 I_s(q_i) = \frac{1}{2 \Delta q_\parallel} 177 \int_{\Delta q_\parallel}^{\Delta q_\parallel} 178 I(q_i + q_\parallel) dq_\parallel 179 180 For a mixture of slit width and height use: 181 182 .. math:: 183 184 I_s(q_i) = \frac{1}{2 \Delta q_\parallel \Delta q_\perp} 185 \int_{\Delta q_\parallel)^{\Delta q_parallel} 186 \int_0^[\Delta q_\perp} 187 I(\sqrt{(q_i + q_\parallel)^2 + q_\perp^2}) 188 dq_\perp dq_\parallel 189 190 191 Algorithm 192  193 194 We are using the midpoint integration rule to assign weights to each 195 element of a weight matrix $W$ so that 196 197 .. math:: 198 199 I_s(q) = W I(q_\text{calc}) 200 201 If *q_calc* is at the midpoint, we can infer the bin edges from the 202 pairwise averages of *q_calc*, adding the missing edges before 203 *q_calc[0]* and after *q_calc[1]*. 204 205 For $q_\parallel = 0$, the smeared value can be computed numerically 206 using the $u$ substitution 207 208 .. math:: 209 210 u_j = \sqrt{q_j^2  q^2} 211 212 This gives 213 214 .. math:: 215 216 I_s(q) \approx \sum_j I(u_j) \Delta u_j 217 218 where $I(u_j)$ is the value at the midpoint, and $\Delta u_j$ is the 219 difference between consecutive edges which have been first converted 220 to $u$. Only $u_j \in [0, \Delta q_\perp]$ are used, which corresponds 221 to $q_j \in [q, \sqrt{q^2 + \Delta q_\perp}]$, so 222 223 .. math:: 224 225 W_{ij} = \frac{1}{\Delta q_\perp} \Delta u_j 226 = \frac{1}{\Delta q_\perp} 227 \sqrt{q_{j+1}^2  q_i^2}  \sqrt{q_j^2  q_i^2} 228 \text{if} q_j \in [q_i, \sqrt{q_i^2 + q_\perp^2}] 229 230 where $I_s(q_i)$ is the theory function being computed and $q_j$ are the 231 midpoints between the calculated values in *q_calc*. We tweak the 232 edges of the initial and final intervals so that they lie on integration 233 limits. 234 235 (To be precise, the transformed midpoint $u(q_j)$ is not necessarily the 236 midpoint of the edges $u((q_{j1}+q_j)/2)$ and $u((q_j + q_{j+1})/2)$, 237 but it is at least in the interval, so the approximation is going to be 238 a little better than the left or right Riemann sum, and should be 239 good enough for our purposes.) 240 241 For $q_\perp = 0$, the $u$ substitution is simpler: 242 243 .. math:: 244 245 u_j = q_j  q 246 247 so 248 249 .. math:: 250 251 W_ij = \frac{1}{2 \Delta q_\parallel} \Delta u_j 252 = \frac{1}{2 \Delta q_\parallel} (q_{j+1}  q_j) 253 \text{if} q_j \in [q\Delta q_\parallel, q+\Delta q_\parallel] 254 255 However, we need to support cases were $u_j < 0$, which means using 256 $2 (q_{j+1}  q_j)$ when $q_j \in [0, q_\parallelq_i]$. This is not 257 an issue for $q_i > q_\parallel$. 258 259 For bot $q_\perp > 0$ and $q_\parallel > 0$ we perform a 2 dimensional 260 integration with 261 262 .. math:: 263 264 u_jk = \sqrt{q_j^2  (q + (k\Delta q_\parallel/L))^2} 265 \text{for} k = L \ldots L 266 267 for $L$ = *n_height*. This gives 268 269 .. math:: 270 271 W_{ij} = \frac{1}{2 \Delta q_\perp q_\parallel} 272 \sum_{k=L}^L \Delta u_jk (\frac{\Delta q_\parallel}{2 L + 1} 273 274 171 275 """ 172 276 #np.set_printoptions(precision=6, linewidth=10000) … … 177 281 weights = np.zeros((len(q), len(q_calc)), 'd') 178 282 283 #print q_calc 179 284 for i, (qi, w, h) in enumerate(zip(q, width, height)): 180 285 if w == 0. and h == 0.: … … 185 290 # distance to the neighbouring points. 186 291 weights[i, :] = (q_calc == qi) 187 elif h == 0 or w > 100.0 * h: 188 # Convert bin edges from q to u 189 u_limit = np.sqrt(qi**2 + w**2) 190 u_edges = q_edges**2  qi**2 191 u_edges[q_edges < qi] = 0. 192 u_edges[q_edges > u_limit] = u_limit**2  qi**2 193 weights[i, :] = np.diff(np.sqrt(u_edges))/w 194 #print "i, qi",i,qi,qi+width 195 #print q_calc 196 #print weights[i,:] 292 elif h == 0: 293 weights[i, :] = _q_perp_weights(q_edges, qi, w) 197 294 elif w == 0: 198 in_x = (q_calc >= qih) * (q_calc <= qi+h) 199 weights[i,:] = in_x * np.diff(q_edges) / (2*h) 295 in_x = 1.0 * ((q_calc >= qih) & (q_calc <= qi+h)) 296 abs_x = 1.0*(q_calc < abs(qi  h)) if qi < h else 0. 297 #print qi  h, qi + h 298 #print in_x + abs_x 299 weights[i,:] = (in_x + abs_x) * np.diff(q_edges) / (2*h) 300 else: 301 L = n_height 302 for k in range(L, L+1): 303 weights[i,:] += _q_perp_weights(q_edges, qi+k*h/L, w) 304 weights[i,:] /= 2*L + 1 200 305 201 306 return weights.T 307 308 309 def _q_perp_weights(q_edges, qi, w): 310 # Convert bin edges from q to u 311 u_limit = np.sqrt(qi**2 + w**2) 312 u_edges = q_edges**2  qi**2 313 u_edges[q_edges < abs(qi)] = 0. 314 u_edges[q_edges > u_limit] = u_limit**2  qi**2 315 weights = np.diff(np.sqrt(u_edges))/w 316 #print "i, qi",i,qi,qi+width 317 #print q_calc 318 #print weights 319 return weights 202 320 203 321 … … 218 336 function. 219 337 """ 220 q_min, q_max = np.min(q), np.max(np.sqrt(q**2 + width**2)) 338 q_min, q_max = np.min(qheight), np.max(np.sqrt((q+height)**2 + width**2)) 339 221 340 return geometric_extrapolation(q, q_min, q_max) 222 341 … … 270 389 q = np.sort(q) 271 390 if q_min < q[0]: 272 if q_min <= 0: q_min = q [0]/10391 if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 273 392 n_low = np.ceil((q[0]q_min) / (q[1]q[0])) if q[1]>q[0] else 15 274 393 q_low = np.linspace(q_min, q[0], n_low+1)[:1] … … 320 439 log_delta_q = log(10.) / points_per_decade 321 440 if q_min < q[0]: 322 if q_min < 0: q_min = q[0] /10441 if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 323 442 n_low = log_delta_q * (log(q[0])log(q_min)) 324 443 q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:1] … … 373 492 _int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars) 374 493 _int_h = lambda h, qi: eval_form(abs(qi+h), form, pars) 375 _int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars) 494 # If both width and height are defined, then it is too slow to use dblquad. 495 # Instead use trapz on a fixed grid, interpolated into the I(Q) for 496 # the extended Q range. 497 #_int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars) 498 q_calc = slit_extend_q(q, np.asarray(width), np.asarray(height)) 499 Iq = eval_form(q_calc, form, pars) 376 500 result = np.empty(len(q)) 377 501 for i, (qi, w, h) in enumerate(zip(q, width, height)): … … 385 509 result[i] = r/(2*h) 386 510 else: 387 r, err = dblquad(_int_wh, h, h, lambda h: 0., lambda h: width, 511 w_grid = np.linspace(0, w, 21)[None,:] 512 h_grid = np.linspace(h, h, 23)[:,None] 513 u = sqrt((qi+h_grid)**2 + w_grid**2) 514 Iu = np.interp(u, q_calc, Iq) 515 #print np.trapz(Iu, w_grid, axis=1) 516 Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:,0]) 517 result[i] = Is / (2*h*w) 518 """ 519 r, err = dblquad(_int_wh, h, h, lambda h: 0., lambda h: w, 388 520 args=(qi,)) 389 521 result[i] = r/(w*2*h) 522 """ 390 523 391 524 # r should be [float, ...], but it is [array([float]), array([float]),...] … … 920 1053 921 1054 def demo_slit_1d(): 922 q = np.logspace(4, np.log10(0.2), 400)1055 q = np.logspace(4, np.log10(0.2), 100) 923 1056 w = h = 0. 924 #w = 0.00 51057 #w = 0.000000277790 925 1058 w = 0.0277790 926 1059 #h = 0.00277790 1060 #h = 0.0277790 927 1061 resolution = Slit1D(q, w, h) 928 1062 _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w,h)) … … 940 1074 941 1075 if __name__ == "__main__": 942 #demo()943 main()1076 demo() 1077 #main()
Note: See TracChangeset
for help on using the changeset viewer.