API reference

PackLab keeps compiled extensions alongside its public Python APIs. The reference below is generated from the installed package, so signatures and docstrings remain aligned with the released interfaces.

Start with one of these public namespaces:

• PackLab.samplers for radius distributions;

• PackLab.monte_carlo for random sequential adsorption (RSA);

• PackLab.analytical for Percus–Yevick mixture calculations;

• PackLab.scattering for the optional PyMieSim integration.

Only public classes and functions are listed below. Use the module-source link next to each object when implementation details are useful.

Radius samplers

Radius sampler classes for RSA simulations with optional radius binning.

class ConstantRadiusSampler

Bases: RadiusSampler

Radius sampler that always returns one radius.

Parameters:

• radius (pint.Quantity) – Particle radius.

• bins (int, default=0) – Optional number of radius bins.

class DiscreteRadiusSampler

Bases: RadiusSampler

Finite set of particle radii with explicit probability weights.

Parameters:

• radii (pint.Quantity, shape (n,)) – Discrete radius values.

• weights (array-like, shape (n,)) – Non-negative relative probabilities. They are normalized internally.

class LogNormalRadiusSampler

Bases: RadiusSampler

Log-normal distribution of strictly positive particle radii.

Parameters:

• median_radius (pint.Quantity) – Median radius.

• geometric_standard_deviation (float) – Multiplicative standard deviation; must exceed one.

• maximum_radius_clip (pint.Quantity) – Upper radius cutoff.

• bins (int, default=0) – Optional number of radius bins.

class NormalRadiusSampler

Bases: RadiusSampler

Clipped normal distribution of particle radii.

Parameters:

• mean (pint.Quantity) – Parameters of the normal distribution.

• standard_deviation (pint.Quantity) – Parameters of the normal distribution.

• maximum_clip (pint.Quantity, optional) – Upper radius cutoff. Defaults to mean + 5 * standard_deviation.

• bins (int, default=0) – Optional number of radius bins.

class RadiusSampler

Bases: pybind11_object

Base class for particle-radius distributions.

Notes

Concrete samplers generate radii in SI units internally. Their Python constructors accept Pint quantities and to_bins returns radii in meters.

property number_of_bins

Number of bins; zero means that binning is disabled.

Type:

int

plot_histogram(self: PackLab.samplers.RadiusSampler, ax: object = None, unit: object = 'nanometer', **kwargs) → object

Plot the binned radius distribution as a probability-mass histogram.

The bars show the representative radius classes and normalized number fractions returned by to_bins(). Continuous samplers therefore require bins to be positive.

Parameters:

• ax (matplotlib.axes.Axes, optional) – Axes receiving the bars. The current axes are used when omitted.

• unit (str or pint.Unit, default="nanometer") – Unit used on the radius axis.

• **kwargs – Additional keyword arguments forwarded to matplotlib.axes.Axes.bar().

Returns:

The axes containing the histogram.

Return type:

matplotlib.axes.Axes

set_number_of_bins(self: PackLab.samplers.RadiusSampler, bins: SupportsInt | SupportsIndex) → None

Set the number of discrete radius bins.

Parameters:

bins (int) – Number of bins. A value of zero disables binning.

to_bins(self: PackLab.samplers.RadiusSampler) → tuple[object, numpy.typing.NDArray[numpy.float64]]

Convert the distribution to representative discrete radius classes.

Returns:

• radii (pint.Quantity, shape (n_bins,)) – Representative radii in meters.

• weights (numpy.ndarray, shape (n_bins,))

• Normalized probability weights.

class UniformRadiusSampler

Bases: RadiusSampler

Uniform distribution of particle radii.

Parameters:

• minimum_radius (pint.Quantity) – Inclusive bounds of the radius distribution.

• maximum_radius (pint.Quantity) – Inclusive bounds of the radius distribution.

• bins (int, default=0) – Optional number of radius bins.

Monte-Carlo hard-sphere workflows

Packing domain

class PackingDomain

Bases: pybind11_object

Three-dimensional rectangular domain for random sequential addition.

Parameters:

• length_x (pint.Quantity) – Positive box lengths, convertible to meters.

• length_y (pint.Quantity) – Positive box lengths, convertible to meters.

• length_z (pint.Quantity) – Positive box lengths, convertible to meters.

• use_periodic_boundaries (bool) – Whether to apply periodic boundary conditions on all three axes.

property length_x

Length of the x axis in meters.

Type:

pint.Quantity

property length_y

Length of the y axis in meters.

Type:

pint.Quantity

