aci-rssa 0.0.9

Fast, streaming reader for MCNP RSSA (surface source) files, with a Rust core and a Polars-based Python API.

RSSA files record every particle that crosses a given surface during an MCNP (or D1S-UNED) simulation, and can easily reach tens or hundreds of gigabytes. aci_rssa memory-maps the file and decodes track records in parallel (via rayon), handing the result to Python as a Polars DataFrame — either all at once, or as a lazy, batch-at-a-time stream for files too large to fit in memory.

Installation

pip install aci_rssa

Prebuilt wheels are published for Linux, macOS, and Windows (x86_64/arm64) on CPython ≥3.10 — no Rust toolchain needed to install.

Quickstart

from aci_rssa import RSSA

rssa = RSSA.read_from_file("surface_source.w")
print(rssa)  # summary: surfaces, track/history counts, source code

rssa.tracks          # polars.DataFrame: history, particle_type, weight, energy,
                      # time, x, y, z, u, v, w, surface_id
rssa.neutron_tracks   # tracks filtered to neutrons
rssa.photon_tracks    # tracks filtered to everything else (typically photons)

Reading only the header (fast, never touches the potentially huge track data):

from aci_rssa import FileParameters

parameters = FileParameters.read_from_file("surface_source.w")
parameters.nrss           # number of tracks recorded
parameters.surfaces       # surfaces tracks were recorded on
parameters.surface_ids    # just their ids

surface = parameters.get_surface(63)
surface.mnemonic          # "cz"
surface.axis, surface.center, surface.radius   # "z", (0.0, 0.0), 164.3
print(surface.describe()) # surface 63: cz (cylinder about the z axis, radius 164.30 cm)

