# Home

## 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

\[f_1(q) = \frac{I_1(q)}{\int I_1(q) {\rm d}q}\]

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

\[P(k; \lambda) = \frac{\lambda^k e^{-\lambda}}{k!}\]

for \(\lambda\) determined by the total thickness divided by the mean free path between scattering, giving

\[I(q) = \sum_{k=0}^\infty P(k; \lambda) f_k(q)\]

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

\[\sum_{k=1}^n \frac{P(k; \lambda)}{1 - P(0; \lambda)} < C\]

Using the convolution theorem, where \(F = \mathcal{F}(f)\) is the Fourier transform,

\[f * g = \mathcal{F}^{-1}\{\mathcal{F}\{f\} \cdot \mathcal{F}\{g\}\}\]

so

\[f * \ldots * f = \mathcal{F}^{-1}\{ F^n \}\]

Since the Fourier transform is a linear operator, we can move the polynomial expression for the convolution into the transform, giving

\[I(q) = \mathcal{F}^{-1}\left\{ \sum_{k=1}^{n} P(k; \lambda) F^k \right\}\]

In the dilute limit \(\lambda \rightarrow 0\) only the \(k=1\) term is active, and so

\[P(1; \lambda) = \lambda e^{-\lambda} = \int I_1(q) {\rm d}q\]

therefore we compute

\[I(q) = \mathcal{F}^{-1}\left\{ \sum_{l=1}^{n} \frac{P(k; \lambda)}{P(1; \lambda))} F^k \right\} \, \int I_1(q) {\rm d}q = \mathcal{F}^{-1}\left\{ \sum_{l=1}^{n} \frac{\lambda^{k-1}}{k!} F^k \right\} \, \int I_1(q) {\rm d}q\]

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 `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: `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: ndarray = None
  q_calc: ndarray = None

class sasmodels.multiscat.NumpyCalculator(dims=None, dtype=dtype('float64'))

  Bases: `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: `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: ModelInfo, opts: 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