Ray Tracer Driver

The RayTracer driver follows deterministic ray paths through the simulation domain (vs. the MCRT driver's stochastic photon transport). Use it for:

• Projections: column density maps, temperature-weighted projections
• Volume rendering: 3D visualization with transfer functions
• Line tracers: quantities along specific sight lines

Architecture

flowchart TD
    subgraph RayTracer["RayTracer Driver"]
        init[Initialize] --> run[run]
        run --> proj[run_projections]
        run --> trace[run_tracers]
        run --> abs[run_absorptiongrids]
        run --> sl[run_sightlines]
    end

    subgraph Operators["Operators"]
        proj --> ProjOp[ProjectionOperator]
        proj --> VolOp[VolumeRenderOperator]
        trace --> TracerOp[TracerOperator]
        abs --> AbsOp[OpticalDepthGridOperator]
        sl --> SlOp[SightlineOperator]
    end

    subgraph Kernel["SYCL Kernel"]
        ProjOp --> kernel[evolve_step]
        VolOp --> kernel
        TracerOp --> kernel
        AbsOp --> kernel
        SlOp --> kernel
        kernel --> field[Field Access]
    end

    field --> smart[SmartFieldAccessor]
    field --> dynamic[DynamicFieldAccessor]

Ray Tracer Parameters

These top-level keys in the raytracer: block are driver-level and apply across all operators:

Operators

ProjectionOperator

Computes 2D projections by integrating field values along rays:

Where:

• \(q(s)\) is the quantity being projected (e.g., Temperature)
• \(w(s)\) is the weight field (typically Density)
• \(L\) is the path length through the domain

Configuration example:

raytracer:
  dist_max: 0.7              # maximum ray distance (box widths)
  max_step: 1e-2             # maximum step length
  outputpath: output.zr
  overwrite: true
  linename: FeXXV            # required for line-specific fields
  populate:                   # map derived fields to loaded dataset fields
    Emissivity: "cloudyemission_Fe25 1.85040A"

  operators:
    projections:
      mode: manual
      view: orthogonal        # orthogonal, perspective, or equirectangular
      position: [0.5, 0.5, 0.0]
      direction: [0.0, 0.0, 1.0]
      up: [0.0, 1.0, 0.0]
      npixels: [512, 512]     # [width, height] in pixels
      widths: [1.0, 1.0]      # field of view [x, y] in box units
      fields:
        - Temperature
        - Density
      use_weighting: true      # weight projections by density
      perform_averaging: true  # normalize by total weight (density-weighted average)

Projection Parameters

View Types

For perspective view, an additional fov_up parameter (in degrees, default 90) controls the vertical field of view.

For healpix, a single camera position emits one ray per HEALPix pixel; the output is a 1D zarr array of length 12 * nside² (RING ordering). Required extra keys:

The angular frame is box z-up: theta is measured from +z, phi increases from +x toward +y, regardless of camera basis. The dataset carries a healpix attribute block (nside, ordering, frame, npix) in its zarr attributes; downstream Python code can call healpy.pix2ang(nside, ipix) directly on the index. Example:

projections:
  mode: manual
  view: healpix
  nside: 16             # 12 * 16² = 3072 pixels (Smith+19 convention)
  position: [0.5, 0.5, 0.5]
  start_dist: 0.005     # optional: skip inner 0.005 box-units (e.g. central SFR core)
  fields: [Density, Temperature]

stereoscopic is not supported in healpix mode.

Camera Modes

The mode parameter controls how rays are generated across frames:

Mode-specific parameters:

• rotate: requires center (point to orbit around)
• traverse: requires travel_distance (total distance to travel along direction)

VolumeRenderOperator

Performs volume rendering with density-weighted field values for 3D visualization.

TracerOperator

Records a per-ray profile along each sight line: instead of collapsing the ray to a single pixel (as ProjectionOperator does), it stores an ordered list of segments with, for each, the segment's path length and an accumulated quantity (e.g. the mean field value over that segment). Output is therefore a ragged, per-ray structure ("skewers"), not a 2D image — so it accepts any ray layout (Grid2D, Healpix, or Scatter ray sources).

The operator is built from four orthogonal, JIT-specialized choices:

1. What ends a segment — the close rule.
2. What is accumulated within a segment — the accumulator.
3. How the ray itself propagates — the propagation mode (straight, or field-aligned / curved).
4. What ends the whole ray — the stop rule (in addition to the always-on dist_max / domain-exit / stuck-ray terminators).

Segment model: close rule + accumulator

Only these (close_rule, accumulator) pairs are supported; any other combination terminates at startup:

The (never, column) pair makes the tracer a column integrator: each ray commits a single entry — the integral Σ Δs · q of each traced field from the seed to the ray's end (domain exit or raytracer.dist_max). Output naming matches weighted_mean (one contributions dataset, or contributions_<field> per field), and lengths holds the ray's total path length. A committed column is in normalized box units × field units; multiply by length_unit_cm from the output store for a physical column. This is the pairing behind the source_points ray source — one obscuration column per star (see Source points).

Curved-ray propagation

By default rays travel straight. Setting propagation_mode: field_aligned makes each ray follow a vector field (direction_field) cell-by-cell — it becomes a streamline. The ray source then only seeds the field-line entry points; the path is determined by the field. This is the same DirectionFunctor abstraction the CoherenceLengthOperator uses, so a field-aligned tracer measures, e.g., the coherence length of a magnetic field line and the mean |B| along it in one pass.

The Coherence Length cookbook page illustrates this segmentation cell by cell (schematics for both angle_reference modes) and applies it to a real snapshot.

Ray-stop rule

The close rule decides where a segment ends; the stop rule decides where the whole ray ends. They are orthogonal. By default (stop_rule: never) a ray runs until one of the always-on terminators fires — raytracer.dist_max, leaving the domain, or the stuck-ray guard. stop_rule: angle_threshold adds an early termination: the ray ends when direction_field bends past stop_angle_threshold.

stop_angle_reference chooses what the rotation is measured against, exactly like angle_reference does for the close rule:

• continuous (default) — compares each cell to the previous one, so the ray stops at the first adjacent cell-to-cell bend exceeding the threshold.
• cumulative — compares each cell to the ray's seed direction, so the ray stops once the field has drifted past the threshold relative to where it started (there is no re-seed — a stop ends the ray).

The angle is measured on the raw field direction — the stop does not apply flip_anti_parallel. A ~180° reversal of the field between cells therefore exceeds any sub-180° threshold and terminates the ray (a field reversal is treated as a genuine coherence break). This matches the close rule's angle convention; note it differs from propagation_mode: field_aligned, which does flip so the ray keeps heading forward through a reversal — so a field-aligned ray can coast through a reversal while the stop rule ends it there.

The canonical pairing is close_rule: every_cell + stop_rule: angle_threshold: it logs every cell (full per-cell fidelity) but only until the field bends, giving a bounded, physically meaningful per-ray length. That bound is what keeps a per-cell trace tractable at full resolution — an unbounded every_cell trace runs to dist_max and emits orders of magnitude more entries.

Tracer Parameters

The camera-grid keys (mode, view, position, direction, up, npixels, widths, …) are shared with the projection operators — see Projection Parameters, View Types, and Camera Modes. The tracer-specific keys:

Tracer Output

Written under the tracers/<field>/ zarr group as flat arrays indexed by a per-ray offsets table. The written arrays are compacted: ray i's entries occupy [offsets[i], offsets[i+1]) exactly (last ray runs to the end of the array), so its segment count is offsets[i+1] - offsets[i] — no padding. Naming note: the array called lengths is not a segment count: under weighted_mean it holds the per-segment path lengths (in normalized box units); under length_only it stays zero and the segment length is written to contributions instead (see the table); under column there is a single segment per ray, so it holds each ray's total path length.

Tracer MPI support

TracerOperator has MPI support for the legacy/default tracer mode: close_rule: every_cell, accumulator: weighted_mean, propagation_mode: straight, and stop_rule: never. This path is covered by the tracer_mpi, tracer_mpi_pbc, and tracer_mpi_pbc_offaxis regression tests, which compare multi-rank output against a one-rank baseline.

The newer tracer modes are single-rank only. propagation_mode: field_aligned, close_rule: angle_threshold, close_rule: never, and stop_rule: angle_threshold terminate at startup under MPI: field-aligned rays bend away from the precomputed per-rank segment table, angle-threshold segmentation would over-split at every rank boundary, a never-close column would fragment across ranks instead of committing once per ray, and the angle stop's per-ray reference cannot follow a ray across a rank handoff. Scatter ray sources (e.g. random_uniform) are likewise single-rank because they have no meaningful camera depth axis to merge.

Tracer Examples

Classic per-cell skewer (default), Density along an orthographic grid:

raytracer:
  operators:
    tracers:
      mode: manual
      view: orthogonal
      position: [0.01, 0.5, 0.5]
      direction: [1.0, 0.0, 0.0]
      npixels: [256, 256]
      widths: [0.98, 0.98]
      fields: [Density]
      # close_rule: every_cell, accumulator: weighted_mean, weight: density (defaults)

Magnetic field-line coherence length and mean |B| per segment (curved ray):

raytracer:
  dist_max: 1.1          # REQUIRED for field_aligned; max path length to trace
  operators:
    tracers:
      mode: manual
      view: orthogonal
      position: [0.01, 0.5, 0.5]
      direction: [1.0, 0.0, 0.0]   # seeds entry pixels only; path follows B
      npixels: [640, 640]
      widths: [0.98, 0.98]
      fields: [MagneticFieldMagnitude]
      close_rule: angle_threshold
      accumulator: weighted_mean
      weight: none                 # path-length-weighted (spatial) mean |B|
      angle_threshold: 90.0
      min_segment_length: 2.71e-8  # ~1 pc: drop sub-resolution micro-segments
      propagation_mode: field_aligned
      direction_field: MagneticField
      flip_anti_parallel: true     # required for magnetic (sign-agnostic) fields
      min_field_magnitude: 1.0e-20

SightlineOperator

Traces a user-supplied list of straight rays — each defined by an explicit start and end point — and dumps the state of every cell each ray intersects. The output schema mirrors Trident's lightray.h5 model: one entry per traversed cell, carrying the path length dl plus per-cell field values. Always-on per-cell fields are density, temperature, and the three velocity components (velocity_x/y/z); extra fields (e.g. ion densities) are requested via fields:.

Optionally, a tau: block enables in-engine optical-depth integration: after the cell-field passes, a Voigt-deposition kernel computes a per-ray optical-depth spectrum τ(λ) directly from the collected per-cell data — no extra mesh queries:

Where:

• \(\rho_c\) is the (per-line) density in cell \(c\), \(dl_c\) the path length through it
• \(\sigma_{\text{Voigt}}\) is the Voigt profile cross-section, with thermal + (optional) turbulent broadening
• \(v_{\text{LOS},c}\) is the line-of-sight velocity, including peculiar motion and optional Hubble flow

Configuration example (cell dump):

raytracer:
  outputpath: output.zr
  overwrite: true

  operators:
    sightlines:
      rays:
        - start: [0.0, 0.5, 0.5]    # box / code units
          end:   [1.0, 0.5, 0.5]
          id: 0                       # optional; defaults to list index
        - start: [0.5, 0.0, 0.5]
          end:   [0.5, 1.0, 0.5]
      fields:                         # extra per-cell fields (optional)
        - HIDensity
      output:
        cells: true                   # write per-cell datasets (default true)

Sightline Parameters

Generated rays (ray_source:):

Instead of listing rays explicitly, you can generate them programmatically with the shared ray_source: factory (the same one the other raytracer operators use). This is the way to produce large ray sets — a grid of parallel sightlines for tomography, an all-sky set from an observer, or a Monte-Carlo sample. Because the generators yield a (position, direction) per ray but no endpoint, sightlines attach a fixed length: (box units): each ray runs from its generated position to position + length · direction.

    sightlines:
      ray_source:
        type: camera_plane       # grid of sightlines across a plane
        npixels: [256, 256]      # -> 65536 rays
        position: [0.5, 0.5, 0.0]
        direction: [0.0, 0.0, 1.0]
        up: [0.0, 1.0, 0.0]
        widths: [1.0, 1.0]
        length: 1.0              # ray length in box units (required)
      # tau:, fields:, output: all work the same as with an explicit rays: list

The per-type knobs (npixels/widths/position/direction/up for camera_plane; nside for camera_healpix; n_rays/direction_distribution/seed for random_uniform) match the Projection Parameters and the other operators' ray_source: blocks. rays: and ray_source: are mutually exclusive — provide exactly one. Generated rays are auto-assigned sequential ids (their index).

Optical-depth integration (tau:):

Add a tau: block to compute τ(λ) spectra. Two mutually exclusive modes, auto-detected by the presence of a lines: sequence:

    sightlines:
      rays: [...]
      tau:
        # --- single-line mode (window relative to line centre) ---
        linename: OVI                 # line identifier
        edge_left: -5.0               # Δλ (Å) blue edge from line centre
        edge_right: 5.0               # Δλ (Å) red edge
        nbins: 200                    # number of wavelength bins

        # --- OR multi-line mode (absolute window) ---
        # lines: [MgII, MgII_2]         # e.g. the MgII 2796/2803 doublet
        # lambda_min: 2790.0            # absolute Å
        # lambda_max: 2810.0
        # nbins: 400

        include_lyman_continuum: false  # add HI Lyman-continuum opacity
        use_vturb: false                # add turbulent broadening (TurbVelocity)
        line_density_fields:            # optional per-line density routing
          OVI: density                  # map a line to its per-cell density field
        ll_density_field: HIDensity     # density field for the LL term

Tau Parameters

gadget:                 # (or unigrid:, …) — the dataset loader computes the ion fields
  fields: [..., cloudydensity_HI, cloudydensity_CIV, cloudydensity_OVI]
...
sightlines:
  tau:
    lines: ["Lya", "CIV 1548", "OVI 1032"]
    line_density_fields:           # auto-collected — no separate sightlines.fields needed
      "Lya": cloudydensity_HI
      "CIV 1548": cloudydensity_CIV
      "OVI 1032": cloudydensity_OVI

Per-species τ (per_species: true). By default all lines are summed into a single tau. With per_species: true, each density-field group (an ion — lines sharing a density field, so a doublet like MgII/MgII_2 stays together) also gets its own tau_<label> dataset, labelled by the group's first line (spaces in NIST labels become underscores, e.g. CIV 1548 → tau_CIV_1548). The combined tau is still written. Unlike Trident — which emits only the summed τ and needs a separate run per species — THOR produces all per-species spectra in one pass.

Each tau_<label> is stored compactly: only over that ion's wavelength window, as a contiguous slice of the shared lam grid. The window is data-driven — the deposited τ already encodes peculiar velocity, Hubble flow (per hubble_mode), and thermal+turbulent broadening, so the window is just its support: the bins where the ion's τ (max over rays) exceeds per_species_window_rtol × peak (default 1e-6), always including the line centres, with a small pad. (A group with no absorption falls back to its lines ± per_species_window_kms.) Placement is recorded in the group attr per_species[<dataset>] = {bin_offset, nbins, lam_min, lam_max}, so a consumer stitches it back with no interpolation: full[bin_offset : bin_offset+nbins] = tau_<label>.

This adapts per ion (a narrow line gets a tight window; a doublet or high-velocity absorber a wider one) and conserves the combined τ's equivalent width to the wing level set by per_species_window_rtol; at the default rtol of 1e-6 the measured EW closure is ≈1e-5. On the OVI/CIV/SiIV AGORA run it is ~1.9× tighter than a fixed ±1500 km/s window and ~2% of a full-grid copy.

    tau:
      lines: [Lya, MgII, MgII_2, OVI]
      per_species: true
      line_density_fields: { Lya: HIDensity, MgII: MgIIDensity, MgII_2: MgIIDensity, OVI: density }
    # writes: tau (combined) + tau_Lya, tau_MgII (doublet summed), tau_OVI

The Hubble flow and cosmology used by the τ kernel are read from the driver-level raytracer: block (not under operators):

Stages

The operator is built in stages. Both are available; the difference is how many MPI ranks the run uses — the YAML and outputs are identical either way.

• Stage 1 — single-rank. Each ray is traced on one rank and its intersected cells emitted directly. Covers the full feature set: inline rays: or generated ray_source:, the always-on core fields + extra fields:, optional tau: integration, and output.cells. v_LOS / redshift_eff are deferred to the Python lightray converter, not written here.

• Stage 2 — MPI multi-rank. When run on >1 rank, each ray is split into per-rank segments (compute_ray_rank_path with per-ray dist_max = |end − start|); τ is integrated on each rank's local per-cell data, then distributed_merge gathers the per-rank dl/field arrays to rank 0 and reconstructs every ray in depth order along its own direction (zero-dl boundary cells are dropped so the MPI trace matches the single-rank one). Rank 0 writes the output.

Ray Sources

Most raytracer examples put camera keys (mode, view, position, direction, npixels, ...) directly under an operator block. Operators can also use a ray_source: block, which selects a named ray generator and makes the ray layout explicit. (SightlineOperator uses this same factory — see its Generated rays note for the sightline-specific length:.)

Supported ray_source.type values:

The camera_* sources are the ray_source: spelling of the manual camera path: type determines the camera view, so mode/view are not needed inside the block. Camera animation modes (rotate, traverse, nframes, and explicit frame lists) remain on the legacy inline camera configuration. stereoscopic is rejected under ray_source: because non-image operators size their buffers from the mono ray layout.

Random uniform rays

random_uniform emits n_rays rays with start positions

x, y, z ~ Uniform([0, 1))

in normalized box coordinates. It sets each ray origin to the sampled position. The output layout is Scatter, so rays have no pixel, image, or HEALPix ordering.

Example with random positions and isotropic random directions:

raytracer:
  operators:
    tracers:
      ray_source:
        type: random_uniform
        n_rays: 10000
        direction_distribution: isotropic
        seed: 12345
      fields: [Density]

Example with random positions and a shared fixed direction:

raytracer:
  operators:
    tracers:
      ray_source:
        type: random_uniform
        n_rays: 10000
        direction_distribution: fixed
        direction: [0.0, 0.0, 1.0]
        seed: 12345
      fields: [Density]

Uniform-grid cell centers

grid_centers puts one ray at the center of each cell of a uniform N³ (or [nx, ny, nz]) lattice — a repeatable, volume-uniform seeding of the box.

Dataset point centers

dataset_points seeds one ray per dataset cell — for a Voronoi mesh, the cell-generating sites — or a random count subsample (0 = all, reproducible from seed). Where the grid is volume-uniform, this samples per cell, so rays concentrate where the cells are smallest, i.e. in the dense gas. Needs a dataset with per-cell positions (Voronoi/particle); others report no points and stop.

Scatter layout is currently useful for TracerOperator on one MPI rank. Multi-rank tracer runs reject Scatter sources because distributed skewer merging needs a meaningful camera depth order. Structured-output operators also reject Scatter: ProjectionOperator requires Grid2D or Healpix, while VolumeRenderOperator, OpticalDepthGridOperator, and CoherenceLengthOperator require Grid2D.

Source points (stellar particles)

source_points seeds one ray per stellar particle, loading the particle list through the same source machinery the MCRT driver uses. The sources: block is a self-contained listsource-shaped body — a loader: key (gadget, tipsy_stars, ramses_part, art_stars, enzo_particles) plus that loader's config section, with the same schema as emissionmodel.listsource (see Configuration → List Source). The gadget loader requires a luminosity_field; the others default it to the Lya_StellarParticle_Luminosity recipe, so set it explicitly to be sure. Stars with non-positive luminosity are dropped before any ray is created.

If the gas dataset uses a zoom_box, the sources: block inherits it automatically so star seeds land in the same frame; set zoom_box inside sources: to override.

Ray order, emitter order, and output row order all match. There is no count subsample (a count key is rejected) — it would break the row ↔ star identity. To keep rows self-identifying, the driver also writes these datasets into the tracers group, in emitter order:

(The tracer's own origins array records where each ray ended up after tracing, not the seeds — use source_positions for the stars.)

source_particle_ids is the stable join key for post-processing — it ties each column back to its snapshot particle regardless of order or filtering, e.g. to write a ParticleID-keyed sidecar field for THOR to re-load and project. It is present only when the loader supplies real IDs (the gadget loader reads PartType<X>/ParticleIDs).

source_points is single-rank only (an MPI run terminates at startup), and pbc should stay at its default false — a column runs to the domain boundary.

The intended use is per-star obscuration for visualization: pair it with close_rule: never + accumulator: column, and point the rays at the observer (a fixed direction of minus the camera direction), so each ray commits one value — the density column between the star and the domain boundary:

raytracer:
  outputpath: star_columns.zr
  overwrite: true
  # pbc defaults to false: each column runs star -> domain boundary.
  # pbc: true would require a finite dist_max.

  operators:
    tracers:
      ray_source:
        type: source_points
        direction_distribution: fixed  # default; isotropic also allowed
        direction: [0.0, 0.0, -1.0]    # toward the observer = minus the camera direction
        sources:                       # listsource-shaped body
          loader: gadget
          gadget:
            path: ./snap_099
            particle_type: "PartType4"
            luminosity_field: Masses   # required; lum <= 0 emitters are dropped
      fields: [Density]
      close_rule: never
      accumulator: column

Recovering per-star columns. A star outside the mesh commits no entry, so read the columns through offsets (row i has offsets[i+1] - offsets[i] entries — 0 or 1 here):

import numpy as np
import zarr

store = zarr.open("star_columns.zr", mode="r")
grp = store["tracers"]

offsets = grp["density/offsets"][:]
contribs = grp["density/contributions"][:]           # 0 or 1 entry per star
counts = np.diff(np.append(offsets, len(contribs)))  # 0 = no column for this star

columns = np.full(offsets.size, np.nan)              # NaN where no column
columns[counts == 1] = contribs[offsets[counts == 1]]

positions = grp["source_positions"][:]               # (n, 3), same row order
luminosities = grp["source_luminosities"][:]         # (n,)

# physical column: box-unit path length -> cm
length_unit_cm = grp["density/contributions"].attrs["Parameters"]["length_unit_cm"]
columns_physical = columns * length_unit_cm          # field units x cm

Field Access

Fields are selected at runtime; see JIT Field Access for the dispatch mechanism.

Output

Results are written to a Zarr store (.zr directory):

• Projections: 2D arrays of shape (npixels_y, npixels_x)
• Volume renders: RGB(A) images
• Tracers: per-ray field profiles
• Sightlines: a single flat group sightlines/default/ (one dataset per quantity, indexed by ray — not one group per ray), containing:
  • ray_id (nray,), start/end (nray, 3) — per-ray metadata
  • cell_offsets (nray+1,) — CSR index; ray i's cells are [cell_offsets[i], cell_offsets[i+1]). Written when output.cells: true
  • dl and one array per emitted field (density, temperature, velocity_x/y/z, extras) — flat per-cell arrays (length = total cells); slice with cell_offsets. Written when output.cells: true
  • tau (nray, nbins) — integrated optical-depth spectra; written when tau: is configured
  • tau_<line> (nray, window_bins) — per-species spectra, compact (each on its ion's own window, a slice of lam); written when tau.per_species: true. Placement is in the group attr per_species[<dataset>] = {bin_offset, nbins, lam_min, lam_max}; stitch with full[bin_offset:bin_offset+nbins] = tau_<line> (no interpolation)
  • lam (nbins,) — bin-centred wavelengths in Å, written once (shared by all rays)
  • group attrs: layout (scatter/grid2d/healpix), nrays, nbins, and grid_shape/nside for structured ray sources

For a grid2d (camera_plane) source, tau reshapes to (ny, nx, nbins) via grid_shape; for healpix, to (npix, nbins) via nside. The per-cell data converts to a Trident-compatible LightRay HDF5 via thor.absorption.sightline.Sightline (slices each ray by cell_offsets, computes velocity_los/redshift_eff); integrated spectra load with thor.absorption.spectrum.Spectrum.from_sightline.

See Raytracer Postprocessing for reading the output.