property length_z

Length of the z axis in meters.

Type:

pint.Quantity

scale(self: PackLab.monte_carlo.domain.PackingDomain, scale_factor: SupportsFloat | SupportsIndex) → None

Scale all box lengths in place.

Parameters:

scale_factor (float) – Positive multiplicative scale factor.

property use_periodic_boundaries

Whether periodic boundary conditions are enabled.

Type:

bool

property volume

Box volume in cubic meters.

Type:

float

RSA simulation

Random sequential addition of non overlapping spheres in a 3D box

class PackingConfiguration

Bases: pybind11_object

Sphere centers, radii, and size-class labels of a packing.

property classes_index

Integer class index for each sphere

property count

Number of spheres in the configuration.

property number_of_classes

Number of distinct particle radius classes

property positions

List of sphere center positions

property radii

List of sphere radii

total_sphere_volume(self: PackLab.monte_carlo.simulator.PackingConfiguration) → object

Compute the total volume occupied by the spheres.

class RSAOptions

Bases: pybind11_object

Stopping criteria and numerical settings for an RSA simulation.

maximum_attempts

Total trial-insertion limit.

Type:

int

maximum_spheres

Sphere-count limit; zero disables this criterion.

Type:

int

target_packing_fraction

Target volume fraction; zero disables this criterion.

Type:

float

class RSASimulator

Bases: pybind11_object

Random Sequential Addition simulator for non-overlapping spheres.

Parameters:

• domain (PackingDomain) – Spatial domain and boundary conditions.

• radius_sampler (RadiusSampler) – Distribution used for candidate sphere radii.

• options (RSAOptions) – Stopping criteria and numerical settings.

Notes

Call run() to generate a PackingResult.

reset(self: PackLab.monte_carlo.simulator.RSASimulator) → None

Reset the simulation to its initial state.

run(self: PackLab.monte_carlo.simulator.RSASimulator) → object

Run the simulation and return a PackingResult.

property sphere_configuration

Current sphere configuration

Metropolis sampling

Equilibrium Metropolis sampling for hard-sphere configurations.

class MetropolisOptions

Bases: pybind11_object

Numerical settings for fixed-volume hard-sphere Metropolis sampling.

The maximum displacement is the maximum absolute displacement proposed along each coordinate in one move. Tune it to obtain a practical acceptance rate.

property maximum_displacement

Maximum absolute displacement along each coordinate.

Type:

pint.Quantity

class MetropolisSimulator

Bases: pybind11_object

Equilibrate a non-overlapping hard-sphere configuration using Metropolis moves.

Parameters:

• domain (PackingDomain) – Fixed simulation domain and boundary conditions.

• initial_configuration (PackingConfiguration) – Valid, non-overlapping configuration, for example from RSASimulator.

• options (MetropolisOptions) – Sweep count, proposal displacement, and seed.

Notes

This sampler holds particle count, radii, and class labels fixed. It is an equilibrium workflow and is not a continuation of RSA deposition.

reset(self: PackLab.monte_carlo.metropolis.MetropolisSimulator) → None

Restore the supplied initial configuration.

run(self: PackLab.monte_carlo.metropolis.MetropolisSimulator) → object

Run all configured sweeps and return a PackingResult.

property sphere_configuration

Current configuration after accepted displacement moves.

property statistics

Metropolis displacement-move statistics.

class MetropolisStatistics

Bases: pybind11_object

Move statistics from an equilibrium hard-sphere Metropolis simulation.

Unlike PackingStatistics, these counters describe displacement moves, not RSA insertion attempts.

Results and visualisation

post_mpl_plot(plot_function)

Return a plotting wrapper with the former show convenience option.

class PackingResult(binding)

Bases: object

Output container for an RSA simulation run.

Holds arrays plus domain metadata, computed statistics, and plotting helpers.

property positions: ndarray

Get the sphere center positions as a NumPy array of shape (N, 3).

Returns:

Array of sphere center positions.

Return type:

np.ndarray

property radii: ndarray

Get the sphere radii as a NumPy array of shape (N,).

Returns:

Array of sphere radii.

Return type:

np.ndarray

property statistics

Get the simulation statistics.

Returns:

The simulation statistics.

Return type:

Statistics

property sphere_configuration

Get the sphere configuration object.

Returns:

The configuration of spheres in the simulation.

Return type:

SphereConfiguration

property domain

Get the simulation domain.

Returns:

The domain (box) dimensions and boundary conditions.

Return type:

Domain

compute_partial_pair_correlation_function(**kwargs) → tuple

