Backtrace (`data.backtrace`) — Experimental¶

data.backtrace is the entry point for integrating particle trajectories backwards in time through EMSES fields. You can compute arrival probabilities (get_probabilities) or trajectories for individual particles (get_backtrace / get_backtraces). Results come back wrapped in dedicated containers that let you chain straight into visualization with shorthand like .vxvz.plot().

Requirements: backtrace relies on the external vdist-solver-fortran package (vdsolverf). Install it with pip install vdist-solver-fortran. Without it, calls to data.backtrace.* raise ImportError.

When to use it¶

You want the phase-space distribution of particles that arrive at a given observation point.
You want to trace back the trajectory of a particle of interest to see where it came from.
You want to draw an energy spectrum of arriving particles.

Backtrace integrates an ODE backwards using data.inp.dt and the saved EMSES fields, so a large max_step can become expensive. If you want to push the work to an HPC node, combine it with the remote-execution backend (see below).

Quick start¶

import emout

data = emout.Emout("output_dir")

# Single particle
position = (20.0, 32.0, 40.0)
velocity = (1.0e5, 0.0, -2.0e5)
bt = data.backtrace.get_backtrace(position, velocity, ispec=0)

bt.tx.plot()      # t vs x trajectory
bt.xvz.plot()     # x vs vz phase space

# Many particles in one call
import numpy as np
positions = np.array([[20, 32, 40], [21, 32, 40], [22, 32, 40]], dtype=float)
velocities = np.zeros_like(positions)
velocities[:, 0] = 1.0e5
many = data.backtrace.get_backtraces(positions, velocities, ispec=0)

many.xz.plot(alpha=0.5)    # overlay all trajectories

# Arrival probability over a 6-D phase-space grid
result = data.backtrace.get_probabilities(
    x=20.0, y=32.0, z=40.0,
    vx=(-3e5, 3e5, 64),
    vy=0.0,
    vz=(-3e5, 3e5, 64),
    ispec=0,
)

result.vxvz.plot(cmap="viridis")      # heatmap in the vx-vz plane
result.plot_energy_spectrum(scale="log")

Single particle: `get_backtrace`¶

get_backtrace(position, velocity, ispec=0, ...) integrates one trajectory and returns a :class:BacktraceResult. The result also supports tuple unpacking.

bt = data.backtrace.get_backtrace(position, velocity, ispec=0, max_step=50000)

ts, prob, positions, velocities = bt   # tuple unpacking
print(bt)                                # <BacktraceResult: n_steps=...>

Attribute	Shape	Meaning
`bt.ts`	`(N,)`	time (EMSES units)
`bt.probability`	`(N,)`	arrival probability per step
`bt.positions`	`(N, 3)`	`[x, y, z]`
`bt.velocities`	`(N, 3)`	`[vx, vy, vz]`

Shorthand plotting¶

bt.pair(var1, var2) selects two variables and returns an :class:XYData. var1 and var2 can each be t, x, y, z, vx, vy, or vz. The concatenated form (bt.tx, bt.xvz, bt.yz, …) is shorthand for the same call.

bt.tx.plot()                 # = bt.pair("t", "x")
bt.xvz.plot()                # = bt.pair("x", "vz")
bt.yz.plot(color="black")    # xy projection of the trajectory

XYData.plot() converts to SI units by default and auto-generates axis labels (use_si=False keeps EMSES units). Passing gap=... inserts NaN breaks where consecutive points are too far apart, which is handy to avoid spurious lines across periodic-boundary jumps.

Many particles: `get_backtraces`¶

get_backtraces(positions, velocities, ispec=0, n_threads=4, ...) returns a :class:MultiBacktraceResult. positions and velocities must be (N, 3) arrays.

ts, probs, pos_list, vel_list, last = many
many.xz.plot(alpha=np.clip(probs, 0, 1))    # alpha weighted by probability
many.sample(50, random_state=0).tvx.plot()   # random 50 trajectories
many.sample(slice(0, 10)).tx.plot()         # first 10 trajectories

Attribute	Shape	Meaning
`ts_list`	`(N_traj, N_steps)`
`probabilities`	`(N_traj,)`	final arrival probability
`positions_list`	`(N_traj, N_steps, 3)`
`velocities_list`	`(N_traj, N_steps, 3)`
`last_indexes`	`(N_traj,)`	end of valid data for each trajectory (padding)

many.pair("t", "x") returns a :class:MultiXYData; .plot() overlays every trajectory. alpha accepts either a scalar or an array of length N_traj, which is convenient for probability weighting.

