source: sasmodels/sasmodels/resolution.py @ 2b2d431

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 2b2d431 was 2b2d431, checked in by Paul Kienzle <pkienzle@…>, 4 years ago

eliminate negative q values from the pinhole resolution distribution

  • Property mode set to 100644
File size: 41.5 KB
Line 
1"""
2Define the resolution functions for the data.
3
4This defines classes for 1D and 2D resolution calculations.
5"""
6from __future__ import division
7
8import unittest
9
10from scipy.special import erf  # type: ignore
11from numpy import sqrt, log, log10, exp, pi  # type: ignore
12import numpy as np  # type: ignore
13
14__all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D",
15           "apply_resolution_matrix", "pinhole_resolution", "slit_resolution",
16           "pinhole_extend_q", "slit_extend_q", "bin_edges",
17           "interpolate", "linear_extrapolation", "geometric_extrapolation",
18          ]
19
20MINIMUM_RESOLUTION = 1e-8
21MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data
22PINHOLE_N_SIGMA = 2.5 # From: Barker & Pedersen 1995 JAC
23
24class Resolution(object):
25    """
26    Abstract base class defining a 1D resolution function.
27
28    *q* is the set of q values at which the data is measured.
29
30    *q_calc* is the set of q values at which the theory needs to be evaluated.
31    This may extend and interpolate the q values.
32
33    *apply* is the method to call with I(q_calc) to compute the resolution
34    smeared theory I(q).
35    """
36    q = None  # type: np.ndarray
37    q_calc = None  # type: np.ndarray
38    def apply(self, theory):
39        """
40        Smear *theory* by the resolution function, returning *Iq*.
41        """
42        raise NotImplementedError("Subclass does not define the apply function")
43
44
45class Perfect1D(Resolution):
46    """
47    Resolution function to use when there is no actual resolution smearing
48    to be applied.  It has the same interface as the other resolution
49    functions, but returns the identity function.
50    """
51    def __init__(self, q):
52        self.q_calc = self.q = q
53
54    def apply(self, theory):
55        return theory
56
57
58class Pinhole1D(Resolution):
59    r"""
60    Pinhole aperture with q-dependent gaussian resolution.
61
62    *q* points at which the data is measured.
63
64    *q_width* gaussian 1-sigma resolution at each data point.
65
66    *q_calc* is the list of points to calculate, or None if this should
67    be estimated from the *q* and *q_width*.
68
69    *nsigma* is the width of the resolution function.  Should be 2.5.
70    See :func:`pinhole_resolution` for details.
71    """
72    def __init__(self, q, q_width, q_calc=None, nsigma=PINHOLE_N_SIGMA):
73        #*min_step* is the minimum point spacing to use when computing the
74        #underlying model.  It should be on the order of
75        #$\tfrac{1}{10}\tfrac{2\pi}{d_\text{max}}$ to make sure that fringes
76        #are computed with sufficient density to avoid aliasing effects.
77
78        # Protect against calls with q_width=0.  The extend_q function will
79        # not extend the q if q_width is 0, but q_width must be non-zero when
80        # constructing the weight matrix to avoid division by zero errors.
81        # In practice this should never be needed, since resolution should
82        # default to Perfect1D if the pinhole geometry is not defined.
83        self.q, self.q_width = q, q_width
84        self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma)
85                       if q_calc is None else np.sort(q_calc))
86
87        # Protect against models which are not defined for very low q.  Limit
88        # the smallest q value evaluated to 0.02*min.  Note that negative q
89        # values are trimmed even for broad resolution.  Although not possible
90        # from the geometry, they may appear since we are using a truncated
91        # gaussian to represent resolution rather than a skew distribution.
92        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
93        self.q_calc = self.q_calc[self.q_calc >= cutoff]
94
95        # Build weight matrix from calculated q values
96        self.weight_matrix = pinhole_resolution(
97            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION),
98            nsigma=nsigma)
99
100    def apply(self, theory):
101        return apply_resolution_matrix(self.weight_matrix, theory)
102
103
104class Slit1D(Resolution):
105    """
106    Slit aperture with resolution function.
107
108    *q* points at which the data is measured.
109
110    *qx_width* slit width in qx
111
112    *qy_width* slit height in qy
113
114    *q_calc* is the list of points to calculate, or None if this should
115    be estimated from the *q* and *q_width*.
116
117    The *weight_matrix* is computed by :func:`slit1d_resolution`
118    """
119    def __init__(self, q, qx_width, qy_width=0., q_calc=None):
120        # Remember what width/dqy was used even though we won't need them
121        # after the weight matrix is constructed
122        self.qx_width, self.qy_width = qx_width, qy_width
123
124        # Allow independent resolution on each point even though it is not
125        # needed in practice.
126        if np.isscalar(qx_width):
127            qx_width = np.ones(len(q))*qx_width
128        else:
129            qx_width = np.asarray(qx_width)
130        if np.isscalar(qy_width):
131            qy_width = np.ones(len(q))*qy_width
132        else:
133            qy_width = np.asarray(qy_width)
134
135        self.q = q.flatten()
136        self.q_calc = slit_extend_q(q, qx_width, qy_width) \
137            if q_calc is None else np.sort(q_calc)
138
139        # Protect against models which are not defined for very low q.  Limit
140        # the smallest q value evaluated (in absolute) to 0.02*min
141        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
142        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff]
143
144        # Build weight matrix from calculated q values
145        self.weight_matrix = \
146            slit_resolution(self.q_calc, self.q, qx_width, qy_width)
147        self.q_calc = abs(self.q_calc)
148
149    def apply(self, theory):
150        return apply_resolution_matrix(self.weight_matrix, theory)
151
152
153def apply_resolution_matrix(weight_matrix, theory):
154    """
155    Apply the resolution weight matrix to the computed theory function.
156    """
157    #print("apply shapes", theory.shape, weight_matrix.shape)
158    Iq = np.dot(theory[None, :], weight_matrix)
159    #print("result shape",Iq.shape)
160    return Iq.flatten()
161
162
163def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA):
164    r"""
165    Compute the convolution matrix *W* for pinhole resolution 1-D data.
166
167    Each row *W[i]* determines the normalized weight that the corresponding
168    points *q_calc* contribute to the resolution smeared point *q[i]*.  Given
169    *W*, the resolution smearing can be computed using *dot(W,q)*.
170
171    Note that resolution is limited to $\pm 2.5 \sigma$.[1]  The true resolution
172    function is a broadened triangle, and does not extend over the entire
173    range $(-\infty, +\infty)$.  It is important to impose this limitation
174    since some models fall so steeply that the weighted value in gaussian
175    tails would otherwise dominate the integral.
176
177    *q_calc* must be increasing.  *q_width* must be greater than zero.
178
179    [1] Barker, J. G., and J. S. Pedersen. 1995. Instrumental Smearing Effects
180    in Radially Symmetric Small-Angle Neutron Scattering by Numerical and
181    Analytical Methods. Journal of Applied Crystallography 28 (2): 105--14.
182    https://doi.org/10.1107/S0021889894010095.
183    """
184    # The current algorithm is a midpoint rectangle rule.  In the test case,
185    # neither trapezoid nor Simpson's rule improved the accuracy.
186    edges = bin_edges(q_calc)
187    #edges[edges < 0.0] = 0.0 # clip edges below zero
188    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :])
189    weights = cdf[1:] - cdf[:-1]
190    # Limit q range to +/- 2.5 sigma
191    weights[q_calc[:, None] < (q - nsigma*q_width)[None, :]] = 0.
192    weights[q_calc[:, None] > (q + nsigma*q_width)[None, :]] = 0.
193    weights /= np.sum(weights, axis=0)[None, :]
194    return weights
195
196
197def slit_resolution(q_calc, q, width, height, n_height=30):
198    r"""
199    Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given
200    $q_\perp$ = *width* and $q_\parallel$ = *height*.  *n_height* is
201    is the number of steps to use in the integration over $q_\parallel$
202    when both $q_\perp$ and $q_\parallel$ are non-zero.
203
204    Each $q$ can have an independent width and height value even though
205    current instruments use the same slit setting for all measured points.
206
207    If slit height is large relative to width, use:
208
209    .. math::
210
211        I_s(q_i) = \frac{1}{\Delta q_\perp}
212            \int_0^{\Delta q_\perp}
213                I\left(\sqrt{q_i^2 + q_\perp^2}\right) \,dq_\perp
214
215    If slit width is large relative to height, use:
216
217    .. math::
218
219        I_s(q_i) = \frac{1}{2 \Delta q_\parallel}
220            \int_{-\Delta q_\parallel}^{\Delta q_\parallel}
221                I\left(|q_i + q_\parallel|\right) \,dq_\parallel
222
223    For a mixture of slit width and height use:
224
225    .. math::
226
227        I_s(q_i) = \frac{1}{2 \Delta q_\parallel \Delta q_\perp}
228            \int_{-\Delta q_\parallel}^{\Delta q_\parallel}
229            \int_0^{\Delta q_\perp}
230                I\left(\sqrt{(q_i + q_\parallel)^2 + q_\perp^2}\right)
231                \,dq_\perp dq_\parallel
232
233    **Definition**
234
235    We are using the mid-point integration rule to assign weights to each
236    element of a weight matrix $W$ so that
237
238    .. math::
239
240        I_s(q) = W\,I(q_\text{calc})
241
242    If *q_calc* is at the mid-point, we can infer the bin edges from the
243    pairwise averages of *q_calc*, adding the missing edges before
244    *q_calc[0]* and after *q_calc[-1]*.
245
246    For $q_\parallel = 0$, the smeared value can be computed numerically
247    using the $u$ substitution
248
249    .. math::
250
251        u_j = \sqrt{q_j^2 - q^2}
252
253    This gives
254
255    .. math::
256
257        I_s(q) \approx \sum_j I(u_j) \Delta u_j
258
259    where $I(u_j)$ is the value at the mid-point, and $\Delta u_j$ is the
260    difference between consecutive edges which have been first converted
261    to $u$.  Only $u_j \in [0, \Delta q_\perp]$ are used, which corresponds
262    to $q_j \in \left[q, \sqrt{q^2 + \Delta q_\perp}\right]$, so
263
264    .. math::
265
266        W_{ij} = \frac{1}{\Delta q_\perp} \Delta u_j
267               = \frac{1}{\Delta q_\perp} \left(
268                    \sqrt{q_{j+1}^2 - q_i^2} - \sqrt{q_j^2 - q_i^2} \right)
269            \ \text{if}\  q_j \in \left[q_i, \sqrt{q_i^2 + q_\perp^2}\right]
270
271    where $I_s(q_i)$ is the theory function being computed and $q_j$ are the
272    mid-points between the calculated values in *q_calc*.  We tweak the
273    edges of the initial and final intervals so that they lie on integration
274    limits.
275
276    (To be precise, the transformed midpoint $u(q_j)$ is not necessarily the
277    midpoint of the edges $u((q_{j-1}+q_j)/2)$ and $u((q_j + q_{j+1})/2)$,
278    but it is at least in the interval, so the approximation is going to be
279    a little better than the left or right Riemann sum, and should be
280    good enough for our purposes.)
281
282    For $q_\perp = 0$, the $u$ substitution is simpler:
283
284    .. math::
285
286        u_j = \left|q_j - q\right|
287
288    so
289
290    .. math::
291
292        W_{ij} = \frac{1}{2 \Delta q_\parallel} \Delta u_j
293            = \frac{1}{2 \Delta q_\parallel} (q_{j+1} - q_j)
294            \ \text{if}\ q_j \in
295                \left[q-\Delta q_\parallel, q+\Delta q_\parallel\right]
296
297    However, we need to support cases were $u_j < 0$, which means using
298    $2 (q_{j+1} - q_j)$ when $q_j \in \left[0, q_\parallel-q_i\right]$.
299    This is not an issue for $q_i > q_\parallel$.
300
301    For both $q_\perp > 0$ and $q_\parallel > 0$ we perform a 2 dimensional
302    integration with
303
304    .. math::
305
306        u_{jk} = \sqrt{q_j^2 - (q + (k\Delta q_\parallel/L))^2}
307            \ \text{for}\ k = -L \ldots L
308
309    for $L$ = *n_height*.  This gives
310
311    .. math::
312
313        W_{ij} = \frac{1}{2 \Delta q_\perp q_\parallel}
314            \sum_{k=-L}^L \Delta u_{jk}
315                \left(\frac{\Delta q_\parallel}{2 L + 1}\right)
316
317
318    """
319    #np.set_printoptions(precision=6, linewidth=10000)
320
321    # The current algorithm is a midpoint rectangle rule.
322    q_edges = bin_edges(q_calc) # Note: requires q > 0
323    #q_edges[q_edges < 0.0] = 0.0 # clip edges below zero
324    weights = np.zeros((len(q), len(q_calc)), 'd')
325
326    #print(q_calc)
327    for i, (qi, w, h) in enumerate(zip(q, width, height)):
328        if w == 0. and h == 0.:
329            # Perfect resolution, so return the theory value directly.
330            # Note: assumes that q is a subset of q_calc.  If qi need not be
331            # in q_calc, then we can do a weighted interpolation by looking
332            # up qi in q_calc, then weighting the result by the relative
333            # distance to the neighbouring points.
334            weights[i, :] = (q_calc == qi)
335        elif h == 0:
336            weights[i, :] = _q_perp_weights(q_edges, qi, w)
337        elif w == 0:
338            in_x = 1.0 * ((q_calc >= qi-h) & (q_calc <= qi+h))
339            abs_x = 1.0*(q_calc < abs(qi - h)) if qi < h else 0.
340            #print(qi - h, qi + h)
341            #print(in_x + abs_x)
342            weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h)
343        else:
344            for k in range(-n_height, n_height+1):
345                weights[i, :] += _q_perp_weights(q_edges, qi+k*h/n_height, w)
346            weights[i, :] /= 2*n_height + 1
347
348    return weights.T
349
350
351def _q_perp_weights(q_edges, qi, w):
352    # Convert bin edges from q to u
353    u_limit = np.sqrt(qi**2 + w**2)
354    u_edges = q_edges**2 - qi**2
355    u_edges[q_edges < abs(qi)] = 0.
356    u_edges[q_edges > u_limit] = u_limit**2 - qi**2
357    weights = np.diff(np.sqrt(u_edges))/w
358    #print("i, qi",i,qi,qi+width)
359    #print(q_calc)
360    #print(weights)
361    return weights
362
363
364def pinhole_extend_q(q, q_width, nsigma=3):
365    """
366    Given *q* and *q_width*, find a set of sampling points *q_calc* so
367    that each point $I(q)$ has sufficient support from the underlying
368    function.
369    """
370    q_min, q_max = np.min(q - nsigma*q_width), np.max(q + nsigma*q_width)
371    return linear_extrapolation(q, q_min, q_max)
372
373
374def slit_extend_q(q, width, height):
375    """
376    Given *q*, *width* and *height*, find a set of sampling points *q_calc* so
377    that each point I(q) has sufficient support from the underlying
378    function.
379    """
380    q_min, q_max = np.min(q-height), np.max(np.sqrt((q+height)**2 + width**2))
381
382    return geometric_extrapolation(q, q_min, q_max)
383
384
385def bin_edges(x):
386    """
387    Determine bin edges from bin centers, assuming that edges are centered
388    between the bins.
389
390    Note: this uses the arithmetic mean, which may not be appropriate for
391    log-scaled data.
392    """
393    if len(x) < 2 or (np.diff(x) < 0).any():
394        raise ValueError("Expected bins to be an increasing set")
395    edges = np.hstack([
396        x[0]  - 0.5*(x[1]  - x[0]),  # first point minus half first interval
397        0.5*(x[1:] + x[:-1]),        # mid points of all central intervals
398        x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval
399        ])
400    return edges
401
402
403def interpolate(q, max_step):
404    """
405    Returns *q_calc* with points spaced at most max_step apart.
406    """
407    step = np.diff(q)
408    index = step > max_step
409    if np.any(index):
410        inserts = []
411        for q_i, step_i in zip(q[:-1][index], step[index]):
412            n = np.ceil(step_i/max_step)
413            inserts.extend(q_i + np.arange(1, n)*(step_i/n))
414        # Extend a couple of fringes beyond the end of the data
415        inserts.extend(q[-1] + np.arange(1, 8)*max_step)
416        q_calc = np.sort(np.hstack((q, inserts)))
417    else:
418        q_calc = q
419    return q_calc
420
421
422def linear_extrapolation(q, q_min, q_max):
423    """
424    Extrapolate *q* out to [*q_min*, *q_max*] using the step size in *q* as
425    a guide.  Extrapolation below uses about the same size as the first
426    interval.  Extrapolation above uses about the same size as the final
427    interval.
428
429    Note that extrapolated values may be negative.
430    """
431    q = np.sort(q)
432    if q_min + 2*MINIMUM_RESOLUTION < q[0]:
433        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15
434        q_low = np.linspace(q_min, q[0], n_low+1)[:-1]
435    else:
436        q_low = []
437    if q_max - 2*MINIMUM_RESOLUTION > q[-1]:
438        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15
439        q_high = np.linspace(q[-1], q_max, n_high+1)[1:]
440    else:
441        q_high = []
442    return np.concatenate([q_low, q, q_high])
443
444
445def geometric_extrapolation(q, q_min, q_max, points_per_decade=None):
446    r"""
447    Extrapolate *q* to [*q_min*, *q_max*] using geometric steps, with the
448    average geometric step size in *q* as the step size.
449
450    if *q_min* is zero or less then *q[0]/10* is used instead.
451
452    *points_per_decade* sets the ratio between consecutive steps such
453    that there will be $n$ points used for every factor of 10 increase
454    in *q*.
455
456    If *points_per_decade* is not given, it will be estimated as follows.
457    Starting at $q_1$ and stepping geometrically by $\Delta q$ to $q_n$
458    in $n$ points gives a geometric average of:
459
460    .. math::
461
462         \log \Delta q = (\log q_n - \log q_1) / (n - 1)
463
464    From this we can compute the number of steps required to extend $q$
465    from $q_n$ to $q_\text{max}$ by $\Delta q$ as:
466
467    .. math::
468
469         n_\text{extend} = (\log q_\text{max} - \log q_n) / \log \Delta q
470
471    Substituting:
472
473    .. math::
474
475         n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n)
476            / (\log q_n - \log q_1)
477    """
478    q = np.sort(q)
479    if points_per_decade is None:
480        log_delta_q = (len(q) - 1) / (log(q[-1]) - log(q[0]))
481    else:
482        log_delta_q = log(10.) / points_per_decade
483    if q_min < q[0]:
484        if q_min < 0:
485            q_min = q[0]*MINIMUM_ABSOLUTE_Q
486        n_low = log_delta_q * (log(q[0])-log(q_min))
487        q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1]
488    else:
489        q_low = []
490    if q_max > q[-1]:
491        n_high = log_delta_q * (log(q_max)-log(q[-1]))
492        q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:]
493    else:
494        q_high = []
495    return np.concatenate([q_low, q, q_high])
496
497
498############################################################################
499# unit tests
500############################################################################
501
502def eval_form(q, form, pars):
503    """
504    Return the SAS model evaluated at *q*.
505
506    *form* is the SAS model returned from :fun:`core.load_model`.
507
508    *pars* are the parameter values to use when evaluating.
509    """
510    from sasmodels import direct_model
511    kernel = form.make_kernel([q])
512    theory = direct_model.call_kernel(kernel, pars)
513    kernel.release()
514    return theory
515
516
517def gaussian(q, q0, dq, nsigma=2.5):
518    """
519    Return the truncated Gaussian resolution function.
520
521    *q0* is the center, *dq* is the width and *q* are the points to evaluate.
522    """
523    # Calculate the density of the tails; the resulting gaussian needs to be
524    # scaled by this amount in ordere to integrate to 1.0
525    two_tail_density = 2 * (1 + erf(-nsigma/sqrt(2)))/2
526    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density)
527
528
529def romberg_slit_1d(q, width, height, form, pars):
530    """
531    Romberg integration for slit resolution.
532
533    This is an adaptive integration technique.  It is called with settings
534    that make it slow to evaluate but give it good accuracy.
535    """
536    from scipy.integrate import romberg  # type: ignore
537
538    par_set = set([p.name for p in form.info.parameters.call_parameters])
539    if any(k not in par_set for k in pars.keys()):
540        extra = set(pars.keys()) - par_set
541        raise ValueError("bad parameters: [%s] not in [%s]"
542                         % (", ".join(sorted(extra)),
543                            ", ".join(sorted(pars.keys()))))
544
545    if np.isscalar(width):
546        width = [width]*len(q)
547    if np.isscalar(height):
548        height = [height]*len(q)
549    _int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars)
550    _int_h = lambda h, qi: eval_form(abs(qi+h), form, pars)
551    # If both width and height are defined, then it is too slow to use dblquad.
552    # Instead use trapz on a fixed grid, interpolated into the I(Q) for
553    # the extended Q range.
554    #_int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars)
555    q_calc = slit_extend_q(q, np.asarray(width), np.asarray(height))
556    Iq = eval_form(q_calc, form, pars)
557    result = np.empty(len(q))
558    for i, (qi, w, h) in enumerate(zip(q, width, height)):
559        if h == 0.:
560            total = romberg(_int_w, 0, w, args=(qi,),
561                            divmax=100, vec_func=True, tol=0, rtol=1e-8)
562            result[i] = total/w
563        elif w == 0.:
564            total = romberg(_int_h, -h, h, args=(qi,),
565                            divmax=100, vec_func=True, tol=0, rtol=1e-8)
566            result[i] = total/(2*h)
567        else:
568            w_grid = np.linspace(0, w, 21)[None, :]
569            h_grid = np.linspace(-h, h, 23)[:, None]
570            u_sub = sqrt((qi+h_grid)**2 + w_grid**2)
571            f_at_u = np.interp(u_sub, q_calc, Iq)
572            #print(np.trapz(Iu, w_grid, axis=1))
573            total = np.trapz(np.trapz(f_at_u, w_grid, axis=1), h_grid[:, 0])
574            result[i] = total / (2*h*w)
575            # from scipy.integrate import dblquad
576            # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w,
577            #                  args=(qi,))
578            # result[i] = r/(w*2*h)
579
580    # r should be [float, ...], but it is [array([float]), array([float]),...]
581    return result
582
583
584def romberg_pinhole_1d(q, q_width, form, pars, nsigma=2.5):
585    """
586    Romberg integration for pinhole resolution.
587
588    This is an adaptive integration technique.  It is called with settings
589    that make it slow to evaluate but give it good accuracy.
590    """
591    from scipy.integrate import romberg  # type: ignore
592
593    par_set = set([p.name for p in form.info.parameters.call_parameters])
594    if any(k not in par_set for k in pars.keys()):
595        extra = set(pars.keys()) - par_set
596        raise ValueError("bad parameters: [%s] not in [%s]"
597                         % (", ".join(sorted(extra)),
598                            ", ".join(sorted(pars.keys()))))
599
600    func = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
601    total = [romberg(func, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi,
602                     args=(qi, dqi), divmax=100, vec_func=True,
603                     tol=0, rtol=1e-8)
604             for qi, dqi in zip(q, q_width)]
605    return np.asarray(total).flatten()
606
607
608class ResolutionTest(unittest.TestCase):
609    """
610    Test the resolution calculations.
611    """
612
613    def setUp(self):
614        self.x = 0.001*np.arange(1, 11)
615        self.y = self.Iq(self.x)
616
617    def Iq(self, q):
618        "Linear function for resolution unit test"
619        return 12.0 - 1000.0*q
620
621    def test_perfect(self):
622        """
623        Perfect resolution and no smearing.
624        """
625        resolution = Perfect1D(self.x)
626        theory = self.Iq(resolution.q_calc)
627        output = resolution.apply(theory)
628        np.testing.assert_equal(output, self.y)
629
630    def test_slit_zero(self):
631        """
632        Slit smearing with perfect resolution.
633        """
634        resolution = Slit1D(self.x, qx_width=0, qy_width=0, q_calc=self.x)
635        theory = self.Iq(resolution.q_calc)
636        output = resolution.apply(theory)
637        np.testing.assert_equal(output, self.y)
638
639    @unittest.skip("not yet supported")
640    def test_slit_high(self):
641        """
642        Slit smearing with height 0.005
643        """
644        resolution = Slit1D(self.x, qx_width=0, qy_width=0.005, q_calc=self.x)
645        theory = self.Iq(resolution.q_calc)
646        output = resolution.apply(theory)
647        answer = [
648            9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
649            5.5555, 4.5584, 3.5606, 2.5623, 2.0000,
650            ]
651        np.testing.assert_allclose(output, answer, atol=1e-4)
652
653    @unittest.skip("not yet supported")
654    def test_slit_both_high(self):
655        """
656        Slit smearing with width < 100*height.
657        """
658        q = np.logspace(-4, -1, 10)
659        resolution = Slit1D(q, qx_width=0.2, qy_width=np.inf)
660        theory = 1000*self.Iq(resolution.q_calc**4)
661        output = resolution.apply(theory)
662        answer = [
663            8.85785, 8.43012, 7.92687, 6.94566, 6.03660,
664            5.40363, 4.40655, 3.40880, 2.41058, 2.00000,
665            ]
666        np.testing.assert_allclose(output, answer, atol=1e-4)
667
668    @unittest.skip("not yet supported")
669    def test_slit_wide(self):
670        """
671        Slit smearing with width 0.0002
672        """
673        resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0, q_calc=self.x)
674        theory = self.Iq(resolution.q_calc)
675        output = resolution.apply(theory)
676        answer = [
677            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
678            ]
679        np.testing.assert_allclose(output, answer, atol=1e-4)
680
681    @unittest.skip("not yet supported")
682    def test_slit_both_wide(self):
683        """
684        Slit smearing with width > 100*height.
685        """
686        resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0.000001,
687                            q_calc=self.x)
688        theory = self.Iq(resolution.q_calc)
689        output = resolution.apply(theory)
690        answer = [
691            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
692            ]
693        np.testing.assert_allclose(output, answer, atol=1e-4)
694
695    def test_pinhole_zero(self):
696        """
697        Pinhole smearing with perfect resolution
698        """
699        resolution = Pinhole1D(self.x, 0.0*self.x)
700        theory = self.Iq(resolution.q_calc)
701        output = resolution.apply(theory)
702        np.testing.assert_equal(output, self.y)
703
704    # TODO: turn pinhole/slit demos into tests
705
706    def test_pinhole(self):
707        """
708        Pinhole smearing with dQ = 0.001 [Note: not dQ/Q = 0.001]
709        """
710        resolution = Pinhole1D(self.x, 0.001*np.ones_like(self.x),
711                               q_calc=self.x)
712        theory = 12.0-1000.0*resolution.q_calc
713        output = resolution.apply(theory)
714        # Note: answer came from output of previous run.  Non-integer
715        # values at ends come from the fact that q_calc does not
716        # extend beyond q, and so the weights don't balance.
717        answer = [
718            10.47037734, 9.86925860,
719            9., 8., 7., 6., 5., 4.,
720            3.13074140, 2.52962266,
721            ]
722        np.testing.assert_allclose(output, answer, atol=1e-8)
723
724
725class IgorComparisonTest(unittest.TestCase):
726    """
727    Test resolution calculations against those returned by Igor.
728    """
729
730    def setUp(self):
731        self.pars = TEST_PARS_PINHOLE_SPHERE
732        from sasmodels import core
733        self.model = core.load_model("sphere", dtype='double')
734
735    def _eval_sphere(self, pars, resolution):
736        from sasmodels import direct_model
737        kernel = self.model.make_kernel([resolution.q_calc])
738        theory = direct_model.call_kernel(kernel, pars)
739        result = resolution.apply(theory)
740        kernel.release()
741        return result
742
743    def _compare(self, q, output, answer, tolerance):
744        #err = (output - answer)/answer
745        #idx = abs(err) >= tolerance
746        #problem = zip(q[idx], output[idx], answer[idx], err[idx])
747        #print("\n".join(str(v) for v in problem))
748        np.testing.assert_allclose(output, answer, rtol=tolerance)
749
750    def test_perfect(self):
751        """
752        Compare sphere model with NIST Igor SANS, no resolution smearing.
753        """
754        pars = TEST_PARS_SLIT_SPHERE
755        data_string = TEST_DATA_SLIT_SPHERE
756
757        data = np.loadtxt(data_string.split('\n')).T
758        q, _, answer, _ = data
759        resolution = Perfect1D(q)
760        output = self._eval_sphere(pars, resolution)
761        self._compare(q, output, answer, 1e-6)
762
763    @unittest.skip("suppress comparison with old version; pinhole calc changed")
764    def test_pinhole(self):
765        """
766        Compare pinhole resolution smearing with NIST Igor SANS
767        """
768        pars = TEST_PARS_PINHOLE_SPHERE
769        data_string = TEST_DATA_PINHOLE_SPHERE
770
771        data = np.loadtxt(data_string.split('\n')).T
772        q, q_width, answer = data
773        resolution = Pinhole1D(q, q_width)
774        output = self._eval_sphere(pars, resolution)
775        # TODO: relative error should be lower
776        self._compare(q, output, answer, 3e-4)
777
778    @unittest.skip("suppress comparison with old version; pinhole calc changed")
779    def test_pinhole_romberg(self):
780        """
781        Compare pinhole resolution smearing with romberg integration result.
782        """
783        pars = TEST_PARS_PINHOLE_SPHERE
784        data_string = TEST_DATA_PINHOLE_SPHERE
785        pars['radius'] *= 5
786
787        data = np.loadtxt(data_string.split('\n')).T
788        q, q_width, answer = data
789        answer = romberg_pinhole_1d(q, q_width, self.model, pars)
790        ## Getting 0.1% requires 5 sigma and 200 points per fringe
791        #q_calc = interpolate(pinhole_extend_q(q, q_width, nsigma=5),
792        #                     2*np.pi/pars['radius']/200)
793        #tol = 0.001
794        ## The default 2.5 sigma and no extra points gets 1%
795        q_calc = None  # type: np.ndarray
796        tol = 0.01
797        resolution = Pinhole1D(q, q_width, q_calc=q_calc)
798        output = self._eval_sphere(pars, resolution)
799        if 0: # debug plot
800            import matplotlib.pyplot as plt  # type: ignore
801            resolution = Perfect1D(q)
802            source = self._eval_sphere(pars, resolution)
803            plt.loglog(q, source, '.')
804            plt.loglog(q, answer, '-', hold=True)
805            plt.loglog(q, output, '-', hold=True)
806            plt.show()
807        self._compare(q, output, answer, tol)
808
809    def test_slit(self):
810        """
811        Compare slit resolution smearing with NIST Igor SANS
812        """
813        pars = TEST_PARS_SLIT_SPHERE
814        data_string = TEST_DATA_SLIT_SPHERE
815
816        data = np.loadtxt(data_string.split('\n')).T
817        q, delta_qv, _, answer = data
818        resolution = Slit1D(q, qx_width=delta_qv, qy_width=0)
819        output = self._eval_sphere(pars, resolution)
820        # TODO: eliminate Igor test since it is too inaccurate to be useful.
821        # This means we can eliminate the test data as well, and instead
822        # use a generated q vector.
823        self._compare(q, output, answer, 0.5)
824
825    def test_slit_romberg(self):
826        """
827        Compare slit resolution smearing with romberg integration result.
828        """
829        pars = TEST_PARS_SLIT_SPHERE
830        data_string = TEST_DATA_SLIT_SPHERE
831
832        data = np.loadtxt(data_string.split('\n')).T
833        q, delta_qv, _, answer = data
834        answer = romberg_slit_1d(q, delta_qv, 0., self.model, pars)
835        q_calc = slit_extend_q(interpolate(q, 2*np.pi/pars['radius']/20),
836                               delta_qv[0], 0.)
837        resolution = Slit1D(q, qx_width=delta_qv, qy_width=0, q_calc=q_calc)
838        output = self._eval_sphere(pars, resolution)
839        # TODO: relative error should be lower
840        self._compare(q, output, answer, 0.025)
841
842    def test_ellipsoid(self):
843        """
844        Compare romberg integration for ellipsoid model.
845        """
846        from .core import load_model
847        pars = {
848            'scale':0.05,
849            'radius_polar':500, 'radius_equatorial':15000,
850            'sld':6, 'sld_solvent': 1,
851            }
852        form = load_model('ellipsoid', dtype='double')
853        q = np.logspace(log10(4e-5), log10(2.5e-2), 68)
854        width, height = 0.117, 0.
855        resolution = Slit1D(q, qx_width=width, qy_width=height)
856        answer = romberg_slit_1d(q, width, height, form, pars)
857        output = resolution.apply(eval_form(resolution.q_calc, form, pars))
858        # TODO: 10% is too much error; use better algorithm
859        #print(np.max(abs(answer-output)/answer))
860        self._compare(q, output, answer, 0.1)
861
862    #TODO: can sas q spacing be too sparse for the resolution calculation?
863    @unittest.skip("suppress sparse data test; not supported by current code")
864    def test_pinhole_sparse(self):
865        """
866        Compare pinhole resolution smearing with NIST Igor SANS on sparse data
867        """
868        pars = TEST_PARS_PINHOLE_SPHERE
869        data_string = TEST_DATA_PINHOLE_SPHERE
870
871        data = np.loadtxt(data_string.split('\n')).T
872        q, q_width, answer = data[:, ::20] # Take every nth point
873        resolution = Pinhole1D(q, q_width)
874        output = self._eval_sphere(pars, resolution)
875        self._compare(q, output, answer, 1e-6)
876
877
878# pinhole sphere parameters
879TEST_PARS_PINHOLE_SPHERE = {
880    'scale': 1.0, 'background': 0.01,
881    'radius': 60.0, 'sld': 1, 'sld_solvent': 6.3,
882    }
883# Q, dQ, I(Q) calculated by NIST Igor SANS package
884TEST_DATA_PINHOLE_SPHERE = """\
8850.001278 0.0002847 2538.41176383
8860.001562 0.0002905 2536.91820405
8870.001846 0.0002956 2535.13182479
8880.002130 0.0003017 2533.06217813
8890.002414 0.0003087 2530.70378586
8900.002698 0.0003165 2528.05024192
8910.002982 0.0003249 2525.10408349
8920.003266 0.0003340 2521.86667499
8930.003550 0.0003437 2518.33907750
8940.003834 0.0003539 2514.52246995
8950.004118 0.0003646 2510.41798319
8960.004402 0.0003757 2506.02690988
8970.004686 0.0003872 2501.35067884
8980.004970 0.0003990 2496.38678318
8990.005253 0.0004112 2491.16237596
9000.005537 0.0004237 2485.63911673
9010.005821 0.0004365 2479.83657083
9020.006105 0.0004495 2473.75676948
9030.006389 0.0004628 2467.40145990
9040.006673 0.0004762 2460.77293372
9050.006957 0.0004899 2453.86724627
9060.007241 0.0005037 2446.69623838
9070.007525 0.0005177 2439.25775219
9080.007809 0.0005318 2431.55421398
9090.008093 0.0005461 2423.58785521
9100.008377 0.0005605 2415.36158137
9110.008661 0.0005750 2406.87009473
9120.008945 0.0005896 2398.12841186
9130.009229 0.0006044 2389.13360806
9140.009513 0.0006192 2379.88958042
9150.009797 0.0006341 2370.39776774
9160.010080 0.0006491 2360.69528793
9170.010360 0.0006641 2350.85169027
9180.010650 0.0006793 2340.42023633
9190.010930 0.0006945 2330.11206013
9200.011220 0.0007097 2319.20109972
9210.011500 0.0007251 2308.43503981
9220.011780 0.0007404 2297.44820179
9230.012070 0.0007558 2285.83853677
9240.012350 0.0007713 2274.41290746
9250.012640 0.0007868 2262.36219581
9260.012920 0.0008024 2250.51169731
9270.013200 0.0008180 2238.45596231
9280.013490 0.0008336 2225.76495666
9290.013770 0.0008493 2213.29618391
9300.014060 0.0008650 2200.19110751
9310.014340 0.0008807 2187.34050325
9320.014620 0.0008965 2174.30529864
9330.014910 0.0009123 2160.61632548
9340.015190 0.0009281 2147.21038112
9350.015470 0.0009440 2133.62023580
9360.015760 0.0009598 2119.37907426
9370.016040 0.0009757 2105.45234903
9380.016330 0.0009916 2090.86319102
9390.016610 0.0010080 2076.60576032
9400.016890 0.0010240 2062.19214565
9410.017180 0.0010390 2047.10550219
9420.017460 0.0010550 2032.38715621
9430.017740 0.0010710 2017.52560123
9440.018030 0.0010880 2001.99124318
9450.018310 0.0011040 1986.84662060
9460.018600 0.0011200 1971.03389745
9470.018880 0.0011360 1955.61395119
9480.019160 0.0011520 1940.08291563
9490.019450 0.0011680 1923.87672225
9500.019730 0.0011840 1908.10656374
9510.020020 0.0012000 1891.66297192
9520.020300 0.0012160 1875.66789021
9530.020580 0.0012320 1859.56357196
9540.020870 0.0012490 1842.79468290
9550.021150 0.0012650 1826.50064489
9560.021430 0.0012810 1810.11533702
9570.021720 0.0012970 1793.06840882
9580.022000 0.0013130 1776.51153580
9590.022280 0.0013290 1759.87201249
9600.022570 0.0013460 1742.57354412
9610.022850 0.0013620 1725.79397319
9620.023140 0.0013780 1708.35831550
9630.023420 0.0013940 1691.45256069
9640.023700 0.0014110 1674.48561783
9650.023990 0.0014270 1656.86525366
9660.024270 0.0014430 1639.79847285
9670.024550 0.0014590 1622.68887088
9680.024840 0.0014760 1604.96421100
9690.025120 0.0014920 1587.85768129
9700.025410 0.0015080 1569.99297335
9710.025690 0.0015240 1552.84580279
9720.025970 0.0015410 1535.54074115
9730.026260 0.0015570 1517.75249337
9740.026540 0.0015730 1500.40115023
9750.026820 0.0015900 1483.03632237
9760.027110 0.0016060 1465.05942429
9770.027390 0.0016220 1447.67682181
9780.027670 0.0016390 1430.46495191
9790.027960 0.0016550 1412.49232282
9800.028240 0.0016710 1395.13182318
9810.028520 0.0016880 1377.93439837
9820.028810 0.0017040 1359.99528971
9830.029090 0.0017200 1342.67274512
9840.029370 0.0017370 1325.55375609
985"""
986
987# Slit sphere parameters
988TEST_PARS_SLIT_SPHERE = {
989    'scale': 0.01, 'background': 0.01,
990    'radius': 60000, 'sld': 1, 'sld_solvent': 4,
991    }
992# Q dQ I(Q) I_smeared(Q)
993TEST_DATA_SLIT_SPHERE = """\
9942.26097e-05 0.117 5.5781372896e+09 1.4626077708e+06
9952.53847e-05 0.117 5.0363141458e+09 1.3117318023e+06
9962.81597e-05 0.117 4.4850108103e+09 1.1594863713e+06
9973.09347e-05 0.117 3.9364658459e+09 1.0094881630e+06
9983.37097e-05 0.117 3.4019975074e+09 8.6518941303e+05
9993.92597e-05 0.117 2.4139519814e+09 6.0232158311e+05
10004.48097e-05 0.117 1.5816877820e+09 3.8739994090e+05
10015.03597e-05 0.117 9.3715407224e+08 2.2745304775e+05
10025.59097e-05 0.117 4.8387917428e+08 1.2101295768e+05
10036.14597e-05 0.117 2.0193586928e+08 6.0055107771e+04
10046.70097e-05 0.117 5.5886110911e+07 3.2749521065e+04
10057.25597e-05 0.117 3.7782348010e+06 2.6350963616e+04
10067.81097e-05 0.117 5.3407817904e+06 2.9624963314e+04
10078.36597e-05 0.117 2.7975485523e+07 3.4403962254e+04
10088.92097e-05 0.117 4.9845448282e+07 3.6130017903e+04
10099.47597e-05 0.117 6.0092588905e+07 3.3495107849e+04
10101.00310e-04 0.117 5.6823430831e+07 2.7475726279e+04
10111.05860e-04 0.117 4.3857024036e+07 2.0144282226e+04
10121.11410e-04 0.117 2.7277144760e+07 1.3647403260e+04
10131.22510e-04 0.117 3.3119334113e+06 6.6519711526e+03
10141.33610e-04 0.117 1.4412859402e+06 6.9726212813e+03
10151.44710e-04 0.117 8.5620162463e+06 8.1441335775e+03
10161.55810e-04 0.117 9.6957429033e+06 6.4559996521e+03
10171.66910e-04 0.117 4.3818341914e+06 3.6252154156e+03
10181.78010e-04 0.117 2.7448997387e+05 2.4006505342e+03
10191.89110e-04 0.117 8.0472009936e+05 2.8187789089e+03
10202.00210e-04 0.117 2.8149552834e+06 3.0915662855e+03
10212.11310e-04 0.117 2.7510907861e+06 2.3722530293e+03
10222.22410e-04 0.117 1.0053133293e+06 1.4473468311e+03
10232.33510e-04 0.117 5.8428305052e+03 1.2048540556e+03
10242.44610e-04 0.117 5.1699305004e+05 1.4625670042e+03
10252.55710e-04 0.117 1.2120227268e+06 1.5010705973e+03
10262.66810e-04 0.117 9.7896842846e+05 1.1336343426e+03
10272.77910e-04 0.117 2.5507264791e+05 8.1848818080e+02
10283.05660e-04 0.117 5.2403101181e+05 7.4913374357e+02
10293.33410e-04 0.117 5.8699343809e+04 4.4669964560e+02
10303.61160e-04 0.117 3.0844327150e+05 4.6774007542e+02
10313.88910e-04 0.117 8.3360142970e+03 2.7169550220e+02
10324.16660e-04 0.117 1.8630080583e+05 3.0710983679e+02
10334.44410e-04 0.117 3.1616804732e-01 1.7959006831e+02
10344.72160e-04 0.117 1.1299016314e+05 2.0763952339e+02
10354.99910e-04 0.117 2.9952522747e+03 1.2536542765e+02
10365.27660e-04 0.117 6.7625695649e+04 1.4013969777e+02
10375.55410e-04 0.117 7.6927460089e+03 8.2145593180e+01
10386.10910e-04 0.117 1.1229057779e+04 8.4519745643e+01
10396.66410e-04 0.117 1.3035567943e+04 8.1554625609e+01
10407.21910e-04 0.117 1.3309931343e+04 7.4437319172e+01
10417.77410e-04 0.117 1.2462626212e+04 6.4697088261e+01
10428.32910e-04 0.117 1.0912927143e+04 5.3773301044e+01
10438.88410e-04 0.117 9.0172597469e+03 4.2843375753e+01
10449.43910e-04 0.117 7.0496495917e+03 3.2771032724e+01
10459.99410e-04 0.117 5.2030483682e+03 2.4113557144e+01
10461.05491e-03 0.117 3.5988976711e+03 1.7160773658e+01
10471.11041e-03 0.117 2.2996060652e+03 1.2016626459e+01
10481.22141e-03 0.117 6.4766590598e+02 6.0373017740e+00
10491.33241e-03 0.117 4.1963483264e+01 4.5215452974e+00
10501.44341e-03 0.117 6.3370708246e+01 5.1054681903e+00
10511.55441e-03 0.117 3.0736750577e+02 5.9176165298e+00
10521.66541e-03 0.117 5.0327682399e+02 5.9815000189e+00
10531.77641e-03 0.117 5.4084331454e+02 5.1634639625e+00
10541.88741e-03 0.117 4.3488671756e+02 3.8535158148e+00
10551.99841e-03 0.117 2.6322287860e+02 2.5824997753e+00
10562.10941e-03 0.117 1.0793633150e+02 1.7315517194e+00
10572.22041e-03 0.117 1.8474448850e+01 1.4077213604e+00
10582.33141e-03 0.117 1.5864062279e+00 1.4771560682e+00
10592.44241e-03 0.117 3.2267213848e+01 1.6916253448e+00
10602.55341e-03 0.117 7.4289116207e+01 1.8274751193e+00
10612.66441e-03 0.117 9.9000521929e+01 1.7706812289e+00
1062"""
1063
1064def main():
1065    """
1066    Run tests given is sys.argv.
1067
1068    Returns 0 if success or 1 if any tests fail.
1069    """
1070    import sys
1071    import xmlrunner  # type: ignore
1072
1073    suite = unittest.TestSuite()
1074    suite.addTest(unittest.defaultTestLoader.loadTestsFromModule(sys.modules[__name__]))
1075
1076    runner = xmlrunner.XMLTestRunner(output='logs')
1077    result = runner.run(suite)
1078    return 1 if result.failures or result.errors else 0
1079
1080
1081############################################################################
1082# usage demo
1083############################################################################
1084
1085def _eval_demo_1d(resolution, title):
1086    import sys
1087    from sasmodels import core
1088    from sasmodels import direct_model
1089    name = sys.argv[1] if len(sys.argv) > 1 else 'cylinder'
1090
1091    if name == 'cylinder':
1092        pars = {'length':210, 'radius':500, 'background': 0}
1093    elif name == 'teubner_strey':
1094        pars = {'a2':0.003, 'c1':-1e4, 'c2':1e10, 'background':0.312643}
1095    elif name == 'sphere' or name == 'spherepy':
1096        pars = TEST_PARS_SLIT_SPHERE
1097    elif name == 'ellipsoid':
1098        pars = {
1099            'scale':0.05, 'background': 0,
1100            'r_polar':500, 'r_equatorial':15000,
1101            'sld':6, 'sld_solvent': 1,
1102            }
1103    else:
1104        pars = {}
1105    model_info = core.load_model_info(name)
1106    model = core.build_model(model_info)
1107
1108    kernel = model.make_kernel([resolution.q_calc])
1109    theory = direct_model.call_kernel(kernel, pars)
1110    Iq = resolution.apply(theory)
1111
1112    if isinstance(resolution, Slit1D):
1113        width, height = resolution.qx_width, resolution.qy_width
1114        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars)
1115    else:
1116        dq = resolution.q_width
1117        Iq_romb = romberg_pinhole_1d(resolution.q, dq, model, pars)
1118
1119    import matplotlib.pyplot as plt  # type: ignore
1120    plt.loglog(resolution.q_calc, theory, label='unsmeared')
1121    plt.loglog(resolution.q, Iq, label='smeared', hold=True)
1122    plt.loglog(resolution.q, Iq_romb, label='romberg smeared', hold=True)
1123    plt.legend()
1124    plt.title(title)
1125    plt.xlabel("Q (1/Ang)")
1126    plt.ylabel("I(Q) (1/cm)")
1127
1128def demo_pinhole_1d():
1129    """
1130    Show example of pinhole smearing.
1131    """
1132    q = np.logspace(-4, np.log10(0.2), 400)
1133    q_width = 0.1*q
1134    resolution = Pinhole1D(q, q_width)
1135    _eval_demo_1d(resolution, title="10% dQ/Q Pinhole Resolution")
1136
1137def demo_slit_1d():
1138    """
1139    Show example of slit smearing.
1140    """
1141    q = np.logspace(-4, np.log10(0.2), 100)
1142    w = h = 0.
1143    #w = 0.000000277790
1144    w = 0.0277790
1145    #h = 0.00277790
1146    #h = 0.0277790
1147    resolution = Slit1D(q, w, h)
1148    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h))
1149
1150def demo():
1151    """
1152    Run the resolution demos.
1153    """
1154    import matplotlib.pyplot as plt  # type: ignore
1155    plt.subplot(121)
1156    demo_pinhole_1d()
1157    #plt.yscale('linear')
1158    plt.subplot(122)
1159    demo_slit_1d()
1160    #plt.yscale('linear')
1161    plt.show()
1162
1163
1164if __name__ == "__main__":
1165    #demo()
1166    main()
Note: See TracBrowser for help on using the repository browser.