Compute the partial pair correlation function g_ij(r).

Parameters:

**kwargs (dict) – Keyword arguments passed to the binding method.

Returns:

A tuple containing: - centers (Quantity): The radial distances with units. - g_ij (np.ndarray): The partial pair correlation values.

Return type:

tuple

property partial_volume_fractions: ndarray

Get the partial volume fractions of each component.

Returns:

Array of partial volume fractions.

Return type:

np.ndarray

property partial_volumes: ndarray

Get the partial volumes of each component.

Returns:

Array of partial volumes.

Return type:

np.ndarray

compute_pair_correlation_function(**kwargs) → None

Compute the total pair correlation function g(r).

Parameters:

**kwargs (dict) – Keyword arguments passed to the binding method.

property pair_correlation_centers: ndarray

Get the radial centers for the pair correlation function.

Returns:

Array of radial centers.

Return type:

np.ndarray

property pair_correlation_values: ndarray

Get the values of the pair correlation function g(r).

Returns:

Array of pair correlation values.

Return type:

np.ndarray

plot_centers_3d(maximum_points_3d: int = 10000) → Figure

Plot the sphere centers in a 3D scatter plot.

Parameters:

maximum_points_3d (int) – Maximum number of points to plot (subsampling if necessary).

Returns:

The matplotlib Figure object containing the 3D scatter plot.

Return type:

plt.Figure

plot_radius_distribution(bins: int = 40, density: bool = True, alpha: float = 0.85) → Figure

Plot the distribution of sphere radii.

Parameters:

• bins (int) – Number of histogram bins.

• density (bool) – Whether to normalize the histogram to form a probability density.

• alpha (float) – Transparency level for the histogram bars.

plot_slice_2d(slice_axis: Literal['x', 'y', 'z'] = 'z', slice_center_fraction: float = 0.5, slice_thickness_fraction: float = 0.08, maximum_circles_in_slice: int = 2500) → Figure

Plot a 2D slice of the sphere configuration.

Parameters:

• slice_axis (Literal["x", "y", "z"]) – Axis along which to take the slice.

• slice_center_fraction (float) – Fractional position along the slice axis where the slice is centered (0.0 to 1.0).

• slice_thickness_fraction (float) – Fractional thickness of the slice relative to the box length along the slice axis (0.0 to 1.0).

• maximum_circles_in_slice (int) – Maximum number of circles to plot in the slice (subsampling if necessary).

plot_pair_correlation(n_bins: int = 80, maximum_pairs: int = 300000) → Figure

Plot the partial pair correlation functions g_ij(r) obtained from the RSA configuration, overlaying all (i, j) curves on a single axis.

Parameters:

• n_bins (int) – Number of radial distance bins.

• maximum_pairs (int) – Number of Monte Carlo sampled pairs used for estimation.

Returns:

Figure containing the overlaid plot of all g_ij(r) curves.

Return type:

matplotlib.figure.Figure

Statistics and ensemble estimates

class PackingStatistics

Bases: pybind11_object

Summary statistics from an RSA simulation.

sphere_count

Number of accepted spheres.

Type:

int

packing_fraction_geometry

Packing fraction calculated from accepted sphere volumes.

Type:

float

total_runtime_seconds

Wall-clock runtime of the simulation.

Type:

float

print(self: PackLab.monte_carlo.statistics.PackingStatistics) → None

Print a human-readable statistics summary.

class EstimatorStatistics

Bases: pybind11_object

Aggregate diagnostics from the most recent PackingEstimator.estimate call.

requested_samples, completed_samples

Requested and successfully completed RSA realisations.

Type:

int

attempted_insertions, accepted_insertions, rejected_insertions

Totals over all completed realisations.

Type:

int

acceptance_rate

Accepted insertions divided by attempted insertions.

Type:

float

mean_sphere_count, mean_packing_fraction

Per-realisation means.

Type:

float

total_runtime_seconds, mean_runtime_seconds

Wall-clock simulation times.

Type:

float

print(self: PackLab.monte_carlo.estimator.EstimatorStatistics) → None

Print a tabular summary of the aggregate diagnostics.

class PackingEstimate

Bases: pybind11_object

Ensemble estimate of partial pair-correlation functions.

centers

Radial bin centers in meters.

Type:

pint.Quantity, shape (n_bins,)

mean_g, std_g

Mean and sample standard deviation of partial pair correlations.

Type:

numpy.ndarray, shape (n_species, n_species, n_bins)

class PackingEstimator

Bases: pybind11_object

