Simulation

Track etching simulator for CR-39 nuclear track detectors.

This module provides TrackSimulator, which computes the 2-D dose map, etch-rate map, and etchant arrival-time map for a single ion track in CR-39, then extracts observable track geometry (radius, depth, contour).

class tracketch.simulation.simulator.TrackSimulator(particle_name: str, start_energy_MeV_u: float, etch_model: EtchRateModel | None = None, RDD_name: str = 'Cucinotta', arrival_time_method_name: str = 'dijkstra_numba', dijkstra_connectivity: int = 8, rz_lims_dict: dict | None = None, r_min_um: float | None = None, r_max_um: float | None = None, z_max_um: float | None = None, n_points_r: int | None = None, n_points_z: int | None = None, material_name: str = 'CR39', stopping_power_source: str = 'SRIM', n_straggling_sigma: int = 1, n_uniform_multiplier: int = 3, n_jobs: int = 1, log_level: int | str | None = None, normalise_to_LET: bool = False, logscale_r: bool = True, theta_deg: float = 0.0)[source]

Bases: object

Simulate chemical etching of a single ion track in CR-39.

The simulation proceeds in three stages:

Dose map – radial dose distribution (RDD) integrated along the ion path, accounting for energy loss and optional range straggling.
Etch-rate map – the calibrated etch_model converts local dose to local etch velocity, with optional debris-damping correction.
Arrival-time map – shortest-time wavefront propagation (Dijkstra or Fast Marching) from the detector surface into the bulk.

From the arrival-time map, iso-time contours, track radii, and track depths are extracted for any desired etching duration.

Parameters:

particle_name (str) – Ion species identifier in the form "<A><symbol>", e.g. "12C", "1H", "56Fe", "238U". When stopping_power_source="SRIM" (default) the particle must be one of tracketch.SRIM_PARTICLES; switch to stopping_power_source="libamtrack" to use any nuclide recognised by libamtrack.
start_energy_MeV_u (float) – Kinetic energy at the detector surface in MeV per nucleon.
etch_model (tracketch.etching.etch_rate_model.EtchRateModel or None) – Calibrated dose -> etch-rate model. If None, a default model is loaded ("Doerschel_etching").
RDD_name (str) – Radial dose distribution model name (default "Cucinotta").
arrival_time_method_name (str) – Wavefront algorithm: "dijkstra" (default, auto-selects fastest available backend), "dijkstra_cpp", "dijkstra_numba", or "fmm".
dijkstra_connectivity (int) – Neighbour connectivity for Dijkstra: 8, 16, or 32.
r_min_um (float or None) – Minimum radial coordinate in um (default 1e-4).
r_max_um (float or None) – Maximum radial coordinate in um (default 20).
z_max_um (float or None) – Maximum depth coordinate in um (default 40).
n_points_r (int or None) – Number of radial grid points (default 400).
n_points_z (int or None) – Number of depth grid points (default 100).
rz_lims_dict (dict or None) – Legacy grid specification dict. Individual keyword arguments above take precedence over keys in this dict when both are supplied.
material_name (str) – Target material – "CR39" (default) or "water".
stopping_power_source (str) – Stopping-power database: "SRIM" (default, tabulated data for common ions) or "libamtrack" (supports any nuclide, uses the PSTAR/ASTAR parametrisation).
n_straggling_sigma (int) – Number of longitudinal-straggling sigma to include (0 = none).
n_uniform_multiplier (int) – Resolution multiplier when regridding to uniform spacing for FMM.
n_jobs (int) – Number of parallel workers for the RDD/dose-map computation. -1 uses all available CPU cores (via joblib). Parallel execution benefits large grids (n_points_r x n_points_z > ~50 000) or stopping_power_source="libamtrack". For the default grid the inter-process overhead outweighs the gain; keep n_jobs=1 (sequential) unless you have increased the grid resolution. Default 1.
log_level (int or str or None) – Set the logging level for the tracketch logger, e.g. logging.DEBUG, logging.INFO, or the string "DEBUG". When None (default) the logger is left untouched.

