-
Notifications
You must be signed in to change notification settings - Fork 30
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
fixed switched nomenclature for length and width #533
base: master
Are you sure you want to change the base?
Changes from 4 commits
ae168e5
fd03bda
a6d8ec0
81ffb2e
23b05a9
f07dc5d
5bdfdce
a2b169e
17010b7
563ad5b
d0dfa73
0e0b63d
aa2b7b2
e57e776
98053f4
2fc80e0
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change | ||
---|---|---|---|---|
|
@@ -157,7 +157,7 @@ def __init__(self, q, q_length=None, q_width=None, q_calc=None): | |||
|
||||
# Build weight matrix from calculated q values | ||||
self.weight_matrix = ( | ||||
slit_resolution(self.q_calc, self.q, q_length, q_width) | ||||
slit_resolution(self.q_calc, self.q, q_width, q_length) | ||||
) | ||||
self.q_calc = abs(self.q_calc) | ||||
|
||||
|
@@ -339,7 +339,6 @@ def slit_resolution(q_calc, q, width, length, n_length=30): | |||
|
||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. On line 85: sasmodels/sasmodels/resolution.py Line 285 in 5bdfdce
the There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I've corrected the There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. At this time I have not done a thorough read of the doc string to check for additional typos in the equations. |
||||
|
||||
""" | ||||
|
||||
# np.set_printoptions(precision=6, linewidth=10000) | ||||
|
||||
# The current algorithm is a midpoint rectangle rule. | ||||
|
@@ -348,37 +347,37 @@ def slit_resolution(q_calc, q, width, length, n_length=30): | |||
weights = np.zeros((len(q), len(q_calc)), 'd') | ||||
|
||||
#print(q_calc) | ||||
for i, (qi, w, l) in enumerate(zip(q, width, length)): | ||||
if w == 0. and l == 0.: | ||||
for i, (qi, wi, li) in enumerate(zip(q, width, length)): | ||||
if wi == 0. and li == 0.: | ||||
# Perfect resolution, so return the theory value directly. | ||||
# Note: assumes that q is a subset of q_calc. If qi need not be | ||||
# in q_calc, then we can do a weighted interpolation by looking | ||||
# up qi in q_calc, then weighting the result by the relative | ||||
# distance to the neighbouring points. | ||||
weights[i, :] = (q_calc == qi) | ||||
elif l == 0: | ||||
weights[i, :] = _q_perp_weights(q_edges, qi, w) | ||||
elif w == 0: | ||||
in_x = 1.0 * ((q_calc >= qi-l) & (q_calc <= qi+l)) | ||||
abs_x = 1.0*(q_calc < abs(qi - l)) if qi < l else 0. | ||||
#print(qi - l, qi + l) | ||||
elif wi == 0: | ||||
weights[i, :] = _q_perp_weights(q_edges, qi, li) | ||||
elif li == 0: | ||||
in_x = 1.0 * ((q_calc >= qi-wi) & (q_calc <= qi+wi)) | ||||
abs_x = 1.0*(q_calc < abs(qi - wi)) if qi < wi else 0. | ||||
#print(qi - w, qi + w) | ||||
#print(in_x + abs_x) | ||||
weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*l) | ||||
weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*wi) | ||||
else: | ||||
for k in range(-n_length, n_length+1): | ||||
weights[i, :] += _q_perp_weights(q_edges, qi+k*l/n_length, w) | ||||
weights[i, :] += _q_perp_weights(q_edges, qi+k*wi/n_length, li) | ||||
weights[i, :] /= 2*n_length + 1 | ||||
|
||||
return weights.T | ||||
|
||||
|
||||
def _q_perp_weights(q_edges, qi, w): | ||||
def _q_perp_weights(q_edges, qi, length): | ||||
# Convert bin edges from q to u | ||||
u_limit = np.sqrt(qi**2 + w**2) | ||||
u_limit = np.sqrt(qi**2 + length**2) | ||||
u_edges = q_edges**2 - qi**2 | ||||
u_edges[q_edges < abs(qi)] = 0. | ||||
u_edges[q_edges > u_limit] = u_limit**2 - qi**2 | ||||
weights = np.diff(np.sqrt(u_edges))/w | ||||
weights = np.diff(np.sqrt(u_edges))/length | ||||
#print("i, qi",i,qi,qi+width) | ||||
#print(q_calc) | ||||
#print(weights) | ||||
|
@@ -581,32 +580,32 @@ def romberg_slit_1d(q, width, length, form, pars): | |||
width = [width]*len(q) | ||||
if np.isscalar(length): | ||||
length = [length]*len(q) | ||||
_int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars) | ||||
_int_l = lambda l, qi: eval_form(abs(qi+l), form, pars) | ||||
_int_w = lambda wi, qi: eval_form(sqrt(qi**2 + wi**2), form, pars) | ||||
_int_l = lambda li, qi: eval_form(abs(qi+li), form, pars) | ||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The I think we should change the There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I updated the romberg_slit_1d function to use length, width rather than width, length and reversed the use of length and width within the function. However, this led me to a similar change for slit_extend_q that was still using width, length. It appeared as though slit_extend_q was called in different ways so I've tried to correct this in the appropriate places. The tests that are still implemented are passing, but this will need careful review. |
||||
# If both width and length are defined, then it is too slow to use dblquad. | ||||
# Instead use trapz on a fixed grid, interpolated into the I(Q) for | ||||
# the extended Q range. | ||||
#_int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars) | ||||
q_calc = slit_extend_q(q, np.asarray(width), np.asarray(length)) | ||||
Iq = eval_form(q_calc, form, pars) | ||||
result = np.empty(len(q)) | ||||
for i, (qi, w, l) in enumerate(zip(q, width, length)): | ||||
if l == 0.: | ||||
total = romberg(_int_w, 0, w, args=(qi,), | ||||
for i, (qi, wi, li) in enumerate(zip(q, width, length)): | ||||
if li == 0.: | ||||
total = romberg(_int_w, 0, wi, args=(qi,), | ||||
divmax=100, vec_func=True, tol=0, rtol=1e-8) | ||||
result[i] = total/w | ||||
elif w == 0.: | ||||
total = romberg(_int_l, -l, l, args=(qi,), | ||||
result[i] = total/wi | ||||
elif wi == 0.: | ||||
total = romberg(_int_l, -li, li, args=(qi,), | ||||
divmax=100, vec_func=True, tol=0, rtol=1e-8) | ||||
result[i] = total/(2*l) | ||||
result[i] = total/(2*li) | ||||
else: | ||||
w_grid = np.linspace(0, w, 21)[None, :] | ||||
l_grid = np.linspace(-l, l, 23)[:, None] | ||||
w_grid = np.linspace(0, wi, 21)[None, :] | ||||
l_grid = np.linspace(-li, li, 23)[:, None] | ||||
u_sub = sqrt((qi+l_grid)**2 + w_grid**2) | ||||
f_at_u = np.interp(u_sub, q_calc, Iq) | ||||
#print(np.trapz(Iu, w_grid, axis=1)) | ||||
total = np.trapz(np.trapz(f_at_u, w_grid, axis=1), l_grid[:, 0]) | ||||
result[i] = total / (2*l*w) | ||||
result[i] = total / (2*li*wi) | ||||
# from scipy.integrate import dblquad | ||||
# r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, | ||||
# args=(qi,)) | ||||
|
@@ -1175,13 +1174,13 @@ def demo_slit_1d(): | |||
Show example of slit smearing. | ||||
""" | ||||
q = np.logspace(-4, np.log10(0.2), 100) | ||||
w = l = 0. | ||||
#w = 0.000000277790 | ||||
w = 0.0277790 | ||||
#l = 0.00277790 | ||||
#l = 0.0277790 | ||||
resolution = Slit1D(q, w, l) | ||||
_eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, l)) | ||||
width = length = 0. | ||||
#width = 0.000000277790 | ||||
width = 0.0277790 | ||||
#length = 0.00277790 | ||||
#length = 0.0277790 | ||||
resolution = Slit1D(q, width, length) | ||||
_eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(width, length)) | ||||
|
||||
def demo(): | ||||
""" | ||||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Rather than reversing the order of the parameters here, please reverse the order of the parameters in the slit_resolution method. That way
Slit1D
(the external API) andslit_resolution
will have the same order and future swapping errors a little less likely.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Switched back to q_length, q_width here and updated slit_resolution