Estimate pair-correlation statistics over repeated RSA simulations.

Parameters:

• domain (PackingDomain) – Simulation domain.

• radius_sampler (RadiusSampler) – Particle-radius distribution.

• options (RSAOptions) – RSA simulation configuration.

• number_of_bins (int) – Number of radial bins used for each estimate.

estimate(self: PackLab.monte_carlo.estimator.PackingEstimator, number_of_samples: SupportsInt | SupportsIndex, maximum_pairs: SupportsInt | SupportsIndex = 0, progress: bool = False, progress_interval: SupportsInt | SupportsIndex = 1) → PackLab.monte_carlo.estimator.PackingEstimate

Estimate partial pair correlations from independent RSA realizations.

Parameters:

• number_of_samples (int) – Number of independent packing realizations.

• maximum_pairs (int, default=0) – Pair-sampling limit per realization; zero selects the native default.

• progress (bool, default=False) – Print completed-sample progress and per-sample insertion diagnostics.

• progress_interval (int, default=1) – Print every this many completed samples when progress=True.

Returns:

Mean and standard deviation of the partial pair-correlation functions.

Return type:

PackingEstimate

print_statistics(self: PackLab.monte_carlo.estimator.PackingEstimator) → None

Print aggregate diagnostics from the most recent estimate.

property statistics

Aggregate diagnostics from the most recent estimate.

Analytical model

Mixture domain

Polydisperse cubic domain (analytical) with unit handling performed in the pybind11 wrapper.

This module exposes a C++ implementation of a cubic simulation domain containing a polydisperse population of spherical particles. The numerical core operates in base SI units (meters), while the Python interface accepts and returns Pint quantities.

Key conventions

• All lengths are stored internally in meters.

• Radii arrays are returned as Pint quantities in meters.

• Volumes are returned as Pint quantities in meter**3.

• Number densities are returned as Pint quantities in 1/meter**3.

• Fractions (number fractions and volume fractions) are dimensionless NumPy arrays.

The domain is specified by: - size: cubic side length - radii: per bin particle radii - volume_fraction: total occupied volume fraction - number_fractions: per bin number fractions (normalized to sum to 1)

The class provides deterministic conversion from fractional counts to integer particle counts using the configured rounding mode.

class PercusYevickDomain

Bases: pybind11_object

Cubic simulation domain containing a polydisperse population of spherical particles.

The numerical core stores all lengths in meters and all fractions as dimensionless floating point values. The wrapper converts Pint quantities to meters on input and reconstructs Pint quantities on output.

Parameters:

• size (pint.Quantity) – Side length of the cubic domain. Converted to meters internally.

• radii (pint.Quantity[ndarray]) – One dimensional array of particle radii defining the bins. Converted to meters internally.

• volume_fraction (float) – Total occupied volume fraction (dimensionless). Typical range is (0, 1].

• number_fractions (array_like of float) – Number fractions per bin. Must be non negative and sum to a positive value. The constructor normalizes them to sum to 1.

• rounding_mode (RoundingMode, default=floor) – Rounding mode used when converting expected counts to integers.

Notes

The model uses number fractions and a total volume fraction as its sole specification mechanism. Alternative specifications should use separate constructors. Returned arrays are NumPy arrays; quantities use SI units (meter, meter**3, 1/meter**3).

bins_table(self: PackLab.analytical.domain.PercusYevickDomain, precision: SupportsInt | SupportsIndex = 6) → str

Return the per bin table as a string.

Parameters:

precision (int, default=6) – Significant digits used for numeric formatting.

Returns:

Formatted table in base units.

Return type:

str

property mean_particle_volume_number_weighted

Number weighted mean particle volume.

Returns:

Mean per particle volume in meter**3 computed using number_fractions as weights.

Return type:

pint.Quantity

property number_fractions

Per bin number fractions.

Returns:

One dimensional array of number fractions (dimensionless). Always normalized to sum to 1.

Return type:

numpy.ndarray

property number_of_particles_per_radius

Per bin integer particle counts.

Returns:

One dimensional integer array of particle counts per radius bin.

Return type:

numpy.ndarray

property number_of_particles_total

Total inferred particle count.

This is computed from: - total_particle_volume = volume_fraction * volume - mean_particle_volume_number_weighted and then rounded according to rounding_mode.

Returns:

Total integer particle count.

Return type:

int

property particle_densities_per_radius

Per bin number densities.

Returns:

One dimensional array of number densities in 1/meter**3.

Return type:

pint.Quantity[ndarray]

property particle_density_total

