26. Multiple scattering support¶
sasmodels.multiscat
¶
Multiple scattering calculator
Calculate multiple scattering using 2D FFT convolution.
Usage:
python -m sasmodels.multiscat [options] model_name model_par=value ...
Options include:
-h, --help: show help and exit
-n, --nq: the number of mesh points (dq = qmax*window/nq)
-o, --outfile: save results to outfile.txt and outfile_powers.txt
-p, --probability: the scattering probability (0.1)
-q, --qmax: that max q that you care about (0.5)
-r, --random: generate a random parameter set
-s, --seed: generate a random parameter set with a given seed
-w, --window: the extension window (q is calculated for qmax*window)
-2, --2d: perform the calculation for an oriented pattern
Assume the probability of scattering is \(p\). After each scattering event, \(1-p\) neutrons will leave the system and go to the detector, and the remaining \(p\) will scatter again.
Let the scattering probability for \(n\) scattering event at \(q\) be \(f_n(q)\), with
for \(I_1(q)\), the single scattering from the system. After two scattering events, the scattering probability will be the convolution of the first scattering and itself, or \(f_2(q) = (f_1*f_1)(q)\). After \(n\) events it will be \(f_n(q) = (f_1 * \cdots * f_1)(q)\). The total scattering is calculated as the weighted sum of \(f_k\), with weights following the Poisson distribution
for \(\lambda\) determined by the total thickness divided by the mean free path between scattering, giving
The \(k=0\) term is ignored since it involves no scattering. We cut the series when cumulative probability is less than cutoff \(C=99\%\), which is \(\max n\) such that
Using the convolution theorem, where \(F = \mathcal{F}(f)\) is the Fourier transform,
so
Since the Fourier transform is a linear operator, we can move the polynomial expression for the convolution into the transform, giving
In the dilute limit \(\lambda \rightarrow 0\) only the \(k=1\) term is active, and so
therefore we compute
For speed we may use the fast fourier transform with a power of two. The resulting \(I(q)\) will be linearly spaced and likely heavily oversampled. The usual pinhole or slit resolution calculation can performed from these calculated values.
By default single precision OpenCL is used for the calculation. Set the environment variable SAS_OPENCL=none to use double precision numpy FFT instead. The OpenCL versions is about 10x faster on an elderly Mac with Intel HD 4000 graphics. The single precision numerical artifacts don’t seem to seriously impact overall accuracy, though they look pretty bad.
- sasmodels.multiscat.Calculator¶
alias of
sasmodels.multiscat.NumpyCalculator
- class sasmodels.multiscat.ICalculator¶
Bases:
object
Multiple scattering calculator
- fft(Iq)¶
Compute the forward FFT for an image, real -> complex.
- ifft(fourier_frame)¶
Compute the inverse FFT for an image, complex -> complex.
- multiple_scattering(Iq, p, coverage=0.99)¶
Compute multiple scattering for I(q) given scattering probability p.
Given a probability p of scattering with the thickness, the expected number of scattering events, \(\lambda = -\log(1 - p)\), giving a Poisson weighted sum of single, double, triple, etc. scattering patterns. The number of patterns used is based on coverage (default 99%).
- class sasmodels.multiscat.MultipleScattering(qmin=None, qmax=None, nq=None, window=2, probability=None, coverage=0.99, is2d=False, resolution=None, dtype=dtype('float64'))¶
Bases:
sasmodels.resolution.Resolution
Compute multiple scattering using Fourier convolution.
The fourier steps are determined by qmax, the maximum \(q\) value desired, nq the number of \(q\) steps and window, the amount of padding around the circular convolution. The \(q\) spacing will be \(\Delta q = 2 q_\mathrm{max} w / n_q\). If nq is not given it will use \(n_q = 2^k\) such that \(\Delta q < q_\mathrm{min}\).
probability is related to the expected number of scattering events in the sample \(\lambda\) as \(p = 1 - e^{-\lambda}\).
coverage determines how many scattering steps to consider. The default is 0.99, which sets \(n\) such that \(1 \ldots n\) covers 99% of the Poisson probability mass function.
is2d is True then 2D scattering is used, otherwise it accepts and returns 1D scattering.
resolution is the resolution function to apply after multiple scattering. If present, then the resolution \(q\) vectors will provide default values for qmin, qmax and nq.
- apply(theory)¶
Smear theory by the resolution function, returning Iq.
- plot_and_save_powers(theory, result, plot=True, outfile='', background=0.0)¶
- radial_profile(Iqxy)¶
Compute that radial profile for the given Iqxy grid. The grid should be defined as for
- q: numpy.ndarray = None¶
- q_calc: numpy.ndarray = None¶
- class sasmodels.multiscat.NumpyCalculator(dims=None, dtype=dtype('float64'))¶
Bases:
sasmodels.multiscat.ICalculator
Multiple scattering calculator using numpy fft.
- fft(Iq)¶
Compute the forward FFT for an image, real -> complex.
- ifft(fourier_frame)¶
Compute the inverse FFT for an image, complex -> complex.
- multiple_scattering(Iq, p, coverage=0.99)¶
Compute multiple scattering for I(q) given scattering probability p.
Given a probability p of scattering with the thickness, the expected number of scattering events, \(\lambda = -\log(1 - p)\), giving a Poisson weighted sum of single, double, triple, etc. scattering patterns. The number of patterns used is based on coverage (default 99%).
- class sasmodels.multiscat.OpenclCalculator(dims, dtype=dtype('float64'))¶
Bases:
sasmodels.multiscat.ICalculator
Multiple scattering calculator using OpenCL via pyfft.
- fft(Iq)¶
Compute the forward FFT for an image, real -> complex.
- ifft(fourier_frame)¶
Compute the inverse FFT for an image, complex -> complex.
- multiple_scattering(Iq, p, coverage=0.99)¶
Compute multiple scattering for I(q) given scattering probability p.
Given a probability p of scattering with the thickness, the expected number of scattering events, \(\lambda = -\log(1 - p)\), giving a Poisson weighted sum of single, double, triple, etc. scattering patterns. The number of patterns used is based on coverage (default 99%).
- polyval1d = None¶
- polyval1f = None¶
- sasmodels.multiscat.annular_average(qxy, Iqxy, qbins)¶
Compute annular average of points in Iqxy at qbins. The \(q_x\), \(q_y\) coordinates for Iqxy are given in qxy.
- sasmodels.multiscat.main()¶
Command line interface to multiple scattering calculator.
- sasmodels.multiscat.parse_pars(model: sasmodels.modelinfo.ModelInfo, opts: argparse.Namespace) Dict[str, float] ¶
Parse par=val arguments from the command line.
- sasmodels.multiscat.plotxy(q, Iq)¶
Plot q, Iq or (qx, qy), Iqxy.
- sasmodels.multiscat.rebin(x, I, xo)¶
Rebin from edges x, bins I into edges xo.
x and xo should be monotonically increasing.
If x has duplicate values, then all corresponding values at I(x) will be effectively summed into the same bin. If xo has duplicate values, the first bin will contain the entire contents and subsequent bins will contain zeros.
- sasmodels.multiscat.scattering_coeffs(p, coverage=0.99)¶
Return the coefficients of the scattering powers for transmission probability p. This is just the corresponding values for the Poisson distribution for \(\lambda = -\ln(1-p)\) such that \(\sum_{k = 0 \ldots n} P(k; \lambda)\) is larger than coverage.
- sasmodels.multiscat.scattering_powers(Iq, n, dtype='f', transform=None)¶
Calculate the scattering powers up to n.
This includes k=1 even though it should just be Iq itself
The frames are unweighted; to weight scale by \(\lambda^k e^{-\lambda}/k!\).
- sasmodels.multiscat.truncated_poisson_invcdf(coverage, L)¶
Return smallest k such that cdf(k; L) > coverage from the truncated Poisson probability excluding k=0