Crate scatrs 

Fast, linearized, rigorously validated single-scattering calculations for spherical and spheroidal particles.

scatrs computes the optical properties of homogeneous spheres by Lorenz–Mie theory, of layered spheres, and of randomly oriented spheroids by the T-matrix method, and the averages of spheres over size distributions and external mixtures, with closed-form derivatives and a declared, tested accuracy. The same crate is the core of the Python package scatrs. This reference covers the Rust API; the scatrs documentation also has a user guide with examples, the theory behind every method, and the validation of the results.

Module	Contents
mie	One sphere: Mie coefficients, efficiencies, amplitudes, phase matrix and Greek coefficients, with their derivatives.
layered	One layered sphere, a core and concentric shells: the same outputs, with derivatives with respect to every layer.
tmatrix	One randomly oriented spheroid by the T-matrix method: efficiencies, Greek coefficients and six-element phase matrix, the T-matrix itself and its convergence, with derivatives.
ensemble	Averages of spheres and layered spheres over size distributions, external mixtures and batches of distributions, with error estimates and derivatives.
special	Gauss–Legendre quadrature and Wigner d-functions.
scalar	The Real trait the numerical kernels are generic over, and its double-double and dual-number types.
error	The Error that every fallible function returns.

One sphere

A Mie workspace solves one sphere at a time; everything else is read from the solution.

use num_complex::Complex;
use scatrs::mie::{Mie, MieOptions, Normalization, PhaseBuffers};

let mut mie = Mie::<f64>::new();
let sphere = mie.compute(10.0, Complex::new(1.33, 1e-8), &MieOptions::new())?;

// Efficiencies (cross sections divided by πr²) and the asymmetry parameter.
let q = sphere.efficiencies();
println!("Q_ext = {}, Q_sca = {}, g = {}", q.q_ext, q.q_sca, q.g);

// The normalized phase matrix at cosines of the scattering angle.
let mu = [1.0, 0.5, 0.0, -0.5, -1.0];
let mut f = [[0.0; 5]; 4];
let [f11, f12, f33, f34] = &mut f;
sphere.phase_matrix(&mu, Normalization::Normalized, PhaseBuffers { f11, f12, f33, f34 })?;

// Its Greek coefficients, every nonzero order: α₁⁰ = 1 and α₁¹ = 3g.
let greek = sphere.greek(None, Normalization::Normalized)?;
assert_eq!(greek.alpha1[0], 1.0);
assert!((greek.alpha1[1] - 3.0 * q.g).abs() < 1e-12);

compute_with_derivatives also gives the derivatives of every output with respect to the size parameter and the real and imaginary parts of the refractive index, in closed form, for 1.2 to 2.2 times the cost of the values:

let mut mie = Mie::<f64>::new();
let sphere = mie.compute_with_derivatives(10.0, Complex::new(1.5, 0.01), &MieOptions::new())?;
let jacobian = sphere.efficiencies_jacobian().expect("computed with derivatives");
println!("dQ_ext/dx = {}, dQ_ext/d(Im m) = {}", jacobian.x.q_ext, jacobian.m_im.q_ext);

A layered sphere

A Layered workspace solves a core in concentric shells, given the size parameters of the interfaces from the core outward and the layers’ refractive indices; the solution has the outputs of a sphere’s.

use num_complex::Complex;
use scatrs::layered::Layered;
use scatrs::mie::MieOptions;

// Black carbon coated in sulfate.
let x = [0.7, 1.3];
let m = [Complex::new(1.95, 0.79), Complex::new(1.43, 1e-8)];
let mut layered = Layered::<f64>::new();
let particle = layered.compute(&x, &m, &MieOptions::new())?;
println!("Q_abs = {}", particle.efficiencies().q_abs);

A spheroid

A TMatrix workspace solves one spheroid, of aspect ratio a/c (oblate above 1) and equal-volume size parameter x, converging its T-matrix to a tolerance; the solution has the orientation-averaged outputs, with F₂₂ and F₄₄ independent of F₁₁ and F₃₃, and an error estimate.

use num_complex::Complex;
use scatrs::mie::Normalization;
use scatrs::tmatrix::{PhaseBuffers, Shape, TMatrix, TMatrixOptions};

let shape = Shape::Spheroid { aspect_ratio: 2.0 };
let mut tm = TMatrix::new();
let s = tm.compute(shape, 5.0, Complex::new(1.53, 0.008), &TMatrixOptions::new())?;
println!("Q_ext = {}, estimated error {:e}", s.efficiencies().q_ext, s.error_estimate());