dose_map

Local dose in Gy.

Type:: ndarray, shape (n_z, n_r)

etch_rate_map

Local etch velocity in um/hr (with debris damping applied, if enabled).

Type:: ndarray, shape (n_z, n_r)

etch_rate_map_nodebris

Local etch velocity in um/hr before any debris damping. Identical to etch_rate_map when debris damping is disabled.

Type:: ndarray, shape (n_z, n_r)

arrival_time_map

Shortest etchant arrival time in hours.

Type:: ndarray, shape (n_z, n_r)

CSDA_range_um

Continuous-slowing-down-approximation range in um.

Type:: float

energy_profile_MeV_u

Kinetic energy vs. depth (with straggling).

Type:: ndarray

LET_profile_keV_um

Linear energy transfer vs. depth (with straggling).

Type:: ndarray

etch_model

The active dose -> etch-rate model.

Type:: tracketch.etching.etch_rate_model.EtchRateModel

Notes

Track radius and depth

After construction, extract the observable track geometry for any etching duration with get_track_radius_um() and get_track_length_um():

sim = TrackSimulator(particle_name="12C", start_energy_MeV_u=270.0)

radius_um = sim.get_track_radius_um(etch_time_h=3.0)
depth_um  = sim.get_track_length_um(etch_time_h=3.0)  # pit depth below surface

Both return float("nan") when the track is not yet open (too short an etching time) or when the ion stops before reaching the detector surface.

Track contour

For the full 2-D pit shape use get_iso_time_contour(), which returns (r, z) arrays in um:

r, z = sim.get_iso_time_contour(etching_time_h=3.0)

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(r, z)
ax.plot(-r, z)   # mirror for full axisymmetric view
ax.invert_yaxis()
ax.set_aspect("equal")

Uncertainties

Uncertainties are propagated from the calibrated v(d) anchor-velocity uncertainties. These must first be estimated from a calibration dataset using calibration.lib_optimiser.estimate_parameter_uncertainties() and then stored on the etch model:

from calibration.lib_optimiser import estimate_parameter_uncertainties

result = estimate_parameter_uncertainties(
    etch_model=sim.etch_model,
    models_dict=my_calibration_data,  # same dict used during fitting
)
sim.etch_model.anchor_velocity_uncertainties_um_h = (
    result["anchor_velocity_uncertainties_um_h"]
)
sim.etch_model.V_max_uncertainty_um_h = result.get("V_max_uncertainty_um_h")

Once the uncertainties are attached, use the _with_uncertainty variants:

radius_um, radius_unc_um = sim.get_track_radius_um_with_uncertainty(
    etch_time_h=3.0,  # etching duration
    n_sigma=1.0,       # propagate +/-1 sigma of v(d) anchors
)

result = sim.get_track_length_um_with_uncertainty(
    etch_time_h=3.0,
    n_sigma=1.0,
)
if result is not None:
    depth_um, depth_unc_um = result

print(f"radius = {radius_um:.2f} +/- {radius_unc_um:.2f} um")
print(f"depth  = {depth_um:.2f} +/- {depth_unc_um:.2f} um")

Both methods return nan uncertainties when a bound contour cannot be extracted. get_track_radius_um_with_uncertainty() also floors the uncertainty at half the minimum radial grid step, since the contour extraction cannot resolve sub-cell differences.

For a full uncertainty band around the track contour see get_iso_time_contour_with_uncertainty().

update_etch_model_and_recalculate(new_etch_model: EtchRateModel | None = None) → None[source]

Replace the etch model and recompute etch-rate & arrival-time maps.

The dose map is not recomputed (it depends only on physics, not the etch model).