Total number density over all bins.

Returns:

Total number density in 1/meter**3.

Return type:

pint.Quantity

property particle_volumes

Per bin particle volumes.

Returns:

One dimensional array of per particle volumes in meter**3 computed as (4/3)π r^3.

Return type:

pint.Quantity[ndarray]

print_bins(self: PackLab.analytical.domain.PercusYevickDomain, precision: SupportsInt | SupportsIndex = 6) → None

Print a per bin table describing the polydisperse population.

The table is printed in base units and includes: - radius (m) - per particle volume (m^3) - number fraction (dimensionless) - volume fraction (dimensionless) - particle count (integer) - number density (1/m^3)

Parameters:

precision (int, default=6) – Significant digits used for numeric formatting.

Return type:

None

property radii

Particle radii bins.

Returns:

One dimensional Pint quantity array of radii in meters.

Return type:

pint.Quantity[ndarray]

sample_radii(self: PackLab.analytical.domain.PercusYevickDomain, number_of_samples: SupportsInt | SupportsIndex, seed: object = None) → object

Sample particle radii according to number_fractions.

Parameters:

• number_of_samples (int) – Number of radii draws.

• seed (Optional[int]) – RNG seed. If None, a non deterministic seed is used.

Returns:

One dimensional Pint quantity array of sampled radii in meters.

Return type:

pint.Quantity[ndarray]

property size

PYDomain side length.

Returns:

The cubic side length as a Pint quantity in meters.

Return type:

pint.Quantity

property total_particle_volume

Total occupied particle volume.

Returns:

Total occupied volume in meter**3 computed as volume_fraction * volume.

Return type:

pint.Quantity

property volume

PYDomain volume.

Returns:

Total cubic volume in meter**3.

Return type:

pint.Quantity

property volume_fraction

Total occupied volume fraction.

Returns:

Dimensionless occupied volume fraction used to infer total particle volume.

Return type:

float

property volume_fraction_per_radius

Per bin occupied volume fractions.

Returns:

One dimensional array of per bin volume fractions (dimensionless).

Return type:

numpy.ndarray

class RoundingMode

Bases: pybind11_object

Rounding mode used when converting expected (non integer) particle counts to integers.

floor always rounds down. round rounds to the nearest integer. The selected mode affects particle counts and densities derived from the supplied number fractions.

Members:

floor

round

RoundingMode.name -> str

Wavenumber grids

Utilities for constructing reciprocal-space grids.

make_wavenumber_grid(radial_resolution: Any, maximum_distance: Any, samples_per_oscillation: int = 12) → Any

Create a reciprocal-space grid for the radial inverse transform.

Parameters:

• radial_resolution (pint.Quantity) – Smallest real-space feature to resolve. It determines the maximum wavenumber through pi / radial_resolution.

• maximum_distance (pint.Quantity) – Largest real-space distance at which the correlation function will be evaluated. It determines the wavenumber spacing.

• samples_per_oscillation (int, default=12) – Samples used for one period of the sinc kernel at maximum_distance. Values below 8 are not recommended.

Returns:

A uniformly spaced, zero-inclusive wavenumber grid with units of inverse length.

Return type:

pint.Quantity

Percus–Yevick solver

Percus Yevick mixture solver (C++ core, Pint handled in wrapper).

class PercusYevickResult

Bases: pybind11_object

Result of a multicomponent Percus-Yevick calculation.

radii

Particle radii in meters.

Type:

pint.Quantity, shape (n_species,)

densities

Number densities in inverse cubic meters.

Type:

pint.Quantity, shape (n_species,)

wavenumber

Wavenumber grid in inverse meters used for the inverse transform.

Type:

pint.Quantity, shape (n_wavenumber,)

g, h, H

Real-space and reciprocal-space correlation functions.

Type:

numpy.ndarray

class PercusYevickSolver

Bases: pybind11_object

Solve the multicomponent Percus-Yevick hard-sphere model.

Parameters:

• densities (pint.Quantity, shape (n_species,)) – Species number densities, convertible to inverse cubic meters.

• radii (pint.Quantity, shape (n_species,)) – Species radii, convertible to meters.

• wavenumber (pint.Quantity or "auto", default="auto") – Strictly increasing wavenumber grid, convertible to inverse meters. "auto" derives a grid when compute() receives the requested distances.

• radial_resolution (pint.Quantity or None, optional) – Smallest real-space feature to resolve when wavenumber="auto". By default this is one twentieth of the smallest particle radius.