Feeding raw Particle objects¶

If you already have vdsolverf.core.Particle instances (for example from ProbabilityResult.particles), use get_backtraces_from_particles(particles, ...):

from vdsolverf.core import Particle

particles = [Particle(p, v) for p, v in zip(positions, velocities)]
many = data.backtrace.get_backtraces_from_particles(particles, ispec=0)

A common pattern is to chain this with get_probabilities — compute the probability grid, then trace back only the particles you care about:

result = data.backtrace.get_probabilities(...)
bt = data.backtrace.get_backtraces_from_particles(result.particles, ispec=0)
bt.xz.plot(alpha=np.clip(result.probabilities, 0, 1))

Arrival probability: `get_probabilities`¶

get_probabilities(x, y, z, vx, vy, vz, ispec=0, ...) builds a 6-D phase-space grid, backtraces a particle from every grid point, and returns the arrival probabilities as a :class:ProbabilityResult.

Each axis accepts:

a tuple (start, stop, n) — equally-spaced grid
an explicit array or list — arbitrary values
a scalar — a size-1 axis (automatically squeezed when you call pair())

result = data.backtrace.get_probabilities(
    x=20.0, y=32.0, z=40.0,     # fixed position
    vx=(-3e5, 3e5, 64),         # scan vx over 64 points
    vy=0.0,
    vz=(-3e5, 3e5, 64),         # scan vz over 64 points
    ispec=0,
    max_step=10000,
    n_threads=8,
)

2-D heatmap projections¶

result.pair(var1, var2) integrates out the four unselected axes (trapezoidal rule) and returns a :class:HeatmapData. The shorthand attribute form works the same way as for BacktraceResult.

result.vxvz.plot(cmap="viridis")   # = result.pair("vx", "vz")
result.xvx.plot()                  # x-vx plane
result.yz.plot(cmap="plasma")      # y-z plane

HeatmapData.plot() draws a pcolormesh with a colour bar and SI-unit labels (use_si=False keeps grid units). Extra keyword arguments are forwarded straight to pcolormesh, so you can use vmin / vmax or norm=LogNorm(...) to control the colour scale. Pass offsets=("center", 0) to centre an axis or apply a numeric shift.

Energy spectrum¶

plot_energy_spectrum(energy_bins=None, scale="log") renders an energy-flux histogram of the arriving particles. energy_bins accepts either an integer (number of bins) or an array of bin edges.

result.plot_energy_spectrum(scale="log", energy_bins=80)

Internally it reads wp (or the photoelectron settings path / curf when nflag_emit == 2) from plasma.inp to compute a reference number density n0, weights each phase-space point by its probability, and integrates.

Raw histogram arrays¶

result.energy_spectrum(energy_bins=...) returns the (hist, bin_edges) tuple directly so you can feed it to custom post-processing or another library (e.g. ax.step, seaborn).

Remote execution integration¶

data.backtrace shares the Emout facade’s remote_open_kwargs, so if an emout server is running the computation automatically runs on the worker and you get back a RemoteProbabilityResult / RemoteBacktraceResult proxy. Because the result is cached on the worker, changing visualisation parameters does not trigger recomputation.

from emout.distributed import remote_figure

result = data.backtrace.get_probabilities(...)   # computed once on the worker

with remote_figure():
    result.vxvz.plot(cmap="viridis")

with remote_figure():
    result.plot_energy_spectrum(scale="log")

result.drop()   # free worker memory when done

If you prefer the explicit remote style, switch to data.remote().backtrace... — it returns the same dedicated proxies:

from emout.distributed import remote_scope, remote_figure

with remote_scope():
    rdata = data.remote()
    bt = rdata.backtrace.get_backtrace(position, velocity, ispec=0)
    result = rdata.backtrace.get_probabilities(...)

    with remote_figure():
        bt.tx.plot()
        result.vxvz.plot()

For the remote-execution mechanics, environment variables, and server management, see the remote execution guide.

`fetch()` for local customisation¶

When you want full matplotlib control (custom annotations, shared colour bars, dropping the heatmap into your own subplot grid), use fetch() to pull the small result arrays back to the client:

heatmap = result.vxvz.fetch()      # local HeatmapData
fig, ax = plt.subplots()
heatmap.plot(ax=ax, cmap="plasma")
ax.axhline(y=0, color="red", linestyle="--")

Backtrace (data.backtrace) — Experimental¶