Surface coefficients are exposed as MCNP stores them (surface.parameters, e.g. a cylinder's squared radius); the properties above do the conversion. Planes give normal and offset instead.

Very large files

For files too large to load eagerly, scan_tracks returns a Polars LazyFrame backed by a Rust reader that decodes the file batch by batch, driven by Polars' query engine — filters and column selection are pushed down so you only decode what you need:

from aci_rssa import scan_tracks
import polars as pl

(
    scan_tracks("surface_source.w")
    .filter(pl.col("energy") > 1.0)
    .select("x", "y", "z", "energy")
    .collect(engine="streaming")
)

RSSAPlot and RSSASpectraPlot accept a LazyFrame directly, so a scan_tracks() result can be plotted without ever loading the file eagerly — build them with RSSAPlot.for_surface() instead of going through RSSA.read_from_file(), which requires an in-memory RSSA:

from aci_rssa import FileParameters, RSSAPlot, neutrons, scan_tracks

parameters = FileParameters.read_from_file("surface_source.w")  # header only
plot = RSSAPlot.for_surface(scan_tracks("surface_source.w"), parameters, 63)

(
    plot.filter(neutrons())
    .with_bin_width(10)
    .particle_flux(1e17)
    .show()
)

Every aggregation (with_bin_width, particle_flux, …) streams through the file batch by batch rather than materializing it.

Plotting

rssa.plot_surface(63).with_bin_width(10).particle_flux(1e17).show()
rssa.plot_spectra().spectrum().show()

There are three steps, and they are separate on purpose: a projection decides how the surface is unrolled into 2D, a plot bins the tracks on it, and a result (CurrentMap, VectorCurrent, Spectrum) is a computed value that knows how to render itself.

Projections

plot_surface() takes the geometry from the file header: it keeps only the tracks that crossed that surface, and picks the projection from the surface's own type — a cylinder is unrolled around its own axis and centre over its full circumference, a px/py/pz plane is plotted on the two axes spanning it, and a tilted p plane gets an orthonormal basis of the plane itself. The surface id may be omitted only when the file recorded a single surface.

This is the only entry point to a current map, because the projection is what makes the map's numbers right: every quantity divides each bin's summed weight by dx * dy, so a projection that doesn't preserve distances on the surface — dropping a tilted plane onto two global axes shrinks one axis by the cosine of the tilt; measuring a cylinder's angle about the coordinate origin instead of its own axis smears it — gives a plausible plot that is wrong by that factor.

The escape hatch is therefore to override the geometry, not to skip it:

from aci_rssa import Cylinder

rssa.plot_surface(63, projection=Cylinder(axis="z", radius=164.3, center=(0.0, 0.0)))

Cylinder, AxisAlignedPlane and TiltedPlane are the built-in projections; a surface type none of them fits (a sphere, a torus) raises rather than being rastered onto two global axes. Unless verify=False, a sample of the tracks is checked against the surface's equation, so a header that means something other than what is assumed here fails loudly instead of producing a plausible, wrong plot.

Selecting tracks

The tracks are a Polars LazyFrame, so selection is done with ordinary Polars expressions via filter(). aci_rssa only supplies the two predicates that need knowledge of the file rather than of Polars — the packed neutron code, and surface ids (matched irrespective of sign, and refused on files that don't record them per track):

import polars as pl
from aci_rssa import neutrons, on_surface

plot.filter(neutrons(), pl.col("energy") > 1.0, pl.col("z").is_between(-600, 800))

filter() returns a new plot, so one object can feed several derived ones.

Results

plot = rssa.plot_surface(63).filter(neutrons()).with_bin_width(10)

flux = plot.particle_flux(1e17)        # CurrentMap: scalar flux, #/cm2/s per bin
current = plot.particle_current(1e17)  # CurrentMap: crossing rate per bin
net = plot.net_current(1e17)           # CurrentMap: net leakage per bin, signed
arrows = plot.vector_current(1e17)     # VectorCurrent: J and anisotropy per bin
angles = plot.angular_distribution(1e17)  # AngularDistribution: ψ(Ω) itself

flux.with_parameters(vmin=1e6, vmax=1e12).show()
flux.ratio_to(other_flux).save("difference.png")

Each result carries its own values, bins and rendering options; with_parameters() takes keyword overrides (keeping the labels and title derived from the surface) or a whole PlotParameters to replace them. figure() hands back the Matplotlib figure and axes for anything the parameters don't cover.

Pass errors=True to any of the maps to estimate its statistical error in the same pass, then plot it with error_map():

flux = plot.particle_flux(1e17, errors=True)
flux.relative_errors      # array, or None when it wasn't asked for
flux.error_map().show()   # the errors as a map of their own

The error is carried on the map rather than computed by a separate call, so it can only ever belong to the quantity it was estimated for — the error of a flux is not the error of a current.

Which quantity do you want?

All four come from one identity: a set of surface crossings samples the angular flux weighted by |μ| — head-on crossings are over-represented, because a patch of surface presents its full area to them. Summing w·f(Ω) over tracks then estimates ∫f(Ω)|μ|ψ(Ω)dΩ, and the choice of f is the choice of quantity:

If you are about to multiply a map by a response function, you want particle_flux(), not particle_current(). They are not interchangeable: for an isotropic field the mean of 1/|μ| over the crossing distribution is exactly 2, so a current map understates the flux by a factor of two where the field is diffuse while agreeing with it where the field is collimated — it distorts the shape of the map, not just its scale. This is MCNP's F1 vs F2 distinction, and aci_rssa computes both from the surface normal the projection already carries.

VectorCurrent.table's anisotropy column is |J|/Φ, the flux-to-current ratio: 1 for a beam, 0.5 for a half-isotropic distribution, ~0 where directions cancel. It also bounds the gap above — since 1/|μ| is convex, Jensen's inequality gives Φ / crossing rate ≥ 1 / anisotropy.

Tracks approaching tangency (|μ| → 0) would send 1/|μ| to infinity. Following MCNP, particle_flux() scores a track below |μ| = 1e-3 as if it were 5e-4, and grazing_fraction() maps how much of each bin's flux came from tracks that were capped. That blow-up is real physics — a shallow track really does deposit more path length in a thin volume at the surface — so it is capped, not clipped away.

Seeing the directions themselves

The four quantities above are all moments of one thing, the angular flux ψ(Ω). angular_distribution() estimates ψ(Ω) directly, which is worth doing wherever a moment is not enough — and one place it never is, is a mean direction: vector_current() cannot tell an isotropic field from two opposing beams, since both average to |J| = 0. Where anisotropy is low, this is what says which of the two you have.

angles = plot.angular_distribution(1e17)
angles.show()        # two equal-area polar disks, outgoing and incoming
angles.scalar_flux   # ∫ψ dΩ, the scalar flux averaged over the patch

Directions are binned on μ and on the azimuth φ in the surface's own tangent plane, aligned with the map's axes: φ = 0 points the way the plot's x axis grows, so a lobe leaning up-and-right in the disk means particles heading up-and-right on the current map — including on a cylinder, where the frame turns from bin to bin. Sectors are centred on the named directions rather than divided by them, so a beam aimed along an axis reads as one petal instead of two halves.

The radius is √(2(1−|μ|)) — head-on at the centre, grazing at the rim. That is the Lambert equal-area radius, chosen so each cell's area on the page is proportional to the solid angle it stands for: plotting the polar angle as the radius instead would inflate the near-normal directions into a wide bullseye, the same area-encoding bias that makes a bar-length rose diagram overstate its longest petal. The μ grid is uniform for the same reason, so every cell subtends the same solid angle and the colours are directly comparable — as are the two disks, which share one scale.

Each track enters its cell as w/|μ|, exactly as in particle_flux(). This is not cosmetic: crossings sample ψ weighted by |μ|, so summing raw weight per direction would draw a dent at grazing angles that is an artefact of how a surface samples a field rather than a feature of the field itself — and grazing angles are precisely where an angular plot gets read.

Unlike the maps, this is a pooled quantity: every track inside the bins contributes equally to one histogram, with no per-bin averaging. The bins therefore only choose the region — filter or narrow them to zoom in — and the area ψ is divided by comes from the extent the tracks span, not from the bin grid, so bin_width does not change the result. (It would otherwise: with_bin_width() pads its last bin past the data, so a 100 cm bin over tracks spanning 60 cm claims a strip two thirds empty and dilutes every value by 1.67. A per-bin map is unaffected — there the padding lands in one edge bin that really does cover surface no track reached.) Pass area= to state the patch explicitly when the region is a deliberate choice rather than wherever the tracks landed:

plot.angular_distribution(1e17, area=200 * 200)   # a defined window, empty parts and all

Migrating from 0.0.7

The default colormap changed from jet to plasma, and relative errors now render on a linear scale from 0 rather than a log one.

Two changes affect numbers, not just names:

• vector_current() now returns the real J. It used to sum w·Ω with no 1/|μ|, which is the first moment of the crossing distribution rather than of the flux — correct only at normal incidence, and short by the incidence cosine everywhere else. Its scalar_current column is now crossing_rate, joined by flux, net_current, and an anisotropy that is |J|/Φ rather than |J|/crossing rate.
• particle_current() is unchanged, but is now one of four quantities rather than the only one. If you were reading it as a flux, switch to particle_flux(); if you were reading it as net leakage, switch to net_current().

Development

This is a Cargo workspace (rssa-core, the pure-Rust parsing/decoding library, and rssa-python, its PyO3 bindings) plus a Python package in python/, built together with maturin.

python -m venv .venv && source .venv/bin/activate
pip install "maturin>=1.9.4,<2.0"
maturin develop --extras test,dev   # builds the extension module, installs test/dev deps

cargo test                          # Rust tests (rssa-core)
pytest                              # Python tests
ty check python/aci_rssa            # type checking

License

Licensed under the European Union Public Licence v1.2 (EUPL-1.2).