• samples_per_oscillation (int, default=12) – Reciprocal-space samples per sinc-kernel oscillation at the largest requested distance when wavenumber="auto". Must be at least 8.

Notes

The inverse radial transform must resolve the sinc kernel at the largest requested distance. compute() emits a RuntimeWarning when this grid has fewer than eight samples per kernel oscillation. Use PackLab.analytical.make_wavenumber_grid to construct a balanced grid.

compute(self: PackLab.analytical.solver.PercusYevickSolver, distances: object) → PackLab.analytical.solver.PercusYevickResult

Compute the correlation functions on a real-space distance grid.

Parameters:

distances (pint.Quantity, shape (n_r,)) – Radial positions, convertible to meters.

Returns:

Computed structural and correlation functions.

Return type:

PercusYevickResult

Scattering tools

Data containers

class ScatteringData(S1: Any, S2: Any, k: Any, Csca: Any, phi: Any)

Bases: object

Far-field scattering data for a single particle diameter.

The quantity-bearing fields intentionally use Any because TypedUnit/Pint quantities are generic at runtime and may wrap scalar or array values.

plot(*, tight_layout: bool = True)

Plot the magnitudes of the two scattering amplitudes.

class ScatteringDataset(iterable=(), /)

Bases: list

Container for multi size scattering data and mixture level post processing.

This class stores ScatteringData instances, one for each diameter. It also provides mixture formulas for a number density distribution over diameters and an inter particle correlation term H(wavenumber).

e_theta, e_phi

Unit basis vectors for the polarization basis used in get_F_matrix. This code treats the scattering amplitude as a 2 component vector in a (theta, phi) polarization basis.

Type:

NDArray, shape (2,)

k

Optical wavenumber in the surrounding medium (1/length).

phi

Polar scattering angle grid in radians. Expected range is [0, pi].

S1, S2

Arrays assembled by process() from each per diameter element.

Csca

Per diameter scattering cross section assembled by process().

Notes

Derived arrays are refreshed automatically by the mixture methods. Calling process() explicitly remains useful when inspecting those arrays.

process()

Stack per diameter fields into dense arrays.

After calling this method, the object provides:

• self.S1: ndarray, shape (N, Pphi)

• self.S2: ndarray, shape (N, Pphi)

• self.Csca: Pint Quantity, shape (N,), units meter**2

where N is the number of diameters (len(self)) and Pphi is the number of angular samples in each per diameter result.

Return type:

None

get_alpha_beta_factor(densities: ndarray[tuple[int, ...], dtype[_ScalarType_co]]) → tuple[ndarray[tuple[int, ...], dtype[_ScalarType_co]], ndarray[tuple[int, ...], dtype[_ScalarType_co]]]

Build number density factors for mixture sums.

Parameters:

densities (NDArray, shape (N,)) – Species number densities n_alpha. Units should be 1/length**3. N must match the number of size classes stored in this object.

Returns:

• n_alpha (NDArray, shape (N,)) – Same array as input (returned for convenience).

• sqrt_alpha_beta (NDArray, shape (N, N)) – Matrix with entries sqrt(n_alpha * n_beta). This is commonly used in symmetric mixture formulations to weight cross correlations.

Notes

This uses the outer product and then takes a square root:

sqrt_alpha_beta[a, b] = sqrt(n_alpha[a] * n_alpha[b])

get_F_matrix(theta_points: int = 100) → tuple

Construct the vector scattering amplitude tensor F for each size, angle, and azimuth.

This builds a 2 component vector amplitude in the (e_theta, e_phi) basis using the convention in your code:

F(θ_azimuth, φ_polar) = (i / k) * [ e_theta * S2(φ) * cos(θ) - e_phi * S1(φ) * sin(θ) ]

Parameters:

theta_points – Number of azimuthal sampling points in [0, 2*pi]. This is a rotation around the incident axis.

Returns:

• F (NDArray, complex, shape (2, N, Pphi, Ptheta)) – Vector amplitude tensor.

  Index meaning: * first axis (size 2): polarization component (theta, phi) * second axis: size class index (N) * third axis: polar scattering angle samples (Pphi), from self.phi * fourth axis: azimuthal angle samples (Ptheta), generated internally

• theta (NDArray, shape (Ptheta,)) – Azimuthal angle grid in radians spanning [0, 2*pi].

Notes

The name theta here refers to an azimuthal angle. The polar scattering angle is stored in self.phi.

get_mu_independant(densities: ndarray[tuple[int, ...], dtype[_ScalarType_co]])

Compute the independent scattering contribution to the attenuation coefficient.