Parameters:: new_etch_model (tracketch.etching.etch_rate_model.EtchRateModel or None) – New model. If None, recomputes using the current model (useful after in-place parameter changes).

recalculate_from_current_etch_model() → None[source]

Recompute maps using the current etch model.

Alias for update_etch_model_and_recalculate(None).

get_iso_time_contour(etching_time_h: float) → tuple[ndarray, ndarray][source]

Extract the iso-time contour from the arrival-time map.

Parameters:: etching_time_h (float) – Etching duration in hours.
Returns:: r_coords, z_coords – Coordinates along the longest contour segment.
Return type:: tuple[ndarray, ndarray]

get_iso_time_contour_with_uncertainty(etching_time_h: float, n_sigma: float = 1.0, n_band_points: int = 200) → tuple[tuple[ndarray, ndarray], ndarray, ndarray, ndarray] | None[source]

Iso-time contour with uncertainty band from v(d) anchor uncertainties.

Computes arrival-time maps for v(d)+/-n_sigma*dV, extracts the three contours, and interpolates the upper/lower bounds onto a common z grid suitable for ax.fill_betweenx.

Parameters:

etching_time_h (float) – Etching duration in hours.
n_sigma (float) – Number of standard deviations to propagate (default 1).
n_band_points (int) – Number of points on the common z grid for the band (default 200).

Returns:

tuple of ((r_central, z_central), z_band, r_lo, r_hi) or None.
Returns *None (and logs a warning) when the etch model carries no*
anchor velocity uncertainties. Assign the output of
estimate_parameter_uncertainties to
etch_model.anchor_velocity_uncertainties_um_h first.

get_iso_time_contour_with_normal_uncertainty(etching_time_h: float, n_sigma: float = 1.0, n_band_points: int = 200, min_uncertainty_um: float = 0.15) → tuple[tuple[ndarray, ndarray], tuple[ndarray, ndarray], tuple[ndarray, ndarray]] | None[source]

Iso-time contour with uncertainty expressed along local contour normals.

Unlike get_iso_time_contour_with_uncertainty (which returns a radial r(z) band for fill_betweenx), this method computes the uncertainty envelope perpendicular to the contour itself. This keeps the band visible and geometrically meaningful even where the contour is horizontal.

Return type:: tuple of ((r_c, z_c), (r_lo, z_lo), (r_hi, z_hi)) or None.

get_track_radius_um(etch_time_h: float, threshold_percent: float = 5.0) → float[source]

Track opening radius from iso-time contour deviation.

Parameters:

etch_time_h (float) – Etching duration in hours.
threshold_percent (float) – Deviation threshold as percentage of bulk depth (default 5 %).

Returns:

Track radius in um, or nan if no contour is found.

Return type:

float

get_iso_time_contour_nodebris(etching_time_h: float) → tuple[ndarray, ndarray][source]

Iso-time contour computed from the etch-rate map without debris damping.

Useful for comparing predicted track shapes with and without the debris-transport correction. The arrival-time map is computed on demand from etch_rate_map_nodebris each time this method is called (not cached).

Parameters:: etching_time_h (float) – Etching duration in hours.
Returns:: r_coords, z_coords – Coordinates along the iso-time contour.
Return type:: tuple[ndarray, ndarray]

get_track_radius_um_nodebris(etch_time_h: float, threshold_percent: float = 5.0) → float[source]

Track opening radius computed without debris damping.

Convenience wrapper around get_iso_time_contour_nodebris().

Parameters:

etch_time_h (float) – Etching duration in hours.
threshold_percent (float) – Deviation threshold as percentage of bulk depth (default 5 %).

Returns:

Track radius in um, or nan if no contour is found.

Return type:

float

get_track_radius_um_with_uncertainty(etch_time_h: float, threshold_percent: float = 5.0, n_sigma: float = 1.0) → tuple[float, float][source]

Track radius with uncertainty propagated from v(d) anchor uncertainties.

