source: sasmodels/sasmodels/resolution.py @ df0d2ca

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

update pinhole resolution tests

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