Parameters:

densities (NDArray, shape (N,)) – Species number densities n_alpha. Units should be 1/length**3.

Returns:

Scalar attenuation coefficient contribution with units 1/length.

Return type:

mu_independant_scattering

Notes

This implements a standard independent scattering approximation term:

mu_s,ind = sum_alpha n_alpha * Csca_alpha

The supplied density vector is validated against the number of species.

get_mu_dependant(densities: ndarray[tuple[int, ...], dtype[_ScalarType_co]], H: ndarray[tuple[int, ...], dtype[_ScalarType_co]], wavenumber: ndarray[tuple[int, ...], dtype[_ScalarType_co]], theta_points: int = 150)

Compute the dependent scattering contribution using inter particle correlations H(wavenumber).

Parameters:

• densities (NDArray, shape (N,)) – Species number densities n_alpha. Units should be 1/length**3.

• H (NDArray, shape (N, N, Pp)) – Mixture total correlation function in reciprocal space, evaluated on the wavenumber grid. This should correspond to the same ordering as size classes.

• wavenumber (NDArray, shape (Pp,)) – Reciprocal space radial grid (1/length). Must cover the range needed to evaluate H at wavenumber = 2k sin(φ/2) for φ in [0, pi].

• theta_points – Number of azimuthal samples for the integral over the rotation angle.

Returns:

Scalar attenuation coefficient contribution. Units depend on the exact normalization conventions of F and H used in the integrand, but it is intended to be in 1/length.

Return type:

mu_dependant_scattering

Notes

The core integrand is:

term[p_index, theta_index] = sum_{a,b} sqrt(n_a n_b) F_a(wavenumber,theta) F*_b(wavenumber,theta) H_ab(wavenumber)

then integrated over polar angle (self.phi) and azimuth (theta).

get_mu(densities: ndarray[tuple[int, ...], dtype[_ScalarType_co]], H: ndarray[tuple[int, ...], dtype[_ScalarType_co]], wavenumber: ndarray[tuple[int, ...], dtype[_ScalarType_co]], theta_points: int = 150) → ndarray[tuple[int, ...], dtype[_ScalarType_co]]

Compute total scattering attenuation coefficient mu_s.

Parameters:

• densities (NDArray, shape (N,)) – Species number densities.

• H (NDArray, shape (N, N, Pp)) – Correlation function in wavenumber space.

• wavenumber (NDArray, shape (Pp,)) – wavenumber space grid (1/length).

• theta_points – Azimuthal sampling count.

Returns:

Total scattering attenuation coefficient (intended units 1/length), computed as:

mu = mu_independant + mu_dependant

Return type:

mu

interpolate_last_axis_linear(H: ndarray, wavenumber: ndarray, evaluation_wavenumber: ndarray) → ndarray

Linear interpolation of H(…, wavenumber) onto evaluation_wavenumber, along the last axis.

Parameters:

• H (ndarray, shape (..., Pp)) – Values sampled on grid wavenumber along the last axis.

• wavenumber (ndarray, shape (Pp,)) – Strictly increasing sample points.

• evaluation_wavenumber (ndarray, shape (Pq,)) – Query points.

Returns:

Interpolated values.

Return type:

ndarray, shape (…, Pq)

get_interpolated_H(H, wavenumber)

Interpolate H(wavenumber) onto the scattering wavevector magnitude wavenumber = 2k sin(φ/2).

Parameters:

• H (array like, shape (N, N, Pp)) – Correlation tensor sampled on the input wavenumber grid.

• wavenumber (array like, shape (Pp,)) – Radial reciprocal space grid (1/length).

Returns:

interpolated_H – H evaluated at the scattering wavenumber for φ in [0, pi], where:

wavenumber(φ) = 2 k sin(φ/2)

Return type:

NDArray, shape (N, N, Pphi)

Notes

This method assumes: * self.phi spans [0, pi] and is the polar scattering angle * wavenumber is provided in increasing order

get_phase_function(densities: ndarray[tuple[int, ...], dtype[_ScalarType_co]], H: ndarray[tuple[int, ...], dtype[_ScalarType_co]], wavenumber: ndarray[tuple[int, ...], dtype[_ScalarType_co]], theta_points: int = 150) → tuple[ndarray[tuple[int, ...], dtype[_ScalarType_co]], ndarray[tuple[int, ...], dtype[_ScalarType_co]], ndarray[tuple[int, ...], dtype[_ScalarType_co]]]

Compute the mixture phase function including dependent scattering corrections.

Parameters:

