diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py index 7662e17a..075705a1 100644 --- a/sasmodels/direct_model.py +++ b/sasmodels/direct_model.py @@ -541,7 +541,7 @@ def main(): q, = values dq = dqw = dql = None #dq = [q*0.05] # 5% pinhole resolution - #dqw, dql = [q*0.05], [1.0] # 5% horizontal slit resolution + dqw, dql = [q*0.05], [1.0] # 5% horizontal slit resolution print(Iq(model, [q], dq=dq, qw=dqw, ql=dql, **pars)[0]) #print(Gxi(model, [q], **pars)[0]) elif len(values) == 2: @@ -554,24 +554,28 @@ def main(): sys.exit(1) def test_simple_interface(): - def near(value, target): + def assert_near(value: np.ndarray, target: list[float]): """Close enough in single precision""" #print(f"value: {value}, target: {target}") - return np.allclose(value, target, rtol=1e-6, atol=0, equal_nan=True) - # Note: target values taken from running main() on parameters. + if not np.allclose(value, target, rtol=1e-6, atol=0, equal_nan=True): + assert value.tolist() == target + #raise ValueError(f"Expected {value} but got {target}") + # Target values taken from adusting main() with target resolution then running with: + # python -m sasmodels.direct_model sphere 0.1 radius=200 background=0 + # python -m sasmodels.direct_model sphere 0.1,0.1 radius=200 background=0 # Resolution was 5% dq/q. pars = dict(radius=200, background=0) # default background=1e-3, scale=1 # simple sphere in 1D (perfect, pinhole, slit) perfect_target = 0.6190146273894904 - assert near(Iq('sphere', [0.1], **pars), [perfect_target]) - assert near(Iq('sphere', [0.1], dq=[0.005], **pars), [2.3009224683980215]) - assert near(Iq('sphere', [0.1], qw=[0.005], ql=[1.0], **pars), [0.3663431784535172]) + assert_near(Iq('sphere', [0.1], **pars), [perfect_target]) + assert_near(Iq('sphere', [0.1], dq=[0.005], **pars), [2.3009224683980215]) + assert_near(Iq('sphere', [0.1], qw=[0.005], ql=[1.0], **pars), [0.1650934496236075]) # simple sphere in 2D (perfect, pinhole) - assert near(Iqxy('sphere', [0.1], [0.1], **pars), [1.1771532874802199]) - assert near(Iqxy('sphere', [0.1], [0.1], dqx=[0.005], dqy=[0.005], **pars), + assert_near(Iqxy('sphere', [0.1], [0.1], **pars), [1.1771532874802199]) + assert_near(Iqxy('sphere', [0.1], [0.1], dqx=[0.005], dqy=[0.005], **pars), [0.8167780778578667]) # sesans (no background or scale) - assert near(Gxi('sphere', [100], **pars), [-0.19146959126623486]) + assert_near(Gxi('sphere', [100], **pars), [-0.19146959126623486]) # Check that single point sesans matches value in an array xi = np.logspace(1, 3, 100) y = Gxi('sphere', xi, **pars) @@ -581,16 +585,16 @@ def near(value, target): assert abs((ysingle-y[k])/y[k]) < 0.1, "SESANS point value not matching vector value within 10%" # magnetic 2D pars = dict(radius=200, sld_M0=3, sld_mtheta=30) - assert near(Iqxy('sphere', [0.1], [0.1], **pars), [1.5577852226925908]) + assert_near(Iqxy('sphere', [0.1], [0.1], **pars), [1.5577852226925908]) # polydisperse 1D pars = dict( radius=200, radius_pd=0.1, radius_pd_n=15, radius_pd_nsigma=2.5, radius_pd_type="uniform") - assert near(Iq('sphere', [0.1], **pars), [2.703169824954617]) + assert_near(Iq('sphere', [0.1], **pars), [2.703169824954617]) # background and scale background, scale = 1e-4, 0.1 pars = dict(radius=200, background=background, scale=scale) - assert near(Iq('sphere', [0.1], **pars), [perfect_target*scale + background]) + assert_near(Iq('sphere', [0.1], **pars), [perfect_target*scale + background]) if __name__ == "__main__": diff --git a/sasmodels/kernelcl.py b/sasmodels/kernelcl.py index 85458fad..76500360 100644 --- a/sasmodels/kernelcl.py +++ b/sasmodels/kernelcl.py @@ -443,7 +443,7 @@ def __setstate__(self, state): def make_kernel(self, q_vectors): # type: (list[np.ndarray]) -> "GpuKernel" - return GpuKernel(self, q_vectors) + return GpuKernel(self, np.asarray(q_vectors, dtype=self.dtype)) def get_function(self, name): # type: (str) -> cl.Kernel diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py index 67e363e4..d5df730c 100644 --- a/sasmodels/resolution.py +++ b/sasmodels/resolution.py @@ -8,6 +8,7 @@ import numpy as np # type: ignore from numpy import exp, log, log10, pi, sqrt # type: ignore +from numpy.typing import NDArray from scipy.special import erf # type: ignore __all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D", @@ -37,8 +38,10 @@ class Resolution: """ - q = None # type: np.ndarray - q_calc = None # type: np.ndarray + q: NDArray + q_calc: NDArray + weight_matrix: NDArray + def apply(self, theory): """ Smear *theory* by the resolution function, returning *Iq*. @@ -54,6 +57,7 @@ class Perfect1D(Resolution): """ def __init__(self, q): self.q_calc = self.q = q + # No weight matrix in Perfect1D def apply(self, theory): return theory @@ -122,7 +126,8 @@ class Slit1D(Resolution): The *weight_matrix* is computed by :func:`slit_resolution` """ - + # TODO: change to (self, q, *, ...) to force use of q_length/q_width keywords + # Note: requires change in SasView def __init__(self, q, q_length=None, q_width=None, q_calc=None): # Remember what width/dqy was used even though we won't need them # after the weight matrix is constructed @@ -145,7 +150,7 @@ def __init__(self, q, q_length=None, q_width=None, q_calc=None): self.q = q.flatten() self.q_calc = ( - slit_extend_q(q, q_width, q_length) + slit_extend_q(q, q_length=q_length, q_width=q_width) if q_calc is None else np.sort(q_calc) ) @@ -156,7 +161,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_length=q_length, q_width=q_width) ) self.q_calc = abs(self.q_calc) @@ -215,15 +220,15 @@ def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA): return weights -def slit_resolution(q_calc, q, width, length, n_length=30): +def slit_resolution(q_calc, q, *, q_length=0., q_width=0., n_length=30): r""" Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given - $q_\perp$ = *width* (in the high-resolution axis) and $q_\parallel$ - = *length* (in the low resolution axis). *n_length* is the number + $q_\perp$ = *q_width* (in the high-resolution axis) and $q_\parallel$ + = *q_length* (in the low resolution axis). *n_length* is the number of steps to use in the integration over $q_\parallel$ when both $q_\perp$ and $q_\parallel$ are non-zero. - Each $q$ can have an independent width and length value even though + Each $q$ can have an independent width and length resolution even though current instruments use the same slit setting for all measured points. If slit length is large relative to width, use: @@ -281,7 +286,7 @@ def slit_resolution(q_calc, q, width, length, n_length=30): where $I(u_j)$ is the value at the mid-point, and $\Delta u_j$ is the difference between consecutive edges which have been first converted to $u$. Only $u_j \in [0, \Delta q_\perp]$ are used, which corresponds - to $q_j \in \left[q, \sqrt{q^2 + \Delta q_\perp}\right]$, so + to $q_j \in \left[q, \sqrt{q^2 + \Delta q_\perp^2}\right]$, so .. math:: @@ -338,7 +343,6 @@ def slit_resolution(q_calc, q, width, length, n_length=30): """ - # np.set_printoptions(precision=6, linewidth=10000) # The current algorithm is a midpoint rectangle rule. @@ -347,40 +351,44 @@ 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, q_width, q_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) + # TODO: need pixel fraction for the boundary pixels + # Diff q_edges in the expression below gives the pixel width in + # q. As wi drops below the width of a single pixel, the + # numerator stays the same but the denominater gets smaller, leading + # to ever increasing weights. We can't see this in SasView since + # it explicitly excludes q_length < 10 q_width. + 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, q_length): + # TODO: likely need pixel fraction for the boundary pixels # Convert bin edges from q to u - u_limit = np.sqrt(qi**2 + w**2) + u_limit = np.sqrt(qi**2 + q_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 - #print("i, qi",i,qi,qi+width) - #print(q_calc) - #print(weights) + weights = np.diff(np.sqrt(u_edges))/q_length return weights @@ -398,13 +406,13 @@ def pinhole_extend_q(q, q_width, nsigma=PINHOLE_N_SIGMA): return linear_extrapolation(q, q_min, q_max) -def slit_extend_q(q, width, length): +def slit_extend_q(q, *, q_length=0, q_width=0): """ - Given *q*, *width* and *length*, find a set of sampling points *q_calc* so + Given *q*, *q_width* and *q_length*, find a set of sampling points *q_calc* so that each point I(q) has sufficient support from the underlying function. """ - q_min, q_max = np.min(q-length), np.max(np.sqrt((q+length)**2 + width**2)) + q_min, q_max = np.min(q-q_width), np.max(np.sqrt((q+q_width)**2 + q_length**2)) return geometric_extrapolation(q, q_min, q_max) @@ -509,7 +517,7 @@ def geometric_extrapolation(q, q_min, q_max, points_per_decade=None): data_min, data_max = q[0], q[-1] if points_per_decade is None: if data_max > data_min: - log_delta_q = (len(q) - 1) / (log(data_max) - log(data_min)) + log_delta_q = (log(data_max) - log(data_min)) / (len(q) - 1) else: log_delta_q = log(10.) / DEFAULT_POINTS_PER_DECADE else: @@ -517,12 +525,12 @@ def geometric_extrapolation(q, q_min, q_max, points_per_decade=None): if q_min < data_min: if q_min < 0: q_min = data_min*MINIMUM_ABSOLUTE_Q - n_low = int(np.ceil(log_delta_q * (log(q[0])-log(q_min)))) + n_low = int(np.ceil((log(q[0])-log(q_min)) / log_delta_q)) q_low = np.logspace(log10(q_min), log10(q[0]), n_low+1)[:-1] else: q_low = [] if q_max > data_max: - n_high = int(np.ceil(log_delta_q * (log(q_max)-log(data_max)))) + n_high = int(np.ceil((log(q_max)-log(data_max)) / log_delta_q)) q_high = np.logspace(log10(data_max), log10(q_max), n_high+1)[1:] else: q_high = [] @@ -542,7 +550,7 @@ def eval_form(q, form, pars): *pars* are the parameter values to use when evaluating. """ from sasmodels import direct_model - kernel = form.make_kernel([q]) + kernel = form.make_kernel(np.array([q])) theory = direct_model.call_kernel(kernel, pars) kernel.release() return theory @@ -560,50 +568,56 @@ def gaussian(q, q0, dq, nsigma=2.5): return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density) -def _quad_slit_1d(q, width, length, form, pars): +def _quad_slit_1d(form, pars, q, *, q_width=0., q_length=0.): """ Scipy integration for slit resolution. This is an adaptive integration technique. It is called with settings that make it slow to evaluate but give it good accuracy. """ - from scipy.integrate import quad # type: ignore par_set = {p.name for p in form.info.parameters.call_parameters} - if any(k not in par_set for k in pars.keys()): + if any(k not in par_set for k in pars): extra = set(pars.keys()) - par_set raise ValueError("bad parameters: [%s] not in [%s]" % (", ".join(sorted(extra)), ", ".join(sorted(pars.keys())))) - if np.isscalar(width): - 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)[0] - _int_l = lambda l, qi: eval_form(abs(qi+l), form, pars)[0] - # 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)) + if np.isscalar(q_width): + q_width = np.full_like(q, q_width) + if np.isscalar(q_length): + q_length = np.full_like(q, q_length) + #def _int_l(lk, qo): return eval_form(sqrt(qo**2 + lk**2), form, pars)[0] + #def _int_w(wk, qo): return eval_form(abs(qo+wk), form, pars)[0] + q_min, q_max = (q-q_width).min(), sqrt((q+q_width).max()**2 + q_length.max()**2) + q_calc = np.linspace(q_min, q_max, 10000) 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, err = quad(_int_w, 0, w, args=(qi,), epsabs=0, epsrel=1e-8) - result[i] = total/w - elif w == 0.: - total, err = quad(_int_l, -l, l, args=(qi,), epsabs=0, epsrel=1e-8) - result[i] = total/(2*l) + nw, nl = 121, 251 + for i, (qi, dw, dl) in enumerate(zip(q, q_width, q_length)): + # print(f"slit {qi} with {dw=}, {dl=}") + if dl == 0.: + #total, err = quad(_int_w, 0, dw, args=(qi,), epsabs=0, epsrel=1e-8) + w_grid = np.linspace(-dw, dw, nw) + u_sub = qi + w_grid + f_at_u = np.interp(u_sub, q_calc, Iq) + total = np.trapezoid(f_at_u, w_grid) + result[i] = total/(2*dw) + elif dw == 0.: + #total, err = quad(_int_l, -dl, dl, args=(qi,), epsabs=0, epsrel=1e-8) + l_grid = np.linspace(-dl, dl, nl) + u_sub = sqrt(qi**2 + l_grid**2) + f_at_u = np.interp(u_sub, q_calc, Iq) + total = np.trapezoid(f_at_u, l_grid) + result[i] = total/(2*dl) else: - w_grid = np.linspace(0, w, 21)[None, :] - l_grid = np.linspace(-l, l, 23)[:, None] - u_sub = sqrt((qi+l_grid)**2 + w_grid**2) + l_grid = np.linspace(-dl, dl, nl)[None, :] + w_grid = np.linspace(-dw, dw, nw)[:, None] + u_sub = sqrt((qi + w_grid)**2 + l_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) + total = np.trapezoid(np.trapezoid(f_at_u, l_grid, axis=1), w_grid[:, 0]) + result[i] = total / (4*dl*dw) # from scipy.integrate import dblquad # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, # args=(qi,)) @@ -613,7 +627,7 @@ def _quad_slit_1d(q, width, length, form, pars): return result -def _quad_pinhole_1d(q, q_width, form, pars, nsigma=2.5): +def _quad_pinhole_1d(form, pars, q, q_width, nsigma=2.5): """ Scipy integration for pinhole resolution. @@ -623,16 +637,17 @@ def _quad_pinhole_1d(q, q_width, form, pars, nsigma=2.5): from scipy.integrate import quad # type: ignore par_set = {p.name for p in form.info.parameters.call_parameters} - if any(k not in par_set for k in pars.keys()): + if any(k not in par_set for k in pars): extra = set(pars.keys()) - par_set raise ValueError("bad parameters: [%s] not in [%s]" % (", ".join(sorted(extra)), ", ".join(sorted(pars.keys())))) - func = lambda q, q0, dq: eval_form(q, form, pars)[0]*gaussian(q, q0, dq)[0] - total = [quad(func, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi, - args=(qi, dqi), epsabs=0, epsrel=1e-8)[0] - for qi, dqi in zip(q, q_width)] + def func(q, q0, dq): return eval_form(q, form, pars)[0]*gaussian(q, q0, dq) + total = [ + quad(func, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi, + args=(qi, dqi), epsabs=0, epsrel=1e-8)[0] + for qi, dqi in zip(q, q_width)] return np.asarray(total).flatten() @@ -670,49 +685,57 @@ def test_slit_zero(self): @unittest.skip("not yet supported") def test_slit_long(self): """ - Slit smearing with length 0.005 + Slit smearing with q_length=0.005 """ resolution = Slit1D(self.x, q_width=0, q_length=0.005, q_calc=self.x) theory = self.Iq(resolution.q_calc) output = resolution.apply(theory) + # TODO calculate expected result rather than using program output answer = [ - 9.0618, 8.6402, 8.1187, 7.1392, 6.1528, - 5.5555, 4.5584, 3.5606, 2.5623, 2.0000, + 9.230181, 8.680682, 7.95333 , 7.167294, 6.28892 , 5.4 , + 4.502882, 3.574456, 2.608276, 1.280625, ] np.testing.assert_allclose(output, answer, atol=1e-4) @unittest.skip("not yet supported") def test_slit_both_high(self): """ - Slit smearing with width < 100*length. + Slit smearing with q_width < 100*q_length. """ q = np.logspace(-4, -1, 10) - resolution = Slit1D(q, q_width=0.2, q_length=np.inf) - theory = 1000*self.Iq(resolution.q_calc**4) + q_width = 0.2*q # 20% ΔQ radial + resolution = Slit1D(q, q_width=q_width, q_length=3.0) + theory = 1000*self.Iq(resolution.q_calc)**-4 output = resolution.apply(theory) + # TODO calculate expected result rather than using program output answer = [ - 8.85785, 8.43012, 7.92687, 6.94566, 6.03660, - 5.40363, 4.40655, 3.40880, 2.41058, 2.00000, - ] + 1.773625e-01, 1.773912e-01, 1.775246e-01, 1.781470e-01, + 1.811408e-01, 1.981946e-01, 2.529219e-01, 1.058098e-03, + 1.386471e-05, 6.579479e-07, + ] np.testing.assert_allclose(output, answer, atol=1e-4) + # Broken code. See sasmodels #697 @unittest.skip("not yet supported") def test_slit_wide(self): """ - Slit smearing with width 0.0002 + Slit smearing with q_width 0.0002 """ - resolution = Slit1D(self.x, q_width=0.0002, q_length=0, q_calc=self.x) + resolution = Slit1D(self.x, q_width=0.000002, q_length=0, q_calc=self.x) theory = self.Iq(resolution.q_calc) output = resolution.apply(theory) + #print("q", resolution.q_calc) + #print("I(q)", theory) + #print("w", resolution.weight_matrix) answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, - ] + ] np.testing.assert_allclose(output, answer, atol=1e-4) @unittest.skip("not yet supported") def test_slit_both_wide(self): """ - Slit smearing with width > 100*length. + Slit smearing with q_width > 100*q_length. """ resolution = Slit1D(self.x, q_width=0.0002, q_length=0.000001, q_calc=self.x) @@ -817,8 +840,8 @@ def test_pinhole_scipy(self): pars['radius'] *= 5 data = np.loadtxt(data_string.split('\n')).T - q, q_width, answer = data - answer = _quad_pinhole_1d(q, q_width, self.model, pars) + q, q_width, _answer = data + answer = _quad_pinhole_1d(self.model, pars, q, q_width) ## Getting 0.1% requires 5 sigma and 200 points per fringe #q_calc = interpolate(pinhole_extend_q(q, q_width, nsigma=5), # 2*np.pi/pars['radius']/200) @@ -833,8 +856,8 @@ def test_pinhole_scipy(self): resolution = Perfect1D(q) source = self._eval_sphere(pars, resolution) plt.loglog(q, source, '.') - plt.loglog(q, answer, '-', hold=True) - plt.loglog(q, output, '-', hold=True) + plt.loglog(q, answer, '-') + plt.loglog(q, output, '-') plt.show() self._compare(q, output, answer, tol) @@ -863,11 +886,12 @@ def test_slit_scipy(self): data_string = TEST_DATA_SLIT_SPHERE data = np.loadtxt(data_string.split('\n')).T - q, delta_qv, _, answer = data - answer = _quad_slit_1d(q, delta_qv, 0., self.model, pars) - q_calc = slit_extend_q(interpolate(q, 2*np.pi/pars['radius']/20), - delta_qv[0], 0.) - resolution = Slit1D(q, q_length=delta_qv, q_width=0, q_calc=q_calc) + q, dqv, _, _answer = data + answer = _quad_slit_1d(self.model, pars, q, q_length=dqv) + q_calc = slit_extend_q( + interpolate(q, 2*np.pi/pars['radius']/20), + q_length=dqv[0], q_width=0.) + resolution = Slit1D(q, q_length=dqv, q_width=0, q_calc=q_calc) output = self._eval_sphere(pars, resolution) # TODO: relative error should be lower self._compare(q, output, answer, 0.025) @@ -885,9 +909,9 @@ def test_slit_ellipsoid_scipy(self): } form = load_model('ellipsoid', dtype='double') q = np.logspace(log10(4e-5), log10(2.5e-2), 68) - width, length = 0.,0.117 - resolution = Slit1D(q, q_length=length, q_width=width) - answer = _quad_slit_1d(q, length, width, form, pars) + q_width, q_length = 0.,0.117 + resolution = Slit1D(q, q_length=q_length, q_width=q_width) + answer = _quad_slit_1d(form, pars, q, q_length=q_length, q_width=q_width) output = resolution.apply(eval_form(resolution.q_calc, form, pars)) # TODO: 10% is too much error; use better algorithm #print(np.max(abs(answer-output)/answer)) @@ -1127,12 +1151,14 @@ def _eval_demo_1d(resolution, title): pars = {'length':210, 'radius':500, 'background': 0} elif name == 'teubner_strey': pars = {'a2':0.003, 'c1':-1e4, 'c2':1e10, 'background':0.312643} - elif name in ('sphere', 'spherepy'): - pars = TEST_PARS_SLIT_SPHERE + elif name == 'sphere': + #pars = TEST_PARS_SLIT_SPHERE + # 20 μm sphere needs lower q range + pars = {'radius': 200, 'background': 0, 'sld': 4, 'sld_solvent': 1, 'scale': 0.01} elif name == 'ellipsoid': pars = { 'scale':0.05, 'background': 0, - 'r_polar':500, 'r_equatorial':15000, + 'radius_polar':500, 'radius_equatorial':15000, 'sld':6, 'sld_solvent': 1, } else: @@ -1145,57 +1171,60 @@ def _eval_demo_1d(resolution, title): Iq = resolution.apply(theory) if isinstance(resolution, Slit1D): - width, length = resolution.q_width, resolution.q_length - Iq_romb = _quad_slit_1d(resolution.q, width, length, model, pars) + q, q_width, q_length = resolution.q, resolution.q_width, resolution.q_length + index = slice(None, None, len(q)//25) + q_target = q[index] + Iq_target = _quad_slit_1d(model, pars, q_target, q_width=q_width[index], q_length=q_length) else: - dq = resolution.q_width - Iq_romb = _quad_pinhole_1d(resolution.q, dq, model, pars) + q, dq = resolution.q, resolution.q_width + index = slice(None, None, len(q)//25) + q_target = q[index] + Iq_target = _quad_pinhole_1d(model, pars, q_target, dq[index]) import matplotlib.pyplot as plt # type: ignore plt.loglog(resolution.q_calc, theory, label='unsmeared') - plt.loglog(resolution.q, Iq, label='smeared', hold=True) - plt.loglog(resolution.q, Iq_romb, label='scipy smeared', hold=True) + plt.loglog(resolution.q, Iq, label='smeared') + plt.loglog(q[index], Iq_target, 'o', label='scipy smeared') plt.legend() plt.title(title) plt.xlabel("Q (1/Ang)") plt.ylabel("I(Q) (1/cm)") -def demo_pinhole_1d(): +def demo_pinhole_1d(dQoQ=0.1): """ Show example of pinhole smearing. """ - q = np.logspace(-4, np.log10(0.2), 400) - q_width = 0.1*q + q = np.logspace(-4, np.log10(0.2), 4000) + q_width = dQoQ*q resolution = Pinhole1D(q, q_width) - _eval_demo_1d(resolution, title="10% dQ/Q Pinhole Resolution") + _eval_demo_1d(resolution, title=f"Pinhole Resolution ({100*dQoQ}% dQ/Q)") -def demo_slit_1d(): +def demo_slit_1d(dQoQ=0., q_length=0.): """ 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)) + q = np.logspace(-4, np.log10(0.2), 4000) + resolution = Slit1D(q, q_width=dQoQ*q*sqrt(3), q_length=q_length) + _eval_demo_1d(resolution, title=f"Slit Resolution (q_length={q_length},q_width={100*dQoQ}%)") def demo(): """ Run the resolution demos. """ import matplotlib.pyplot as plt # type: ignore - plt.subplot(121) - demo_pinhole_1d() + plt.subplot(221) + demo_pinhole_1d(dQoQ=0.05) #plt.yscale('linear') - plt.subplot(122) - demo_slit_1d() + plt.subplot(222) + demo_slit_1d(q_length=0.03) + plt.subplot(223) + demo_slit_1d(dQoQ=0.05) + plt.subplot(224) + demo_slit_1d(q_length=0.03, dQoQ=0.05) #plt.yscale('linear') plt.show() if __name__ == "__main__": - #demo() - main() + demo() + #main() # main() runs the unit tests