Evaluates the etch-rate and arrival-time maps for v(d)+n_sigma*dV and v(d)-n_sigma*dV, extracts a radius from each, and returns the central radius together with half the spread as the uncertainty estimate.

Parameters:

etch_time_h (float) – Etching duration in hours.
threshold_percent (float) – Deviation threshold as percentage of bulk depth (default 5 %).
n_sigma (float) – Number of standard deviations to propagate (default 1).

Returns:

radius_um (float) – Central track radius in um.
radius_uncertainty_um (float) – Half-spread of the upper/lower bound radii (1-sigma estimate).

Raises:

ValueError – If the etch model carries no anchor velocity uncertainties.

get_track_length_um_with_uncertainty(etch_time_h: float, relative_to_surface: bool = True, n_sigma: float = 1.0) → tuple[float, float] | None[source]

Track depth with uncertainty propagated from v(d) anchor uncertainties.

Evaluates the arrival-time map for v(d)+n_sigma*dV and v(d)-n_sigma*dV, reads the depth at r = 0 for each, and returns the central depth together with half the spread as the uncertainty estimate.

Parameters:

etch_time_h (float) – Etching duration in hours.
relative_to_surface (bool) – If True (default), subtract the bulk-etched surface depth.
n_sigma (float) – Number of standard deviations to propagate (default 1).

Returns:

tuple of (depth_um, depth_uncertainty_um) or None.
Returns *None (and logs a warning) when the etch model carries no*
anchor velocity uncertainties.

get_track_length_um(etch_time_h: float | list, relative_to_surface: bool = True) → float | ndarray[source]

Track depth at r = 0 for a given etching time.

Parameters:

etch_time_h (float or list[float]) – Etching duration(s) in hours.
relative_to_surface (bool) – If True (default), subtract the bulk-etched surface position so the result is the pit depth below the current surface.

Returns:

Depth in um, or nan if the track does not reach r = 0.

Return type:

float or ndarray

get_track_detectability(etch_time_h: float, min_radius_um: float = 0.5, min_depth_um: float = 0.3, min_cone_half_angle_deg: float = 2.0, threshold_percent: float = 5.0) → dict[source]

Geometric detectability metrics for a given etching time.

A track is considered detectable when all three criteria are met:

radius_um >= min_radius_um
depth_um >= min_depth_um (pit depth below bulk surface)
cone_half_angle_deg >= min_cone_half_angle_deg

Returns:: Keys: detected, radius_um, depth_um, cone_half_angle_deg, radius_ok, depth_ok, angle_ok.
Return type:: dict

is_track_detected(etch_time_h: float, min_radius_um: float = 0.5, min_depth_um: float = 0.3, min_cone_half_angle_deg: float = 2.0, threshold_percent: float = 5.0) → bool[source]

Return True if the track is detectable at this etching time.

Convenience wrapper around get_track_detectability().

plot_iso_time_contour(ax: Axes, etching_time_h: float) → tuple[ndarray, ndarray][source]

Plot the iso-time contour on ax.

Returns:: r_coords, z_coords
Return type:: tuple[ndarray, ndarray]

plot_track_radius(etch_time_h: float, ax: Axes | None = None, show_contour: bool = True, threshold_percent: float = 5.0) → tuple[source]

Plot iso-time contour with track-radius indicators.

Return type:: fig, ax

plot_map(name: str = 'dose', grayscale_mode: bool = False, plot_contours: bool = True, annotate_figure: bool = True) → tuple[source]: Plot the dose, etch-rate, or arrival-time map.

plot_LET_energy_profiles() → tuple[source]: Plot energy and LET profiles vs. depth.

calculate_RDD_Gy(E_MeV_u: float) → ndarray[source]

Calculate the radial dose distribution at a given energy.

Parameters:: E_MeV_u (float) – Kinetic energy in MeV per nucleon.
Returns:: Dose in Gy at each radial grid point.
Return type:: ndarray