• densities (NDArray, shape (N,)) – Species number densities n_alpha.

• H (NDArray, shape (N, N, Pp)) – Correlation tensor in reciprocal space on the wavenumber grid.

• wavenumber (NDArray, shape (Pp,)) – wavenumber space grid (1/length).

• theta_points – Number of azimuthal samples.

Returns:

• phi (ndarray, shape (Pphi,)) – Polar angle grid in radians.

• theta (ndarray, shape (Ptheta,)) – Azimuthal angle grid in radians.

• phase_function (NDArray, shape (Pphi, Ptheta)) – Phase function indexed first by polar angle, then by azimuth.

Notes

The phase function is assembled as:

• Independent term: sum of n |F|^2 over species.

• Dependent term: weighted cross terms involving F, its conjugate, and the correlation tensor H.

Amplitude generation

compute_scattering_amplitudes(wavelength: Length, diameters: Length, material: RefractiveIndex, medium: RefractiveIndex, phi: Angle, plot: bool = False, polarization: float = <Quantity(0, 'degree')>, debug_mode: bool = False) → ScatteringDataset

Compute far field amplitude scattering functions S1 and S2 for a set of sphere diameters.

This helper constructs a Gaussian source and a Sphere scatterer for each diameter, calls Setup.get_s1s2(…), and stores the resulting objects in a ScatteringDataset container.

Parameters:

• phi – Angle grid relative to the transverse plane. The amplitudes are evaluated at phi + pi / 2 and the returned dataset stores this physical polar-angle grid in the interval [0, pi].

• diameters – Iterable of particle diameters. Each element is expected to be a Pint Quantity with length units (for example 6 * ureg.nanometer).

• plot – If True, call data.plot(…) for each diameter.

Returns:

A list like container where each element is the return value of Setup.get_s1s2(angles=…). Additional attributes are attached:

• datas.k: optical wavenumber of the source (1/length)

• datas.phi: angular sampling grid taken from the first result (radians)

• after calling datas.process(), arrays S1, S2, Csca are available

Return type:

ScatteringDataset

Notes

This function sets two extra attributes on each data element:

• data.k is set to source.wavenumber

• data.Csca is set to scatterer.Csca

The ScatteringDataset instance also receives k and phi from the last and first element respectively. If you want stricter correctness, you should assert that all returned data.phi grids are identical.

Phase-function plots

plot_phase_function_3d(phi: ndarray, theta: ndarray, phase_function: ndarray, *, mode: Literal['spherical', 'surface'] = 'spherical', normalize: bool = True, use_magnitude: bool = True) → Figure

Plot a 3D representation of the phase function P(phi, theta).

Parameters:

• phi – Polar scattering angle grid in radians, expected in [0, pi].

• theta – Azimuthal angle grid in radians, expected in [0, 2*pi].

• phase_function – Phase function sampled on (phi, theta) or (theta, phi).

• mode – “spherical”: map intensity to radius on a sphere and render a 3D surface. “surface”: plot a 3D surface with axes (theta, phi, P).

• normalize – If True, normalize the plotted values by their maximum for visual clarity.

• use_magnitude – If True, plot abs(P). If False, plot real(P).

Returns:

Matplotlib figure containing the 3D plot.

Return type:

figure

Notes

In “spherical” mode, the surface is parameterized as:

x = r * sin(phi) * cos(theta)
y = r * sin(phi) * sin(theta)
z = r * cos(phi)

where r is the (optionally normalized) phase function.

plot_phase_function_2d_projection(phi: ndarray, theta: ndarray, phase_function: ndarray, *, projection: Literal['azimuth_average', 'heatmap'] = 'azimuth_average', use_magnitude: bool = True, normalize: bool = False) → Figure

Plot a 2D representation of the phase function.

Parameters:

• phi – Polar scattering angle grid in radians, expected in [0, pi].

• theta – Azimuthal angle grid in radians, expected in [0, 2*pi].

• phase_function – Phase function sampled on (phi, theta) or (theta, phi).

• projection – “azimuth_average”: plot P_avg(phi) = (1/2pi) * integral P(phi,theta) dtheta. “heatmap”: plot a 2D image of P(phi,theta) with axes theta and phi.

• use_magnitude – If True, plot abs(P). If False, plot real(P).

• normalize – If True, normalize plotted values by their maximum.

Returns:

Matplotlib figure containing the 2D plot.

Return type:

figure

Notes

The azimuth averaged projection is the closest analogue to typical S1 and S2 plots because it reduces the 2D angular dependence into a 1D curve versus phi.