// Backscattering depolarizes: F₂₂ < F₁₁ at 180°, where a sphere has F₂₂ = F₁₁.
let mu = [-1.0];
let mut f = [[0.0; 1]; 6];
let [f11, f12, f22, f33, f34, f44] = &mut f;
s.phase_matrix(&mu, Normalization::Normalized, PhaseBuffers { f11, f12, f22, f33, f34, f44 })?;
assert!(f[2][0] < f[0][0]);

A size distribution

ensemble::compute averages over a size distribution at each wavelength, in parallel. Radii and wavelengths share any length unit, and cross sections are per particle, in that unit squared. Every output has an error estimate, and converged says whether the size integration met its tolerance.

use num_complex::Complex;
use scatrs::ensemble::{
    EnsembleOptions, EnsembleParam, ParamKind, SizeDistribution, compute,
};

let sulfate = SizeDistribution::LogNormal { median_radius: 0.08, width: 1.6 }; // µm
let wavelengths = [0.35, 0.55, 1.0]; // µm
let m = [Complex::new(1.43, 1e-8); 3]; // one refractive index per wavelength
let median_radius = EnsembleParam::new(0, None, ParamKind::Distribution("median_radius"));
let options = EnsembleOptions::new().with_derivatives(vec![median_radius]);
for r in compute(&sulfate, &wavelengths, &m, &options)? {
    assert!(r.converged);
    let dc_ext = r.derivatives[0].c_ext; // the derivatives follow the order of r.params
    println!(
        "C_ext = {} ± {} µm², ω = {}, dC_ext/dr_g = {dc_ext} µm",
        r.values.c_ext, r.errors.c_ext, r.values.ssa,
    );
}

Conventions

• Refractive index. m = n + ik, with k ≥ 0 for absorbing particles (Bohren & Huffman, time dependence e^(−iωt)), relative to the surrounding medium. A negative imaginary part is an error, never silently flipped. Codes that write m = n − ik, such as MIEV0 and sasktran2, need the complex conjugate, and report the complex conjugates of the coefficients and amplitudes.
• Size. mie takes the size parameter x = 2πr/λ and returns efficiencies, layered that of each interface (the outer one normalizing), and tmatrix that of the sphere of equal volume; ensemble takes radii and wavelengths and returns cross sections.
• Angles. Functions of angle take μ = cos Θ, with Θ the scattering angle.
• Phase matrix. Normalized, ½∫F₁₁ dμ = 1; the degree of linear polarization is −F₁₂/F₁₁.
• Greek coefficients. The expansion of Mishchenko, Travis & Lacis (2002) in Wigner d-functions; normalized, α₁⁰ = 1 and α₁¹ = 3g.

Accuracy

Every Mie output has a declared accuracy inside the validated domain, 10⁻⁶ ≤ x ≤ 2·10⁴, 0 < Re m ≤ 10, 0 ≤ Im m ≤ 10, |m| ≥ 0.1 (see VALIDATED_DOMAIN), checked in the tests against independent arbitrary-precision reference values. Outside it, calls return Error::OutsideValidatedDomain unless allow_unvalidated is set. Layered spheres and spheroids have validated domains of their own (layered::VALIDATED_DOMAIN, tmatrix::VALIDATED_DOMAIN); a spheroid’s error is relative to its error estimate as well as to its tolerance. The ensembles estimate the error of their size integration instead.

Precision and derivatives through generic code

The numerical kernels are generic over Real: besides f64, they run in double-double precision (DoubleDouble, about 32 significant digits) and with forward-mode derivatives (Dual). The ensembles are f64 only.

use num_complex::Complex;
use scatrs::mie::{Mie, MieOptions};
use scatrs::scalar::DoubleDouble;

let dd = DoubleDouble::from_f64;
let mut mie = Mie::<DoubleDouble>::new();
let sphere = mie.compute(dd(10.0), Complex::new(dd(1.33), dd(1e-8)), &MieOptions::new())?;
println!("Q_ext = {:?}", sphere.efficiencies().q_ext);

Threads and determinism

A Mie workspace serves one thread; make one per thread. The ensembles run on rayon’s global thread pool, or on a pool the caller installs, and give bitwise identical results for any number of threads.

Re-exports

pub use error::Error;

pub use error::Result;

pub use scalar::Real;

Modules

ensemble

Ensembles of particles: size distributions, adaptive size integration and external mixtures.

error

Error type shared by every scatrs computation.

layered

Scattering by layered spheres: a core and any number of concentric shells.

mie

Lorenz–Mie scattering by homogeneous spheres.

scalar

Scalar types that every numerical kernel is generic over.

special

Special functions and quadrature rules shared by all scattering methods.

tmatrix

Scattering by randomly oriented spheroids, by the T-matrix method.