Plotting utilities for TrackSimulator.

tracketch.simulation.plots.plot_map(etcher: tracketch.simulation.simulator.TrackSimulator, name: str = 'dose', grayscale_mode: bool = False, plot_contours: bool = True, annotate_figure: bool = True) → tuple[source]

Plot the dose, etch-rate, or arrival-time map.

Parameters:

etcher (tracketch.simulation.simulator.TrackSimulator) – Simulator instance.
name (str) – 'dose', 'etch', or 'arrival'.
grayscale_mode (bool) – Use greyscale colour map.
plot_contours (bool) – Overlay contour lines.
annotate_figure (bool) – Add CSDA-range annotation.

Returns:

fig (Figure)
ax (Axes)

tracketch.simulation.plots.plot_LET_energy_profiles(etcher: tracketch.simulation.simulator.TrackSimulator) → tuple[source]

Plot kinetic-energy and LET depth profiles.

Parameters:

etcher (tracketch.simulation.simulator.TrackSimulator) – Simulator instance.

Returns:

fig (Figure)
axes (ndarray of Axes)

Simulation grid construction and coordinate-transform helpers.

tracketch.simulation.utils.convert_map_logspace_to_linspace(r_log: ndarray, z_linear: ndarray, data_map: ndarray, n_uniform_multiplier: int = 3) → tuple[ndarray, ndarray, ndarray][source]

Interpolate a 2-D map from log-spaced r to uniform linear r.

Parameters:

r_log (ndarray) – Log-spaced radial coordinates (1-D).
z_linear (ndarray) – Linear depth coordinates (1-D).
data_map (ndarray, shape (n_z, n_r)) – Data on the (z_linear, r_log) grid.
n_uniform_multiplier (int) – Resolution multiplier for the uniform grid.

Returns:

r_uniform (ndarray)
z_linear (ndarray)
data_uniform (ndarray)

tracketch.simulation.utils.convert_map_linspace_to_logspace(r_uniform: ndarray, z_linear: ndarray, data_uniform: ndarray, r_log: ndarray) → ndarray[source]

Interpolate a 2-D map from uniform linear r back to log-spaced r.

Parameters:

r_uniform (ndarray) – Uniform radial coordinates.
z_linear (ndarray) – Depth coordinates.
data_uniform (ndarray, shape (n_z, n_r_uniform)) – Data on the uniform grid.
r_log (ndarray) – Target log-spaced radial coordinates.

Return type:

ndarray, shape (n_z, n_r_log)

tracketch.simulation.utils.create_simulation_grid(rz_lims_dict: dict) → tuple[ndarray, ndarray, tuple[float, float], tuple[float, float]][source]

Create a 2-D (r, z) simulation grid.

The radial axis is log-spaced (necessary for resolving narrow RDDs).

Parameters:

rz_lims_dict (dict) – Grid parameters. Recognised keys: r_min_um, r_max_um, z_min_um, z_max_um, n_points_r, n_points_z.

Returns:

r_grid_um (ndarray)
z_grid_um (ndarray)
r_bounds (tuple[float, float])
z_bounds (tuple[float, float])

Raises:

TypeError – If a recognised key maps to a non-numeric value.

tracketch.simulation.utils.get_track_radius_from_contour(r_contour: ndarray, z_contour: ndarray, V_bulk_um_h: float, etch_time_h: float, threshold_percent: float = 5.0) → float[source]

Find track radius where the iso-time contour deviates from bulk etching.

Parameters:

r_contour (ndarray) – Radial coordinates of iso-time contour (um).
z_contour (ndarray) – Depth coordinates of iso-time contour (um).
V_bulk_um_h (float) – Bulk etch rate in um/hr.
etch_time_h (float) – Etching time in hours.
threshold_percent (float) – Deviation threshold (percentage of bulk depth).

Returns:

Track radius in um, or nan if not found.

Return type:

float