API Reference

ViSAGE exposes a Python API for scripting, integration with notebooks, and extending the viewer.

IO layer

visage.io.par_reader

visage.io.par_reader.parse_par(par_path)

Parse a SAGE .par file into a SimConfig.

Handles both absolute and relative paths. Relative paths in the par file are resolved relative to the par file's parent directory (the SAGE root). On Windows, /mnt/c/... paths (written for SAGE running inside WSL) are read as C:... .

Source code in visage/io/par_reader.py

def parse_par(par_path: str | Path) -> SimConfig:
    """Parse a SAGE .par file into a SimConfig.

    Handles both absolute and relative paths. Relative paths in the par file
    are resolved relative to the par file's parent directory (the SAGE root).
    On Windows, /mnt/c/... paths (written for SAGE running inside WSL) are
    read as C:\\... .
    """
    par_path = Path(par_path).resolve()
    # Paths in the par file are relative to the SAGE root (parent of input/)
    root = par_path.parent.parent

    raw: dict[str, str] = {}
    with open(par_path, encoding="utf-8", errors="replace") as f:
        for line in f:
            line = line.strip()
            # Strip inline comments — SAGE par files use either '%' (old
            # style) or ';' (newer style) as the comment marker.
            for marker in ("%", ";", "#"):
                if marker in line:
                    line = line[: line.index(marker)].strip()
            if not line:
                continue
            # Skip snapshot list arrow
            if line.startswith("->"):
                continue
            parts = line.split(None, 1)
            if len(parts) == 2:
                raw[parts[0]] = parts[1].strip()

    def _path(key: str, default: str) -> Path:
        p = native_path(raw.get(key, default))
        return p if p.is_absolute() else root / p

    def _int(key: str, default: int) -> int:
        return int(raw.get(key, default))

    def _float(key: str, default: float) -> float:
        return float(raw.get(key, default))

    def _str(key: str, default: str) -> str:
        return raw.get(key, default)

    known_keys = {
        "OutputDir",
        "FileNameGalaxies",
        "FirstFile",
        "LastFile",
        "OutputFormat",
        "TreeName",
        "TreeType",
        "SimulationDir",
        "FileWithSnapList",
        "LastSnapShotNr",
        "NumSimulationTreeFiles",
        "Omega",
        "OmegaLambda",
        "Hubble_h",
        "BoxSize",
        "PartMass",
        "NumOutputs",
    }
    extra = {k: v for k, v in raw.items() if k not in known_keys}

    snap_list_default = "input/millennium/trees/millennium.a_list"

    cfg = SimConfig(
        par_path=par_path,
        output_dir=_path("OutputDir", "output"),
        file_name_galaxies=_str("FileNameGalaxies", "model"),
        first_file=_int("FirstFile", 0),
        last_file=_int("LastFile", 7),
        output_format=_str("OutputFormat", "sage_hdf5"),
        tree_name=_str("TreeName", "trees_063"),
        tree_type=_str("TreeType", "lhalo_binary"),
        simulation_dir=_path("SimulationDir", "input/millennium/trees"),
        snap_list_path=_path("FileWithSnapList", snap_list_default),
        last_snapshot_nr=_int("LastSnapShotNr", 63),
        num_sim_tree_files=_int("NumSimulationTreeFiles", 8),
        omega=_float("Omega", 0.25),
        omega_lambda=_float("OmegaLambda", 0.75),
        hubble_h=_float("Hubble_h", 0.73),
        box_size=_float("BoxSize", 62.5),
        part_mass=_float("PartMass", 0.086),
        extra=extra,
    )
    return cfg

visage.io.snapshot_table

visage.io.snapshot_table.SnapshotTable

Maps snapshot indices ↔ redshifts ↔ scale factors.

Reads a plain-text file of scale factor values (one per line), such as millennium.a_list or Uchuu100_scalefactor.txt.

Source code in visage/io/snapshot_table.py

class SnapshotTable:
    """Maps snapshot indices ↔ redshifts ↔ scale factors.

    Reads a plain-text file of scale factor values (one per line), such as
    millennium.a_list or Uchuu100_scalefactor.txt.
    """

    def __init__(self, a_list_path: str | Path) -> None:
        self._path = Path(a_list_path)
        self._a = self._load(self._path)
        self._z = 1.0 / self._a - 1.0

    @staticmethod
    def _load(path: Path) -> np.ndarray:
        values = []
        with open(path) as f:
            for line in f:
                line = line.strip()
                if not line or line.startswith("#"):
                    continue
                # Some files have multiple whitespace-separated values per line
                for token in line.split():
                    try:
                        values.append(float(token))
                    except ValueError:
                        pass
        return np.array(values, dtype=np.float64)

    @property
    def count(self) -> int:
        return len(self._a)

    @property
    def scale_factors(self) -> np.ndarray:
        return self._a.copy()

    @property
    def redshifts(self) -> np.ndarray:
        return self._z.copy()

    def snap_to_a(self, snap: int) -> float:
        return float(self._a[snap])

    def snap_to_z(self, snap: int) -> float:
        return float(self._z[snap])

    def z_to_snap(self, z: float) -> int:
        """Return the snapshot index closest to redshift z."""
        return int(np.argmin(np.abs(self._z - z)))

    def a_to_snap(self, a: float) -> int:
        """Return the snapshot index closest to scale factor a."""
        return int(np.argmin(np.abs(self._a - a)))

    def label(self, snap: int) -> str:
        z = self.snap_to_z(snap)
        a = self.snap_to_a(snap)
        return f"Snap {snap:02d}  |  z = {z:.2f}  |  a = {a:.4f}"

a_to_snap(a)

Return the snapshot index closest to scale factor a.

Source code in visage/io/snapshot_table.py

def a_to_snap(self, a: float) -> int:
    """Return the snapshot index closest to scale factor a."""
    return int(np.argmin(np.abs(self._a - a)))

z_to_snap(z)

Return the snapshot index closest to redshift z.

Source code in visage/io/snapshot_table.py

def z_to_snap(self, z: float) -> int:
    """Return the snapshot index closest to redshift z."""
    return int(np.argmin(np.abs(self._z - z)))

visage.io.halo_reader

visage.io.halo_reader.load_halo_snapshot(tree_dir, tree_name, snap_num, first_file=0, last_file=7, mass_cut=10000000000.0, max_halos=None, hubble_h=0.73, n_jobs=-1, box_size=0.0, tree_type='')

Load halo positions and masses for one snapshot from lhalo_binary tree files.

Parameters:

Name	Type	Description	Default
tree_dir	str | Path		required
tree_name	str		required
snap_num	int		required
mass_cut	float		10000000000.0
max_halos	int | None	loads every halo — no thinning of the field.	None
n_jobs	int		-1
The			required
sets			required
load			required

Source code in visage/io/halo_reader.py

def load_halo_snapshot(
    tree_dir: str | Path,
    tree_name: str,
    snap_num: int,
    first_file: int = 0,
    last_file: int = 7,
    mass_cut: float = 1.0e10,
    max_halos: int | None = None,
    hubble_h: float = 0.73,
    n_jobs: int = -1,
    box_size: float = 0.0,
    tree_type: str = "",
) -> HaloSnapshot:
    """Load halo positions and masses for one snapshot from lhalo_binary tree files.

    Parameters
    ----------
    tree_dir:  directory containing the tree files
    tree_name: base name (e.g. 'trees_063'); files are tree_name.{first_file..last_file}
    snap_num:  snapshot index to extract
    mass_cut:  minimum halo mass in Msun (after h correction)
    max_halos: optional random downsample ceiling.  None (the default)
               loads every halo — no thinning of the field.
    n_jobs:    joblib parallel workers (-1 = all CPUs)

    The mass is read from the first populated column of MASS_FIELDS, so tree
    sets that leave Mvir at zero and carry the mass in M_TopHat / M_Mean200
    load correctly instead of coming back empty.
    """
    tree_dir = Path(tree_dir)
    tree_files = _resolve_tree_files(
        tree_dir, tree_name, first_file, last_file, tree_type
    )
    if not tree_files:
        # Worth its own message: an unreadable/misnamed tree path used to look
        # exactly like a mass cut that filtered everything out.
        if VERBOSE:
            print(
                f"  Haloes: no readable tree files for "
                f"{tree_dir / tree_name}"
                + (f"  (TreeType {tree_type})" if tree_type else "")
            )
            if "ascii" in tree_type.lower():
                # Say why, rather than leaving an empty box to puzzle over.
                print(
                    "  Haloes: ASCII merger trees are not read directly — "
                    "ViSAGE needs the lhalo_binary (or HDF5) trees.  Point "
                    "SimulationDir at the converted set, or keep it beside "
                    "the ASCII file and it will be found automatically."
                )
        return HaloSnapshot.empty(snap_num)

    n_files = len(tree_files)
    if VERBOSE:
        print(
            f"  Haloes: reading {n_files} tree file(s) in parallel (snap {snap_num})..."
        )

    hdf5_files = [tf for tf in tree_files if h5py.is_hdf5(tf)]
    if hdf5_files:
        # HDF5 trees (SAGE's lhalo_hdf5): read serially — h5py serialises the
        # per-tree reads anyway, and each file is scanned once then cached.
        from visage.io import hdf5_tree_reader

        results = [
            hdf5_tree_reader.read_snapshot(
                tf, snap_num, mass_cut, hubble_h, box_size
            )
            for tf in hdf5_files
        ]
    else:
        # prefer="threads": file I/O releases the GIL so threads are fully parallel
        # and avoid the semaphore / mmap leak that loky process pools produce
        results = Parallel(n_jobs=n_jobs, prefer="threads")(
            delayed(_read_tree_file)(
                tf, snap_num, mass_cut, hubble_h, box_size
            )
            for tf in tree_files
        )

    results = [r for r in results if len(r[0]) > 0]
    if not results:
        if VERBOSE:
            print(f"  Haloes: none found above mass cut ({mass_cut:.1e} Msun)")
        return HaloSnapshot.empty(snap_num)

    positions = np.vstack([r[0] for r in results])
    masses = np.concatenate([r[1] for r in results])
    vmax = np.concatenate([r[2] for r in results])
    rvir = np.concatenate([r[3] for r in results])
    vvir = np.concatenate([r[4] for r in results])

    if max_halos is not None and len(positions) > max_halos:
        rng = np.random.default_rng(42)
        idx = rng.choice(len(positions), max_halos, replace=False)
        positions, masses, vmax, rvir, vvir = (
            positions[idx],
            masses[idx],
            vmax[idx],
            rvir[idx],
            vvir[idx],
        )

    if VERBOSE:
        # Say so when the trees don't carry the mass in Mvir, since it changes
        # what "halo mass" means for colouring, filtering and Rvir/Vvir.
        used = sorted({r[5] for r in results} - {"", MASS_FIELDS[0]})
        note = f" (halo mass from {', '.join(used)})" if used else ""
        print(f"  Haloes: {len(positions):,} loaded{note}")
    return HaloSnapshot(
        positions=positions,
        masses=masses,
        vmax=vmax,
        rvir=rvir,
        vvir=vvir,
        snap_num=snap_num,
    )

visage.io.halo_reader.HaloSnapshot dataclass

Source code in visage/io/halo_reader.py

@dataclass
class HaloSnapshot:
    positions: np.ndarray  # (N, 3) float32, Mpc/h
    masses: np.ndarray  # (N,)   float32, Msun
    vmax: np.ndarray  # (N,)   float32, km/s
    rvir: np.ndarray  # (N,)   float32, Mpc/h  (computed from Mvir)
    vvir: np.ndarray  # (N,)   float32, km/s   (computed from Rvir)
    snap_num: int

    @property
    def count(self) -> int:
        return len(self.positions)

    @classmethod
    def empty(cls, snap_num: int) -> HaloSnapshot:
        z = np.empty(0, dtype=np.float32)
        return cls(
            positions=np.empty((0, 3), dtype=np.float32),
            masses=z,
            vmax=z,
            rvir=z,
            vvir=z,
            snap_num=snap_num,
        )

visage.io.galaxy_reader

visage.io.galaxy_reader.load_galaxy_snapshot(hdf5_path, snap_num, min_stellar_mass=100000000.0, max_galaxies=None, hubble_h=None, scale_factors=None, omega_m=None, omega_l=None, min_particles=None, particle_mass=None)

Source code in visage/io/galaxy_reader.py

def load_galaxy_snapshot(
    hdf5_path: str | Path,
    snap_num: int,
    min_stellar_mass: float = 1.0e8,
    max_galaxies: int | None = None,
    hubble_h: float | None = None,
    scale_factors: np.ndarray | None = None,
    omega_m: float | None = None,
    omega_l: float | None = None,
    min_particles: int | None = None,
    particle_mass: float | None = None,
) -> GalaxySnapshot:
    # Prefer the cosmology + scale factors that are written into the HDF5
    # header.  Arguments are now overrides / fallbacks only.
    from visage.io.sage_header import read_sage_header

    hdr = read_sage_header(hdf5_path)
    if hubble_h is None:
        hubble_h = hdr.get("hubble_h", 0.73)
    if omega_m is None:
        omega_m = hdr.get("omega_m", 0.315)
    if omega_l is None:
        omega_l = hdr.get("omega_l", 0.685)
    if scale_factors is None:
        scale_factors = hdr.get("scale_factors", None)
    """Load galaxy positions and properties for one snapshot from SAGE HDF5 output.

    Parameters
    ----------
    hdf5_path:        path to SAGE HDF5 output (model_0.hdf5; an MPI run's
                      other rank files model_1.hdf5 … are read with it)
    snap_num:         snapshot index
    min_stellar_mass: minimum stellar mass in Msun (after h correction)
    min_particles:    with `particle_mass` (10^10 Msun/h), cut instead at
                      the stellar mass where this snapshot's haloes are
                      resolved with that many particles
                      (visage.utils.resolution)
    max_galaxies:     random downsample if more galaxies than this are found
    hubble_h:         Hubble parameter h (from par file)
    """
    hdf5_path = Path(hdf5_path)
    group_key = f"Snap_{snap_num}"

    if VERBOSE:
        n_files = len(galaxy_files(hdf5_path))
        which = (
            hdf5_path.name
            if n_files == 1
            else f"{hdf5_path.name} + {n_files - 1} more MPI rank files"
        )
        print(f"  Galaxies: reading {which} / {group_key}...")
    with open_snapshot(hdf5_path, snap_num) as grp:
        if grp is None:
            if VERBOSE:
                print(
                    f"  Galaxies: group {group_key} not found — empty snapshot"
                )
            return GalaxySnapshot.empty(snap_num)

        def _get(field: str) -> np.ndarray:
            return np.array(grp[field])

        try:
            posx = _get("Posx")
            posy = _get("Posy")
            posz = _get("Posz")
            stellar_raw = _get("StellarMass")
            bulge_raw = _get("BulgeMass")
            cold_gas_raw = _get("ColdGas")
            mvir_raw = _get("Mvir")
            sfr_disk_arr = _get("SfrDisk")
            sfr_bulge_arr = _get("SfrBulge")
            gal_type = _get("Type").astype(np.int32)
        except KeyError as e:
            raise KeyError(
                f"Missing field {e} in {hdf5_path}:{group_key}"
            ) from e

        # Optional fields — present in current SAGE outputs but tolerate absence
        def _opt(field: str, default_dtype) -> np.ndarray:
            if field in grp:
                return np.array(grp[field])
            return np.zeros(len(posx), dtype=default_dtype)

        bh_raw = _opt("BlackHoleMass", np.float32)
        ics_raw = _opt("IntraClusterStars", np.float32)
        ffb_regime_raw = _opt("FFBRegime", np.int32).astype(np.int32)
        cgm_regime_raw = _opt("Regime", np.int32).astype(np.int32)
        cmvir_raw = _opt("CentralMvir", np.float32)
        rvir_raw = _opt("Rvir", np.float32)  # Mpc/h, no unit conversion
        h2_raw = _opt("H2gas", np.float32)
        cgm_gas_raw = _opt("CGMgas", np.float32)
        hot_gas_raw = _opt("HotGas", np.float32)
        galid_raw = _opt("GalaxyIndex", np.int64).astype(np.int64)
        cid_raw = _opt("CentralGalaxyIndex", np.int64).astype(np.int64)
        tinfall_raw = _opt("TimeOfInfall", np.int32).astype(np.int32)
        # Halo structural
        len_raw = _opt("Len", np.int32).astype(np.int32)
        vmax_raw = _opt("Vmax", np.float32)
        conc_raw = _opt("Concentration", np.float32)
        spin_raw = _opt("Spin", np.float32)
        # Galaxy structural
        disk_rad_raw = _opt("DiskRadius", np.float32)
        bulge_rad_raw = _opt("BulgeRadius", np.float32)
        mbm_raw = _opt("MergerBulgeMass", np.float32)
        mbr_raw = _opt("MergerBulgeRadius", np.float32)
        ibm_raw = _opt("InstabilityBulgeMass", np.float32)
        ibr_raw = _opt("InstabilityBulgeRadius", np.float32)
        # Gas / outflows
        h1_raw = _opt("H1gas", np.float32)
        ejected_raw = _opt("EjectedMass", np.float32)
        outflow_raw = _opt("OutflowRate", np.float32)
        massload_raw = _opt("MassLoading", np.float32)
        cooling_raw = _opt("Cooling", np.float32)
        heating_raw = _opt("Heating", np.float32)
        # SFR components (metallicity fields are dimensionless fractions)
        sfr_bulge_z_raw = _opt("SfrBulgeZ", np.float32)
        sfr_disk_z_raw = _opt("SfrDiskZ", np.float32)
        # Metals
        mcg_raw = _opt("MetalsColdGas", np.float32)
        msm_raw = _opt("MetalsStellarMass", np.float32)
        mbm_met_raw = _opt("MetalsBulgeMass", np.float32)
        mhg_raw = _opt("MetalsHotGas", np.float32)
        mem_raw = _opt("MetalsEjectedMass", np.float32)
        misc_raw = _opt("MetalsIntraClusterStars", np.float32)
        mcgm_raw = _opt("MetalsCGMgas", np.float32)

        # SFH arrays (one row per galaxy × ~64 time bins).  Only read if
        # present; the age computation tolerates absent SFH gracefully.
        sfh_disk_name, sfh_bulge_name = find_sfh_datasets(grp)
        has_sfh = sfh_disk_name is not None and sfh_bulge_name is not None
        sfh_disk_raw = np.asarray(grp[sfh_disk_name]) if has_sfh else None
        sfh_bulge_raw = np.asarray(grp[sfh_bulge_name]) if has_sfh else None

        # Anything else SAGE wrote — read as-is, so a new output field is
        # available to the UI without a code change here.
        extra_raw = {
            name: np.asarray(grp[name])
            for name in discover_extra_fields(grp, len(posx))
        }

    # All mass fields stored as 10^10 Msun/h → convert to Msun
    f = 1.0e10 / hubble_h
    stellar_mass = stellar_raw.astype(np.float32) * f
    bulge_mass = bulge_raw.astype(np.float32) * f
    cold_gas = cold_gas_raw.astype(np.float32) * f
    bh_mass = bh_raw.astype(np.float32) * f
    ics_mass = ics_raw.astype(np.float32) * f
    central_mvir = cmvir_raw.astype(np.float32) * f
    h2_mass = h2_raw.astype(np.float32) * f
    cgm_gas = cgm_gas_raw.astype(np.float32) * f
    hot_gas = hot_gas_raw.astype(np.float32) * f
    sfr = (sfr_disk_arr + sfr_bulge_arr).astype(np.float32)
    sfr_bulge = sfr_bulge_arr.astype(np.float32)
    sfr_disk = sfr_disk_arr.astype(np.float32)
    ssfr = sfr / np.where(stellar_mass > 0, stellar_mass, np.inf)
    # New mass fields (10^10 Msun/h → Msun)
    h1_gas = h1_raw.astype(np.float32) * f
    ejected_mass = ejected_raw.astype(np.float32) * f
    merger_bulge_mass = mbm_raw.astype(np.float32) * f
    instability_bulge_mass = ibm_raw.astype(np.float32) * f
    metals_cold_gas = mcg_raw.astype(np.float32) * f
    metals_stellar_mass = msm_raw.astype(np.float32) * f
    metals_bulge_mass = mbm_met_raw.astype(np.float32) * f
    metals_hot_gas = mhg_raw.astype(np.float32) * f
    metals_ejected_mass = mem_raw.astype(np.float32) * f
    metals_ics = misc_raw.astype(np.float32) * f
    metals_cgm_gas = mcgm_raw.astype(np.float32) * f
    # Fields already in correct units (no conversion)
    disk_radius = disk_rad_raw.astype(np.float32)  # Mpc/h
    bulge_radius = bulge_rad_raw.astype(np.float32)  # Mpc/h
    merger_bulge_radius = mbr_raw.astype(np.float32)  # Mpc/h
    instability_bulge_radius = ibr_raw.astype(np.float32)  # Mpc/h
    vmax = vmax_raw.astype(np.float32)  # km/s
    concentration = conc_raw.astype(np.float32)  # dimensionless
    spin = spin_raw.astype(np.float32)  # dimensionless
    outflow_rate = outflow_raw.astype(np.float32)  # Msun/yr
    mass_loading = massload_raw.astype(np.float32)  # dimensionless
    cooling = cooling_raw.astype(np.float32)  # SAGE units
    heating = heating_raw.astype(np.float32)  # SAGE units
    sfr_bulge_z = sfr_bulge_z_raw.astype(np.float32)  # dimensionless
    sfr_disk_z = sfr_disk_z_raw.astype(np.float32)  # dimensionless

    if min_particles and particle_mass:
        from visage.utils.resolution import resolved_stellar_mass

        cut = resolved_stellar_mass(
            stellar_mass,
            mvir_raw / particle_mass,
            gal_type == 0,
            min_particles,
        )
        min_stellar_mass = cut if cut is not None else 0.0
        if VERBOSE:
            print(
                f"  Galaxies: resolved (haloes >= {min_particles} particles)"
                f" from M* = {min_stellar_mass:.2e} Msun"
            )
    mask = (stellar_mass > min_stellar_mass) & (mvir_raw > 0)
    indices = np.where(mask)[0]

    # max_galaxies=None → load every galaxy (no downsample). All snapshots are
    # preloaded, so there's no display cap.
    if max_galaxies is not None and len(indices) > max_galaxies:
        rng = np.random.default_rng(42)
        indices = rng.choice(indices, max_galaxies, replace=False)

    if VERBOSE:
        print(f"  Galaxies: {len(indices):,} loaded")
    positions = np.column_stack(
        [posx[indices], posy[indices], posz[indices]]
    ).astype(np.float32)

    # Mass-weighted stellar age, Gyr (one value per galaxy).  Needs SFH arrays
    # and scale factors; otherwise filled with zeros.
    mean_age = _compute_mean_ages(
        sfh_disk=sfh_disk_raw,
        sfh_bulge=sfh_bulge_raw,
        indices=indices,
        snap_num=snap_num,
        scale_factors=scale_factors,
        hubble_h=hubble_h,
        omega_m=omega_m,
        omega_l=omega_l,
    )

    return GalaxySnapshot(
        positions=positions,
        stellar_mass=stellar_mass[indices],
        mvir=mvir_raw[indices].astype(np.float32),
        rvir=rvir_raw[indices].astype(np.float32),
        sfr=sfr[indices],
        ssfr=ssfr[indices],
        cold_gas=cold_gas[indices],
        bulge_mass=bulge_mass[indices],
        gal_type=gal_type[indices],
        bh_mass=bh_mass[indices],
        ics_mass=ics_mass[indices],
        ffb_regime=ffb_regime_raw[indices],
        cgm_regime=cgm_regime_raw[indices],
        central_mvir=central_mvir[indices],
        h2_mass=h2_mass[indices],
        cgm_gas=cgm_gas[indices],
        hot_gas=hot_gas[indices],
        galaxy_id=galid_raw[indices],
        central_id=cid_raw[indices],
        time_of_infall=tinfall_raw[indices],
        mean_age=mean_age,
        len_particles=len_raw[indices],
        vmax=vmax[indices],
        concentration=concentration[indices],
        spin=spin[indices],
        disk_radius=disk_radius[indices],
        bulge_radius=bulge_radius[indices],
        merger_bulge_mass=merger_bulge_mass[indices],
        merger_bulge_radius=merger_bulge_radius[indices],
        instability_bulge_mass=instability_bulge_mass[indices],
        instability_bulge_radius=instability_bulge_radius[indices],
        h1_gas=h1_gas[indices],
        ejected_mass=ejected_mass[indices],
        outflow_rate=outflow_rate[indices],
        mass_loading=mass_loading[indices],
        cooling=cooling[indices],
        heating=heating[indices],
        sfr_bulge=sfr_bulge[indices],
        sfr_disk=sfr_disk[indices],
        sfr_bulge_z=sfr_bulge_z[indices],
        sfr_disk_z=sfr_disk_z[indices],
        metals_cold_gas=metals_cold_gas[indices],
        metals_stellar_mass=metals_stellar_mass[indices],
        metals_bulge_mass=metals_bulge_mass[indices],
        metals_hot_gas=metals_hot_gas[indices],
        metals_ejected_mass=metals_ejected_mass[indices],
        metals_ics=metals_ics[indices],
        metals_cgm_gas=metals_cgm_gas[indices],
        extra={k: v[indices] for k, v in extra_raw.items()},
        sage_indices=indices.astype(np.int64),
        snap_num=snap_num,
        source_path=str(hdf5_path),
        has_sfh=sfh_disk_raw is not None and sfh_bulge_raw is not None,
    )

visage.io.galaxy_reader.GalaxySnapshot dataclass

Source code in visage/io/galaxy_reader.py

@dataclass
class GalaxySnapshot:
    positions: np.ndarray  # (N, 3) float32, Mpc/h
    stellar_mass: np.ndarray  # (N,)   float32, Msun
    mvir: np.ndarray  # (N,)   float32, 10^10 Msun/h (raw)
    rvir: np.ndarray  # (N,)   float32, Mpc/h (subhalo virial radius)
    sfr: np.ndarray  # (N,)   float32, Msun/yr
    ssfr: np.ndarray  # (N,)   float32, yr^-1
    cold_gas: np.ndarray  # (N,)   float32, Msun
    bulge_mass: np.ndarray  # (N,)   float32, Msun
    gal_type: np.ndarray  # (N,)   int32, 0=central 1+=satellite
    bh_mass: np.ndarray  # (N,)   float32, Msun
    ics_mass: np.ndarray  # (N,)   float32, Msun (intra-cluster stars)
    ffb_regime: np.ndarray  # (N,)   int32, FFB regime flag
    cgm_regime: np.ndarray  # (N,)   int32, 0=cold 1=hot (Regime field)
    central_mvir: np.ndarray  # (N,)   float32, Msun (host FOF Mvir)
    h2_mass: np.ndarray  # (N,)   float32, Msun
    cgm_gas: np.ndarray  # (N,)   float32, Msun
    hot_gas: np.ndarray  # (N,)   float32, Msun
    galaxy_id: np.ndarray  # (N,)   int64
    central_id: np.ndarray  # (N,)   int64
    time_of_infall: np.ndarray  # (N,)   int32
    mean_age: np.ndarray  # (N,)   float32, Gyr
    # ── Halo structural (per-galaxy from merger trees) ────────────────
    len_particles: np.ndarray  # (N,)   int32,   DM particle count
    vmax: np.ndarray  # (N,)   float32, km/s
    concentration: np.ndarray  # (N,)   float32, NFW concentration
    spin: np.ndarray  # (N,)   float32, dimensionless
    # ── Galaxy structural ─────────────────────────────────────────────
    disk_radius: np.ndarray  # (N,) float32, Mpc/h
    bulge_radius: np.ndarray  # (N,) float32, Mpc/h
    merger_bulge_mass: np.ndarray  # (N,) float32, Msun
    merger_bulge_radius: np.ndarray  # (N,) float32, Mpc/h
    instability_bulge_mass: np.ndarray  # (N,) float32, Msun
    instability_bulge_radius: np.ndarray  # (N,) float32, Mpc/h
    # ── Gas / outflows ────────────────────────────────────────────────
    h1_gas: np.ndarray  # (N,) float32, Msun  (HI gas)
    ejected_mass: np.ndarray  # (N,) float32, Msun
    outflow_rate: np.ndarray  # (N,) float32, Msun/yr
    mass_loading: np.ndarray  # (N,) float32, dimensionless
    cooling: np.ndarray  # (N,) float32, internal SAGE units
    heating: np.ndarray  # (N,) float32, internal SAGE units
    # ── SFR components ───────────────────────────────────────────────
    sfr_bulge: np.ndarray  # (N,) float32, Msun/yr (bulge SFR)
    sfr_disk: np.ndarray  # (N,) float32, Msun/yr (disk SFR)
    sfr_bulge_z: (
        np.ndarray
    )  # (N,) float32, dimensionless (bulge SFR metallicity)
    sfr_disk_z: (
        np.ndarray
    )  # (N,) float32, dimensionless (disk SFR metallicity)
    # ── Metals ───────────────────────────────────────────────────────
    metals_cold_gas: np.ndarray  # (N,) float32, Msun
    metals_stellar_mass: np.ndarray  # (N,) float32, Msun
    metals_bulge_mass: np.ndarray  # (N,) float32, Msun
    metals_hot_gas: np.ndarray  # (N,) float32, Msun
    metals_ejected_mass: np.ndarray  # (N,) float32, Msun
    metals_ics: np.ndarray  # (N,) float32, Msun
    metals_cgm_gas: np.ndarray  # (N,) float32, Msun
    # ─────────────────────────────────────────────────────────────────
    sage_indices: np.ndarray  # (N,) int64 — row indices in raw HDF5 snap group
    snap_num: int
    # ── Synthetic photometry (LightSAGE lightcones only) ────────────────
    # key -> (N,) float32 AB magnitude, e.g. "mag_rest_sdss_g",
    # "mag_obs_sdss_r". Empty for every non-SED-enabled snapshot/model.
    sed_mags: dict = field(default_factory=dict)
    # ── Auto-discovered properties ──────────────────────────────────────
    # Every other per-galaxy dataset found in the snapshot group, keyed by
    # its HDF5 name.  This is what keeps the UI honest when SAGE grows a
    # new output field: it shows up as a filter and a Colour-by mode with
    # no code change here.  Values are raw (no unit conversion is possible
    # without knowing what the field means).
    extra: dict = field(default_factory=dict)
    # ── Rendering ───────────────────────────────────────────────────────
    # How many galaxy glows overlap at each galaxy (utils.sizing.
    # glow_crowding), filled in by the loader's background thread.
    crowding: np.ndarray | None = None
    # ── Star-formation history, summarised ─────────────────────────────
    # Stellar mass formed (raw, 10^10 Msun/h) in each lookback-time bin of
    # SFH_AGE_EDGES, for the disk ([:, 0]) and the bulge ([:, 1]): shape
    # (N, 2, len(SFH_AGE_EDGES) - 1).  None without SFHMassDisk/Bulge
    # (SAGE's SaveFullSFH off), and None by default with them too: HD reads
    # the few rows it needs on demand (sfh_bins_for) rather than every
    # snapshot load binning every galaxy.
    sfh_bins: np.ndarray | None = None
    # Where the snapshot came from, and whether it has SFH arrays.
    source_path: str = ""
    has_sfh: bool = False
    # ── Lightcones: each galaxy is from its own snapshot ────────────────
    # (N,) source snapshot and cosmological redshift per galaxy; None in a
    # box, where every galaxy is at `snap_num`.
    snap_nums: np.ndarray | None = None
    redshifts: np.ndarray | None = None
    # The cone's PhotometryHeader (how its synthetic photometry was made;
    # empty without photometry).
    photometry: dict = field(default_factory=dict)

    def sed_dust(self) -> tuple[str | None, float]:
        """The dust in this catalogue's synthetic photometry: ("charlot_fall",
        0) for each galaxy's own (TauV_ISM, CosInclination), ("calzetti",
        tau_V) for one screen for all, ("none", 0), or (None, 0) without
        photometry."""
        desc = str(self.photometry.get("Dust", "")) if self.photometry else ""
        if not self.photometry:
            return None, 0.0
        low = desc.lower()
        if low.startswith("charlot"):
            return "charlot_fall", 0.0
        if low.startswith("calzetti"):
            import re

            m = re.search(r"tau_v\s*=\s*([0-9.eE+-]+)", desc, re.I)
            return "calzetti", float(m.group(1)) if m else 0.0
        if low == "none":
            return "none", 0.0
        return None, 0.0

    @property
    def is_cone(self) -> bool:
        return self.snap_nums is not None

    def snap_of(self, i: int) -> int:
        """The snapshot galaxy `i` is from."""
        if self.snap_nums is not None:
            return int(self.snap_nums[i])
        return int(self.snap_num)

    @property
    def count(self) -> int:
        return len(self.positions)

    @classmethod
    def empty(cls, snap_num: int) -> GalaxySnapshot:
        z = np.empty(0, dtype=np.float32)
        zi = np.empty(0, dtype=np.int32)
        zi64 = np.empty(0, dtype=np.int64)
        return cls(
            positions=np.empty((0, 3), dtype=np.float32),
            stellar_mass=z,
            mvir=z,
            rvir=z,
            sfr=z,
            ssfr=z,
            cold_gas=z,
            bulge_mass=z,
            gal_type=zi,
            bh_mass=z,
            ics_mass=z,
            central_mvir=z,
            ffb_regime=zi,
            cgm_regime=zi,
            h2_mass=z,
            cgm_gas=z,
            hot_gas=z,
            galaxy_id=zi64,
            central_id=zi64,
            time_of_infall=zi,
            mean_age=z,
            len_particles=zi,
            vmax=z,
            concentration=z,
            spin=z,
            disk_radius=z,
            bulge_radius=z,
            merger_bulge_mass=z,
            merger_bulge_radius=z,
            instability_bulge_mass=z,
            instability_bulge_radius=z,
            h1_gas=z,
            ejected_mass=z,
            outflow_rate=z,
            mass_loading=z,
            cooling=z,
            heating=z,
            sfr_bulge=z,
            sfr_disk=z,
            sfr_bulge_z=z,
            sfr_disk_z=z,
            metals_cold_gas=z,
            metals_stellar_mass=z,
            metals_bulge_mass=z,
            metals_hot_gas=z,
            metals_ejected_mass=z,
            metals_ics=z,
            metals_cgm_gas=z,
            sage_indices=np.empty(0, dtype=np.int64),
            snap_num=snap_num,
        )

sed_dust()

The dust in this catalogue's synthetic photometry: ("charlot_fall", 0) for each galaxy's own (TauV_ISM, CosInclination), ("calzetti", tau_V) for one screen for all, ("none", 0), or (None, 0) without photometry.

Source code in visage/io/galaxy_reader.py

def sed_dust(self) -> tuple[str | None, float]:
    """The dust in this catalogue's synthetic photometry: ("charlot_fall",
    0) for each galaxy's own (TauV_ISM, CosInclination), ("calzetti",
    tau_V) for one screen for all, ("none", 0), or (None, 0) without
    photometry."""
    desc = str(self.photometry.get("Dust", "")) if self.photometry else ""
    if not self.photometry:
        return None, 0.0
    low = desc.lower()
    if low.startswith("charlot"):
        return "charlot_fall", 0.0
    if low.startswith("calzetti"):
        import re

        m = re.search(r"tau_v\s*=\s*([0-9.eE+-]+)", desc, re.I)
        return "calzetti", float(m.group(1)) if m else 0.0
    if low == "none":
        return "none", 0.0
    return None, 0.0

snap_of(i)

The snapshot galaxy i is from.

Source code in visage/io/galaxy_reader.py

def snap_of(self, i: int) -> int:
    """The snapshot galaxy `i` is from."""
    if self.snap_nums is not None:
        return int(self.snap_nums[i])
    return int(self.snap_num)

Scene layer

visage.scene.scene

visage.scene.scene.Scene

Owner of the PyVista plotter and a dict of Models.

One model is the primary (its layers are exposed via scene.halo_layer / scene.galaxy_layer so all existing UI code keeps working). Additional models can be loaded and made visible as overlays if they share the same box size and snapshot count, or placed side-by-side as independent adjacent boxes with their own snapshot, filters and render settings.

Playback animation is driven externally by the Trame async event loop (see toolbar.py) so that all VTK calls stay on the main thread.

Source code in visage/scene/scene.py

class Scene:
    """Owner of the PyVista plotter and a dict of Models.

    One model is the *primary* (its layers are exposed via `scene.halo_layer` /
    `scene.galaxy_layer` so all existing UI code keeps working). Additional
    models can be loaded and made visible as overlays if they share the same
    box size and snapshot count, or placed side-by-side as independent
    adjacent boxes with their own snapshot, filters and render settings.

    Playback animation is driven externally by the Trame async event loop
    (see toolbar.py) so that all VTK calls stay on the main thread.
    """

    def __init__(
        self,
        primary_par_path: str | Path,
        off_screen: bool = False,
        initial_snap: int | None = None,
        n_jobs: int = -1,
        min_halo_mass: float = 1.0e10,
        min_stellar_mass: float = 1.0e8,
        max_halos: int | None = None,
        max_galaxies: int | None = None,
        cache_bytes: int | None = None,
        lightcone_path: str | Path | None = None,
        min_particles: int | None = None,
    ) -> None:
        self._plotter = pv.Plotter(
            off_screen=off_screen, window_size=[1600, 900]
        )
        # A hint of blue at the bottom rather than flat black: the box
        # then sits in something instead of floating in a void.  Subtle
        # enough not to tint the data.
        self._plotter.set_background("#080b14", top="#000000")
        # FXAA: a post-process pass, so it smooths the splat edges and the
        # box labels without the 4x cost of supersampling — and without
        # touching the render resolution the remote view ships.
        try:
            self._plotter.enable_anti_aliasing("fxaa")
        except Exception:
            # Some VTK builds / headless GL stacks have no FXAA pass.
            pass
        self._plotter.renderer.SetNearClippingPlaneTolerance(0.00001)

        # Loader kwargs reused when adding additional models later
        self._loader_kwargs = dict(
            n_jobs=n_jobs,
            min_halo_mass=min_halo_mass,
            min_stellar_mass=min_stellar_mass,
            max_halos=max_halos,
            max_galaxies=max_galaxies,
            cache_bytes=cache_bytes,
            min_particles=min_particles,
        )

        self._models: dict[str, Model] = {}
        self._primary_name: str = ""
        self._is_lightcone: bool = lightcone_path is not None

        # Build the primary model — a lightcone (flat cone) or a SAGE box.
        if self._is_lightcone:
            from visage.scene.lightcone_model import LightconeModel

            primary = LightconeModel(
                lightcone_path, self._plotter, self._loader_kwargs
            )
        else:
            primary = Model(
                primary_par_path, self._plotter, self._loader_kwargs
            )
        self._models[primary.name] = primary
        self._primary_name = primary.name

        # Adjacent-box state
        self._adjacent_order: list[str] = []  # names in placement order
        self._active_box_name: str = primary.name

        # Label actors keyed by model name (only shown when multiple boxes are loaded)
        self._label_actors: dict[str, list] = {}
        # When False, the per-box name/redshift labels are suppressed entirely
        # (e.g. Story Mode prefers a clean view). Restored on exit.
        self._labels_enabled: bool = True

        self._camera = CameraController(self._plotter, primary.box_size)
        if self._is_lightcone:
            # Frame the actual cone bounds, not a cube.
            self._camera.set_lightcone_bounds(getattr(primary, "bounds", None))

        if initial_snap is not None:
            self._current_snap: int = initial_snap
        elif self._is_lightcone:
            # Full cone by default: the near-side cut is "off" at snap_max.
            self._current_snap = getattr(
                primary, "snap_max", primary.snap_count - 1
            )
        else:
            self._current_snap = primary.snap_count - 1

        self._on_snap_change: list[Callable[[int], None]] = []
        self._on_model_change: list[Callable[[], None]] = []
        self._focus_region: dict | None = None

        self.set_snapshot(self._current_snap)
        self._camera.reset()

        # Start background preloading of all snapshots immediately so
        # navigation is stall-free by the time the browser opens.
        # Suppress per-snapshot log lines so they don't bury the server URL.
        from visage.io import halo_reader, galaxy_reader

        halo_reader.VERBOSE = False
        galaxy_reader.VERBOSE = False
        # A lightcone is a single static cloud — nothing to preload.
        if not self._is_lightcone:
            self.primary.loader.preload_all()

    # ------------------------------------------------------------------
    # Primary model & layer access
    # ------------------------------------------------------------------

    @property
    def is_lightcone(self) -> bool:
        return self._is_lightcone

    @property
    def primary(self) -> Model:
        return self._models[self._primary_name]

    @property
    def primary_name(self) -> str:
        return self._primary_name

    @property
    def active_model(self) -> Model:
        """The model whose settings the UI panel currently controls."""
        return self._models.get(self._active_box_name, self.primary)

    # halo_layer / galaxy_layer route to the ACTIVE model so all existing
    # filter and rendering handlers automatically target the right box.
    @property
    def halo_layer(self) -> HaloLayer:
        return self.active_model.halo_layer

    @property
    def galaxy_layer(self) -> GalaxyLayer:
        return self.active_model.galaxy_layer

    @property
    def camera(self) -> CameraController:
        return self._camera

    @property
    def plotter(self) -> pv.Plotter:
        return self._plotter

    @property
    def current_snap(self) -> int:
        return self.active_model.current_snap

    @property
    def snap_label(self) -> str:
        m = self.active_model
        return m.snap_table.label(m.current_snap)

    # Backward-compat shims used by various callers
    @property
    def _cfg(self):
        return self.primary.cfg

    @property
    def _snap_table(self):
        return self.primary.snap_table

    @property
    def _loader(self):
        return self.primary.loader

    @property
    def snap_count(self) -> int:
        return self.primary.snap_count

    # ------------------------------------------------------------------
    # Active-box management
    # ------------------------------------------------------------------

    @property
    def active_box_name(self) -> str:
        return self._active_box_name

    def set_active_box(self, name: str) -> None:
        """Switch which box the UI controls.
        Profile save/restore is handled by the app layer."""
        if name not in self._models:
            return
        self._active_box_name = name
        self._update_labels()
        self._plotter.render()

    # ------------------------------------------------------------------
    # Adjacent-box management
    # ------------------------------------------------------------------

    def is_adjacent(self, name: str) -> bool:
        return name in self._adjacent_order

    def toggle_adjacent(self, par_path: str | Path) -> tuple[bool, str | None]:
        """Add or remove a model as an adjacent side-by-side box.

        Returns (is_now_adjacent, error_message_or_None).
        """
        # Load the model if not already known
        if not self.has_model_by_path(par_path):
            model = Model(par_path, self._plotter, self._loader_kwargs)
            self._models[model.name] = model
        else:
            model = self._model_by_path(par_path)

        name = model.name
        if name == self._primary_name:
            return False, "Cannot place the primary model as adjacent."

        if name in self._adjacent_order:
            # Remove it
            self._adjacent_order.remove(name)
            if self._active_box_name == name:
                self._active_box_name = self._primary_name
            model.offset = np.zeros(3)
            model.visible = False
            self._remove_label(name)
            self._recompute_offsets()
            self._update_labels()
            self._plotter.render()
            for cb in self._on_model_change:
                cb()
            return False, None
        else:
            # Add it
            self._adjacent_order.append(name)
            self._recompute_offsets()
            # Always start adjacent boxes with default color/colormap so they
            # don't inherit whatever the primary box had been set to.
            model.halo_layer.color_mode = "mvir"
            model.halo_layer.colormap = "viridis"
            model.galaxy_layer.color_mode = "structure"
            model.galaxy_layer.colormap = "plasma"
            model.set_snapshot(model.snap_count - 1)
            model.visible = True
            self._update_labels()
            self._plotter.render()
            for cb in self._on_model_change:
                cb()
            return True, None

    def _recompute_offsets(self) -> None:
        """Lay adjacent boxes out along +X with a small gap."""
        gap = self.primary.box_size * _BOX_GAP_FRACTION
        x = self.primary.box_size + gap
        for name in self._adjacent_order:
            m = self._models.get(name)
            if m is None:
                continue
            m.offset = np.array([x, 0.0, 0.0])
            x += m.box_size + gap

    def has_model_by_path(self, par_path: str | Path) -> bool:
        p = Path(par_path)
        return any(m.path == p for m in self._models.values())

    def _model_by_path(self, par_path: str | Path) -> Model:
        p = Path(par_path)
        for m in self._models.values():
            if m.path == p:
                return m
        raise KeyError(par_path)

    # ------------------------------------------------------------------
    # 3-D text labels (model name + redshift, below each box)
    # Shown only when multiple boxes are loaded.
    # ------------------------------------------------------------------

    def _box_center_xz(self, name: str) -> tuple[float, float, float]:
        """Return (cx, 0, cz) — the XZ centre of the named box."""
        m = self._models[name]
        off = m.offset
        bs = m.box_size
        return float(off[0] + bs / 2), 0.0, float(off[2] + bs / 2)

    def _label_position(self, name: str) -> tuple[float, float, float]:
        cx, _, cz = self._box_center_xz(name)
        m = self._models[name]
        return cx, -m.box_size * 0.22, cz

    def _label_text(self, name: str) -> str:
        m = self._models[name]
        snap = max(0, m.current_snap)
        z = m.snap_table.snap_to_z(snap)
        return f"{m.name}  z={z:.2f}"

    def _remove_label(self, name: str) -> None:
        for actor in self._label_actors.pop(name, []):
            self._plotter.renderer.RemoveActor(actor)

    def _add_label(self, name: str) -> None:
        from vtkmodules.vtkRenderingCore import vtkBillboardTextActor3D

        cx, y, cz = self._label_position(name)
        text = self._label_text(name)

        actor = vtkBillboardTextActor3D()
        actor.SetInput(text)
        actor.SetPosition(cx, y, cz)

        tp = actor.GetTextProperty()
        if name == self._active_box_name:
            tp.SetColor(0.2, 1.0, 0.4)  # green for active
        else:
            tp.SetColor(1.0, 1.0, 1.0)  # white for idle
        tp.SetFontSize(14)
        tp.SetJustificationToCentered()
        tp.SetVerticalJustificationToCentered()
        tp.SetBackgroundOpacity(0.0)
        tp.SetBold(False)
        tp.SetShadow(False)

        self._plotter.renderer.AddActor(actor)
        self._label_actors[name] = [actor]

    def _update_labels(self) -> None:
        """Rebuild labels for all boxes — only when more than one box is loaded."""
        if not self._labels_enabled:
            for name in list(self._label_actors):
                self._remove_label(name)
            self._plotter.render()
            return
        if len(self._adjacent_order) == 0:
            # Back to a single box: remove any lingering labels
            for name in list(self._label_actors):
                self._remove_label(name)
            self._plotter.render()
            return
        label_names = [self._primary_name] + list(self._adjacent_order)
        for name in list(self._label_actors):
            if name not in label_names:
                self._remove_label(name)
        for name in label_names:
            self._remove_label(name)
            self._add_label(name)
        self._plotter.render()

    def set_labels_enabled(self, enabled: bool) -> None:
        """Enable/disable the per-box name+redshift labels and rebuild/clear."""
        enabled = bool(enabled)
        if enabled == self._labels_enabled:
            return
        self._labels_enabled = enabled
        self._update_labels()

    def refresh_label(self, name: str) -> None:
        """Refresh just one label (e.g. after a snapshot change)."""
        if not self._labels_enabled:
            return
        if len(self._adjacent_order) == 0:
            return
        if name in self._label_actors or name in (
            [self._primary_name] + self._adjacent_order
        ):
            self._remove_label(name)
            self._add_label(name)
            self._plotter.render()

    # ------------------------------------------------------------------
    # Overlay model management (unchanged from before)
    # ------------------------------------------------------------------

    def list_models(self) -> list[Model]:
        return list(self._models.values())

    def has_model(self, name: str) -> bool:
        return name in self._models

    def add_model(self, par_path: str | Path) -> Model:
        """Load a new model alongside the primary. Starts hidden by default."""
        model = Model(par_path, self._plotter, self._loader_kwargs)
        if model.name in self._models:
            return self._models[model.name]
        model.set_snapshot(self._current_snap)
        model.visible = False
        self._models[model.name] = model
        for cb in self._on_model_change:
            cb()
        return model

    def remove_model(self, name: str) -> None:
        if name == self._primary_name or name not in self._models:
            return
        model = self._models.pop(name)
        model.visible = False
        model.halo_layer._clear_actors()
        model.galaxy_layer._clear_actors()
        model.shutdown()
        if name in self._adjacent_order:
            self._adjacent_order.remove(name)
        self._remove_label(name)
        for cb in self._on_model_change:
            cb()

    def switch_primary(self, name: str) -> None:
        """Switch the active primary model.

        Adjacent boxes stay in place; their offsets are recomputed relative
        to the new primary.  If the active box was the old primary, the
        active box switches to the new primary.
        """
        if name == self._primary_name or name not in self._models:
            return
        old_primary_box = self.primary.box_size
        self._models[self._primary_name].visible = False

        # If old primary was active, transfer focus to new primary
        if self._active_box_name == self._primary_name:
            self._active_box_name = name

        self._primary_name = name

        # Hide overlays that are incompatible with the new primary
        for other_name, m in self._models.items():
            if other_name == self._primary_name:
                continue
            if other_name in self._adjacent_order:
                continue  # adjacent boxes are not overlay-checked
            if m.visible and not self.is_compatible_for_overlay(other_name):
                m.visible = False

        self.primary.visible = True
        snap = self.primary.snap_count - 1
        self._current_snap = snap
        self.primary.set_snapshot(snap)

        new_box = self.primary.box_size
        self._camera._box_size = new_box
        if abs(new_box - old_primary_box) > 1e-3:
            self._camera.reset()

        # Recompute adjacent offsets relative to new primary
        self._recompute_offsets()
        self._update_labels()

        if self._focus_region is not None:
            halos, galaxies = self.primary.loader.get(self._current_snap)
            self._apply_focus_masks_for_layer(
                self.primary.halo_layer,
                self.primary.galaxy_layer,
                halos.positions,
                galaxies.positions,
            )
        for cb in self._on_snap_change:
            cb(self._current_snap)
        for cb in self._on_model_change:
            cb()

    def set_overlay_visible(self, name: str, vis: bool) -> str | None:
        if name == self._primary_name or name not in self._models:
            return None
        if name in self._adjacent_order:
            return "Use the Side-by-Side toggle to manage adjacent boxes."
        m = self._models[name]
        if vis:
            if not self.is_compatible_for_overlay(name):
                primary = self.primary
                cand = self._models[name]
                if abs(primary.box_size - cand.box_size) > 1e-3:
                    detail = (
                        f"box size {cand.box_size:.1f} ≠ "
                        f"{primary.box_size:.1f} Mpc/h"
                    )
                elif primary.snap_count != cand.snap_count:
                    detail = (
                        f"snapshot count {cand.snap_count} ≠ "
                        f"{primary.snap_count}"
                    )
                else:
                    detail = "incompatible simulation parameters"
                return (
                    f"Can't overlay '{name}' on '{primary.name}': {detail}. "
                    f"Use Side-by-Side instead."
                )
            m.set_snapshot(self._current_snap)
        m.visible = vis
        for cb in self._on_model_change:
            cb()
        return None

    def is_compatible_for_overlay(self, name: str) -> bool:
        if name == self._primary_name or name not in self._models:
            return False
        primary = self.primary
        candidate = self._models[name]
        return (
            abs(primary.box_size - candidate.box_size) < 1e-3
            and primary.snap_count == candidate.snap_count
        )

    # ------------------------------------------------------------------
    # Snapshot control
    # ------------------------------------------------------------------

    def set_snapshot(self, snap_num: int) -> None:
        """Update snapshot.

        When the active box is an adjacent model, only that model's snapshot
        changes.  When the primary is active, the original behaviour applies
        (primary + compatible overlays all update together).
        """
        if (
            self._active_box_name != self._primary_name
            and self._active_box_name in self._adjacent_order
        ):
            active = self.active_model
            snap_num = max(0, min(int(snap_num), active.snap_count - 1))
            active.set_snapshot(snap_num)
            halos, galaxies = active.loader.get(snap_num)
            off = active.offset.astype(np.float32)
            self._camera.update_halo_index(
                halos.positions + off, tree=active.loader.get_tree(snap_num)
            )
            self._camera.update_galaxy_positions(galaxies.positions + off)
            self.refresh_label(self._active_box_name)
            for cb in self._on_snap_change:
                cb(snap_num)
            return

        # Primary is active: original behaviour
        snap_num = max(0, min(int(snap_num), self.primary.snap_count - 1))
        self._current_snap = snap_num
        self.primary.set_snapshot(snap_num)

        # Update compatible overlays (not adjacent boxes — they're independent)
        for name, m in self._models.items():
            if name == self._primary_name:
                continue
            if name in self._adjacent_order:
                continue
            if m.visible:
                m.set_snapshot(min(snap_num, m.snap_count - 1))

        halos, galaxies = self.primary.loader.get(snap_num)
        self._camera.update_halo_index(
            halos.positions, tree=self.primary.loader.get_tree(snap_num)
        )
        self._camera.update_galaxy_positions(galaxies.positions)

        if self._focus_region is not None:
            self._apply_focus_masks_for_layer(
                self.primary.halo_layer,
                self.primary.galaxy_layer,
                halos.positions,
                galaxies.positions,
            )

        self.refresh_label(self._primary_name)
        for cb in self._on_snap_change:
            cb(snap_num)

    # ------------------------------------------------------------------
    # Click-to-activate: determine which box a 3-D click point belongs to
    # ------------------------------------------------------------------

    def box_name_at(self, world_x: float) -> str:
        """Return the model name whose X range contains *world_x*.

        Falls back to the primary if no adjacent box matches.
        """
        for name in reversed(self._adjacent_order):
            m = self._models[name]
            x0 = float(m.offset[0])
            x1 = x0 + m.box_size
            if x0 <= world_x <= x1:
                return name
        return self._primary_name

    # ------------------------------------------------------------------
    # Focus / spatial masking (applies to active model only)
    # ------------------------------------------------------------------

    def set_focus_box(
        self,
        xmin: float,
        xmax: float,
        ymin: float,
        ymax: float,
        zmin: float,
        zmax: float,
    ) -> None:
        self._focus_region = dict(
            type="box",
            xmin=xmin,
            xmax=xmax,
            ymin=ymin,
            ymax=ymax,
            zmin=zmin,
            zmax=zmax,
        )
        active = self.active_model
        halos, galaxies = active.loader.get(active.current_snap)
        self._apply_focus_masks_for_layer(
            active.halo_layer,
            active.galaxy_layer,
            halos.positions,
            galaxies.positions,
        )
        self._refit_clipping()

    def set_focus_sphere(
        self,
        center: tuple[float, float, float],
        radius: float,
    ) -> None:
        self._focus_region = dict(type="sphere", center=center, radius=radius)
        active = self.active_model
        halos, galaxies = active.loader.get(active.current_snap)
        off = active.offset.astype(np.float32)
        self._apply_focus_masks_for_layer(
            active.halo_layer,
            active.galaxy_layer,
            halos.positions + off,
            galaxies.positions + off,
        )
        self._refit_clipping()

    def set_focus_indices(
        self,
        halo_indices=None,
        galaxy_indices=None,
    ) -> None:
        """Isolate specific halos/galaxies — everything else is hidden.

        Each argument is a sequence of indices into the current snapshot
        (None or empty hides that population entirely).
        """
        self._focus_region = dict(
            type="indices",
            halo_indices=np.asarray(
                [] if halo_indices is None else halo_indices, dtype=np.int64
            ),
            galaxy_indices=np.asarray(
                [] if galaxy_indices is None else galaxy_indices,
                dtype=np.int64,
            ),
        )
        active = self.active_model
        halos, galaxies = active.loader.get(active.current_snap)
        self._apply_focus_masks_for_layer(
            active.halo_layer,
            active.galaxy_layer,
            halos.positions,
            galaxies.positions,
        )
        self._refit_clipping()

    def clear_focus(self) -> None:
        self._focus_region = None
        for m in self._models.values():
            m.halo_layer.set_mask(None)
            m.galaxy_layer.set_mask(None)
        self._refit_clipping()

    def _refit_clipping(self) -> None:
        """Fit the camera's depth range to what is now visible.  A fly-to
        sets the camera while the old focus region is still all that is
        visible; if that lies behind the new view, VTK collapses the range
        to ~1e-18 and nothing draws until it is refitted."""
        try:
            self.plotter.renderer.ResetCameraClippingRange()
        except Exception:
            pass

    def _apply_focus_masks_for_layer(
        self,
        halo_layer: HaloLayer,
        gal_layer: GalaxyLayer,
        halo_pos: np.ndarray,
        gal_pos: np.ndarray,
    ) -> None:
        r = self._focus_region
        if r is None:
            return
        if r["type"] == "box":

            def _box_mask(pos):
                if len(pos) == 0:
                    return np.array([], dtype=bool)
                return (
                    (pos[:, 0] >= r["xmin"])
                    & (pos[:, 0] <= r["xmax"])
                    & (pos[:, 1] >= r["ymin"])
                    & (pos[:, 1] <= r["ymax"])
                    & (pos[:, 2] >= r["zmin"])
                    & (pos[:, 2] <= r["zmax"])
                )

            halo_layer.set_mask(_box_mask(halo_pos))
            gal_layer.set_mask(_box_mask(gal_pos))
        elif r["type"] == "indices":

            def _index_mask(pos, idxs):
                mask = np.zeros(len(pos), dtype=bool)
                if len(idxs):
                    keep = idxs[(idxs >= 0) & (idxs < len(pos))]
                    mask[keep] = True
                return mask

            halo_layer.set_mask(_index_mask(halo_pos, r["halo_indices"]))
            gal_layer.set_mask(_index_mask(gal_pos, r["galaxy_indices"]))
        elif r["type"] == "sphere":
            cx, cy, cz = r["center"]
            rad = r["radius"]

            def _sphere_mask(pos):
                if len(pos) == 0:
                    return np.array([], dtype=bool)
                return (
                    np.linalg.norm(pos - np.array([cx, cy, cz]), axis=1) <= rad
                )

            halo_layer.set_mask(_sphere_mask(halo_pos))
            gal_layer.set_mask(_sphere_mask(gal_pos))

    def _apply_focus_masks(self, halo_pos, gal_pos) -> None:
        self._apply_focus_masks_for_layer(
            self.primary.halo_layer,
            self.primary.galaxy_layer,
            halo_pos,
            gal_pos,
        )

    def next_snap_num(self) -> int:
        return (
            self.active_model.current_snap + 1
        ) % self.active_model.snap_count

    def prev_snap_num(self) -> int:
        return (
            self.active_model.current_snap - 1
        ) % self.active_model.snap_count

    # ------------------------------------------------------------------
    # Callbacks
    # ------------------------------------------------------------------

    def register_snap_change_callback(self, cb: Callable[[int], None]) -> None:
        self._on_snap_change.append(cb)

    def register_model_change_callback(self, cb: Callable[[], None]) -> None:
        self._on_model_change.append(cb)

    # ------------------------------------------------------------------
    # Cleanup
    # ------------------------------------------------------------------

    def close(self) -> None:
        for m in self._models.values():
            m.shutdown()
        self._plotter.close()

active_model property

The model whose settings the UI panel currently controls.

add_model(par_path)

Load a new model alongside the primary. Starts hidden by default.

Source code in visage/scene/scene.py

def add_model(self, par_path: str | Path) -> Model:
    """Load a new model alongside the primary. Starts hidden by default."""
    model = Model(par_path, self._plotter, self._loader_kwargs)
    if model.name in self._models:
        return self._models[model.name]
    model.set_snapshot(self._current_snap)
    model.visible = False
    self._models[model.name] = model
    for cb in self._on_model_change:
        cb()
    return model

box_name_at(world_x)

Return the model name whose X range contains world_x.

Falls back to the primary if no adjacent box matches.

Source code in visage/scene/scene.py

def box_name_at(self, world_x: float) -> str:
    """Return the model name whose X range contains *world_x*.

    Falls back to the primary if no adjacent box matches.
    """
    for name in reversed(self._adjacent_order):
        m = self._models[name]
        x0 = float(m.offset[0])
        x1 = x0 + m.box_size
        if x0 <= world_x <= x1:
            return name
    return self._primary_name

refresh_label(name)

Refresh just one label (e.g. after a snapshot change).

Source code in visage/scene/scene.py

def refresh_label(self, name: str) -> None:
    """Refresh just one label (e.g. after a snapshot change)."""
    if not self._labels_enabled:
        return
    if len(self._adjacent_order) == 0:
        return
    if name in self._label_actors or name in (
        [self._primary_name] + self._adjacent_order
    ):
        self._remove_label(name)
        self._add_label(name)
        self._plotter.render()

set_active_box(name)

Switch which box the UI controls. Profile save/restore is handled by the app layer.

Source code in visage/scene/scene.py

def set_active_box(self, name: str) -> None:
    """Switch which box the UI controls.
    Profile save/restore is handled by the app layer."""
    if name not in self._models:
        return
    self._active_box_name = name
    self._update_labels()
    self._plotter.render()

set_focus_indices(halo_indices=None, galaxy_indices=None)

Isolate specific halos/galaxies — everything else is hidden.

Each argument is a sequence of indices into the current snapshot (None or empty hides that population entirely).

Source code in visage/scene/scene.py

def set_focus_indices(
    self,
    halo_indices=None,
    galaxy_indices=None,
) -> None:
    """Isolate specific halos/galaxies — everything else is hidden.

    Each argument is a sequence of indices into the current snapshot
    (None or empty hides that population entirely).
    """
    self._focus_region = dict(
        type="indices",
        halo_indices=np.asarray(
            [] if halo_indices is None else halo_indices, dtype=np.int64
        ),
        galaxy_indices=np.asarray(
            [] if galaxy_indices is None else galaxy_indices,
            dtype=np.int64,
        ),
    )
    active = self.active_model
    halos, galaxies = active.loader.get(active.current_snap)
    self._apply_focus_masks_for_layer(
        active.halo_layer,
        active.galaxy_layer,
        halos.positions,
        galaxies.positions,
    )
    self._refit_clipping()

set_labels_enabled(enabled)

Enable/disable the per-box name+redshift labels and rebuild/clear.

Source code in visage/scene/scene.py

def set_labels_enabled(self, enabled: bool) -> None:
    """Enable/disable the per-box name+redshift labels and rebuild/clear."""
    enabled = bool(enabled)
    if enabled == self._labels_enabled:
        return
    self._labels_enabled = enabled
    self._update_labels()

set_snapshot(snap_num)

Update snapshot.

When the active box is an adjacent model, only that model's snapshot changes. When the primary is active, the original behaviour applies (primary + compatible overlays all update together).

Source code in visage/scene/scene.py

def set_snapshot(self, snap_num: int) -> None:
    """Update snapshot.

    When the active box is an adjacent model, only that model's snapshot
    changes.  When the primary is active, the original behaviour applies
    (primary + compatible overlays all update together).
    """
    if (
        self._active_box_name != self._primary_name
        and self._active_box_name in self._adjacent_order
    ):
        active = self.active_model
        snap_num = max(0, min(int(snap_num), active.snap_count - 1))
        active.set_snapshot(snap_num)
        halos, galaxies = active.loader.get(snap_num)
        off = active.offset.astype(np.float32)
        self._camera.update_halo_index(
            halos.positions + off, tree=active.loader.get_tree(snap_num)
        )
        self._camera.update_galaxy_positions(galaxies.positions + off)
        self.refresh_label(self._active_box_name)
        for cb in self._on_snap_change:
            cb(snap_num)
        return

    # Primary is active: original behaviour
    snap_num = max(0, min(int(snap_num), self.primary.snap_count - 1))
    self._current_snap = snap_num
    self.primary.set_snapshot(snap_num)

    # Update compatible overlays (not adjacent boxes — they're independent)
    for name, m in self._models.items():
        if name == self._primary_name:
            continue
        if name in self._adjacent_order:
            continue
        if m.visible:
            m.set_snapshot(min(snap_num, m.snap_count - 1))

    halos, galaxies = self.primary.loader.get(snap_num)
    self._camera.update_halo_index(
        halos.positions, tree=self.primary.loader.get_tree(snap_num)
    )
    self._camera.update_galaxy_positions(galaxies.positions)

    if self._focus_region is not None:
        self._apply_focus_masks_for_layer(
            self.primary.halo_layer,
            self.primary.galaxy_layer,
            halos.positions,
            galaxies.positions,
        )

    self.refresh_label(self._primary_name)
    for cb in self._on_snap_change:
        cb(snap_num)

switch_primary(name)

Switch the active primary model.

Adjacent boxes stay in place; their offsets are recomputed relative to the new primary. If the active box was the old primary, the active box switches to the new primary.

Source code in visage/scene/scene.py

def switch_primary(self, name: str) -> None:
    """Switch the active primary model.

    Adjacent boxes stay in place; their offsets are recomputed relative
    to the new primary.  If the active box was the old primary, the
    active box switches to the new primary.
    """
    if name == self._primary_name or name not in self._models:
        return
    old_primary_box = self.primary.box_size
    self._models[self._primary_name].visible = False

    # If old primary was active, transfer focus to new primary
    if self._active_box_name == self._primary_name:
        self._active_box_name = name

    self._primary_name = name

    # Hide overlays that are incompatible with the new primary
    for other_name, m in self._models.items():
        if other_name == self._primary_name:
            continue
        if other_name in self._adjacent_order:
            continue  # adjacent boxes are not overlay-checked
        if m.visible and not self.is_compatible_for_overlay(other_name):
            m.visible = False

    self.primary.visible = True
    snap = self.primary.snap_count - 1
    self._current_snap = snap
    self.primary.set_snapshot(snap)

    new_box = self.primary.box_size
    self._camera._box_size = new_box
    if abs(new_box - old_primary_box) > 1e-3:
        self._camera.reset()

    # Recompute adjacent offsets relative to new primary
    self._recompute_offsets()
    self._update_labels()

    if self._focus_region is not None:
        halos, galaxies = self.primary.loader.get(self._current_snap)
        self._apply_focus_masks_for_layer(
            self.primary.halo_layer,
            self.primary.galaxy_layer,
            halos.positions,
            galaxies.positions,
        )
    for cb in self._on_snap_change:
        cb(self._current_snap)
    for cb in self._on_model_change:
        cb()

toggle_adjacent(par_path)

Add or remove a model as an adjacent side-by-side box.

Returns (is_now_adjacent, error_message_or_None).

Source code in visage/scene/scene.py

def toggle_adjacent(self, par_path: str | Path) -> tuple[bool, str | None]:
    """Add or remove a model as an adjacent side-by-side box.

    Returns (is_now_adjacent, error_message_or_None).
    """
    # Load the model if not already known
    if not self.has_model_by_path(par_path):
        model = Model(par_path, self._plotter, self._loader_kwargs)
        self._models[model.name] = model
    else:
        model = self._model_by_path(par_path)

    name = model.name
    if name == self._primary_name:
        return False, "Cannot place the primary model as adjacent."

    if name in self._adjacent_order:
        # Remove it
        self._adjacent_order.remove(name)
        if self._active_box_name == name:
            self._active_box_name = self._primary_name
        model.offset = np.zeros(3)
        model.visible = False
        self._remove_label(name)
        self._recompute_offsets()
        self._update_labels()
        self._plotter.render()
        for cb in self._on_model_change:
            cb()
        return False, None
    else:
        # Add it
        self._adjacent_order.append(name)
        self._recompute_offsets()
        # Always start adjacent boxes with default color/colormap so they
        # don't inherit whatever the primary box had been set to.
        model.halo_layer.color_mode = "mvir"
        model.halo_layer.colormap = "viridis"
        model.galaxy_layer.color_mode = "structure"
        model.galaxy_layer.colormap = "plasma"
        model.set_snapshot(model.snap_count - 1)
        model.visible = True
        self._update_labels()
        self._plotter.render()
        for cb in self._on_model_change:
            cb()
        return True, None

visage.scene.camera

visage.scene.camera.CameraController

High-level camera operations on top of a PyVista Plotter.

All coordinate arguments are in Mpc/h, matching the simulation box units.

Source code in visage/scene/camera.py

class CameraController:
    """High-level camera operations on top of a PyVista Plotter.

    All coordinate arguments are in Mpc/h, matching the simulation box units.
    """

    def __init__(self, plotter: pv.Plotter, box_size: float = 62.5) -> None:
        self._pl = plotter
        self._box_size = box_size
        self._halo_index: NearestHaloIndex = NearestHaloIndex()
        self._galaxy_positions: np.ndarray | None = None
        self._indicator_actor = None
        # The selection ring (target + label + group members).
        from visage.scene.selection_marker import SelectionMarker

        self.marker = SelectionMarker(plotter)
        self._member_actors: list = (
            []
        )  # splats for FOF group members (one per regime colour)
        self._selected_actors: list = (
            []
        )  # splats for selected galaxy: [white border, regime fill]
        self._central_actor = (
            None  # thin white outline marking the FOF central
        )
        self._group_ring_actor = None  # red ring sized to enclose the group
        # Star field — a backdrop drawn behind an isolated object.  Kept
        # anchored to the camera, so it needs the state below to follow it.
        self._sparkle_actors: list = []
        self._sparkle_state: dict | None = None
        self._sparkle_obs: int | None = None
        # Lightcone mode: frame the actual point-cloud bounds instead of a
        # cubic box (a cone is not a box). None => normal box framing.
        self._lc_bounds: np.ndarray | None = None

    def set_lightcone_bounds(self, bounds: np.ndarray | None) -> None:
        """Enable lightcone framing: bounds = [[xmin,ymin,zmin],[xmax,ymax,zmax]]."""
        self._lc_bounds = (
            None if bounds is None else np.asarray(bounds, dtype=float)
        )

    # ------------------------------------------------------------------
    # Index updates (called by Scene on snapshot change)
    # ------------------------------------------------------------------

    def update_halo_index(self, positions: np.ndarray, tree=None) -> None:
        if len(positions) > 0:
            self._halo_index.update(positions, tree=tree)

    def update_galaxy_positions(self, positions: np.ndarray) -> None:
        self._galaxy_positions = positions

    # ------------------------------------------------------------------
    # Sparkle field (Isolate view)
    # ------------------------------------------------------------------

    @property
    def has_sparkles(self) -> bool:
        return bool(self._sparkle_actors)

    # How far out the backdrop sits, in units of the isolate framing
    # radius: far enough to read as sky, near enough to stay inside the
    # far clipping plane.
    _STAR_DISTANCE = 25.0

    def clear_sparkles(self) -> None:
        if self._sparkle_obs is not None:
            try:
                self._pl.renderer.GetActiveCamera().RemoveObserver(
                    self._sparkle_obs
                )
            except Exception:
                pass
            self._sparkle_obs = None
        self._sparkle_state = None
        for a in self._sparkle_actors:
            self._pl.remove_actor(a, render=False)
        self._sparkle_actors.clear()

    # Stellar colours by spectral class: O/B blue, A white-blue, F white,
    # G yellow-white, K orange, M red.
    _STAR_COLORS = np.array(
        [
            [155, 176, 255],  # O/B  blue
            [202, 215, 255],  # A    blue-white
            [248, 247, 255],  # F    white
            [255, 244, 234],  # G    yellow-white
            [255, 204, 140],  # K    orange
            [255, 156, 84],  # M    deep orange
            [255, 108, 72],  # M    red
            [246, 74, 60],  # M    deep red
        ],
        dtype=np.uint8,
    )
    _STAR_WEIGHTS = np.array([0.06, 0.10, 0.15, 0.18, 0.17, 0.14, 0.12, 0.08])

    def add_sparkles(
        self,
        center: tuple[float, float, float],
        radius: float,
        n_stars: int = 90_000,
        seed: int = 12345,
        point_size: float = 2.0,
    ) -> None:
        """A star backdrop behind the isolated object.

        `n_stars` covers the whole sphere around the camera, of which the
        field of view shows about 3-4%, so the count is far higher than the
        number on screen at any moment (~3,000 of 90,000).

        The stars are anchored to the camera rather than scattered around
        the object: a shell centred on the object puts half of its stars
        between the viewer and the object, where they read as specks over
        the galaxies instead of sky behind them.  Anchored to the camera
        they sit at a fixed distance in every direction — the parallax-free
        backdrop a real star field has — and they are repositioned as the
        camera moves.  Round points of a fixed pixel size, so they look the
        same however the user zooms.
        """
        self.clear_sparkles()
        if radius <= 0 or n_stars <= 0:
            return
        rng = np.random.default_rng(seed)

        dirs = rng.normal(size=(n_stars, 3))
        dirs /= np.maximum(np.linalg.norm(dirs, axis=1)[:, None], 1e-12)
        # A little depth variation so the field doesn't look like a decal.
        spread = 1.0 + 0.35 * rng.random(n_stars)

        # Colours down the spectral sequence, O/B through M, with the
        # white/yellow middle most common — real stars, not a rainbow.
        colors = self._STAR_COLORS[
            rng.choice(
                len(self._STAR_COLORS), size=n_stars, p=self._STAR_WEIGHTS
            )
        ]
        cloud = pv.PolyData(np.zeros((n_stars, 3)))

        self._sparkle_state = {
            "dirs": dirs,
            "spread": spread,
            "distance": max(float(radius), 1e-6) * self._STAR_DISTANCE,
            "center": np.asarray(center, dtype=float),
            "radius": max(float(radius), 1e-6),
            "rgba": np.concatenate(
                [colors, np.full((n_stars, 1), 255, np.uint8)], axis=1
            ),
            "cloud": cloud,
            "cam_key": None,
        }
        self._reposition_sparkles()

        self._sparkle_actors.append(
            self._pl.add_mesh(
                cloud,
                scalars="rgba",
                rgb=True,
                point_size=point_size,
                render_points_as_spheres=True,
                opacity=0.9,
                # Flat: a two-pixel dot has no shading to show, and a
                # white specular highlight would wash the reds out.
                lighting=False,
                show_scalar_bar=False,
                render=False,
                reset_camera=False,
            )
        )
        try:
            self._sparkle_obs = (
                self._pl.renderer.GetActiveCamera().AddObserver(
                    "ModifiedEvent", self._on_camera_moved
                )
            )
        except Exception:
            self._sparkle_obs = None
        # The backdrop sits well beyond the isolated object, so the far
        # clipping plane has to be re-fitted or the stars are clipped away.
        self._pl.renderer.ResetCameraClippingRange()
        # Show them straight away — without this the stars only turn up on
        # the next render the view happens to do (e.g. a camera move).
        self._pl.render()

    def _camera_key(self) -> tuple:
        cam = self._pl.camera
        return (
            tuple(round(float(v), 6) for v in cam.position),
            tuple(round(float(v), 6) for v in cam.focal_point),
        )

    def _reposition_sparkles(self) -> None:
        """Park the backdrop around the camera, beyond what it is looking at.

        The distance tracks the camera-to-focal-point distance as well as
        the framing radius: a fixed distance would end up in front of the
        scene as soon as the viewer pulled back.
        """
        st = self._sparkle_state
        if st is None:
            return
        cam = np.asarray(self._pl.camera.position, dtype=float)
        focal = np.asarray(self._pl.camera.focal_point, dtype=float)
        to_focal = float(np.linalg.norm(focal - cam))
        distance = max(st["distance"], 2.5 * to_focal)
        st["cloud"].points = (
            cam + st["dirs"] * (distance * st["spread"])[:, None]
        )

        # Emissive galaxies add light rather than blocking it, so a star
        # behind one shines straight through the glow.  Rather than empty
        # the object's line of sight — which leaves a hard-edged hole in
        # the field — the stars fade out as they approach that direction,
        # so the field simply thins where the object is.
        to_obj = st["center"] - cam
        obj_dist = float(np.linalg.norm(to_obj))
        rgba = st["rgba"].copy()
        if obj_dist > 1e-9:
            axis = to_obj / obj_dist
            # Tight: the object itself, then full sky just beyond it.  A
            # wide ramp dims most of the frame at the isolate framing,
            # where the object already subtends a good part of the view.
            inner = np.arctan2(st["radius"] * 1.0, obj_dist)
            outer = np.arctan2(st["radius"] * 2.2, obj_dist)
            # Element-wise rather than a matmul: Accelerate raises a
            # spurious divide-by-zero warning for this shape, and this
            # runs on every camera move.
            cos_to_axis = (st["dirs"] * axis).sum(axis=1)
            angle = np.arccos(np.clip(cos_to_axis, -1.0, 1.0))
            fade = np.clip(
                (angle - inner) / max(outer - inner, 1e-9), 0.0, 1.0
            )
            # Smoothstep, so there is no visible edge where it starts.
            fade = fade * fade * (3.0 - 2.0 * fade)
            rgba[:, 3] = (rgba[:, 3] * fade).astype(np.uint8)
        st["cloud"].point_data["rgba"] = rgba
        st["cam_key"] = self._camera_key()
        # Keep the far plane out past the backdrop.  This fires the
        # camera's ModifiedEvent, but the key above has already been
        # updated, so the observer returns immediately instead of looping.
        try:
            self._pl.renderer.ResetCameraClippingRange()
        except Exception:
            pass

    def _on_camera_moved(self, *_args) -> None:
        # The camera fires ModifiedEvent for things that do not move it —
        # a clipping-range update during a render above all — and
        # rewriting the points mid-render costs that frame.
        st = self._sparkle_state
        if st is None:
            return
        if self._camera_key() == st["cam_key"]:
            return
        self._reposition_sparkles()

    # ------------------------------------------------------------------
    # Zoom indicators
    # ------------------------------------------------------------------

    def _clear_indicator(self, keep_marker: bool = False) -> None:
        if not keep_marker:
            self.marker.clear()
        if self._indicator_actor is None:
            return
        actors = (
            self._indicator_actor
            if isinstance(self._indicator_actor, list)
            else [self._indicator_actor]
        )
        for a in actors:
            if a is not None:
                self._pl.remove_actor(a, render=False)
        self._indicator_actor = None

    # Regime colours: cold CGM = dodger blue, hot = tomato, unknown = cyan
    _REGIME_COLORS = {0: "dodgerblue", 1: "tomato", -1: "cyan"}

    # The FOF central.  Green rather than the old gold: against the members'
    # blue and red, and over the galaxies' own warm splats, gold was the one
    # marker that disappeared.
    _CENTRAL_COLOR = "springgreen"

    # Markers are lit spheres: the diffuse falloff and specular highlight
    # give them shape, so they read as balls rather than flat discs.
    _MARKER_SHADING = dict(
        ambient=0.22, diffuse=0.85, specular=1.0, specular_power=30
    )

    def _add_marker_points(
        self,
        positions: np.ndarray,
        color: str,
        size: float,
        opacity: float,
        store: list,
    ) -> None:
        """Shaded marker spheres at `positions`."""
        cloud = pv.PolyData(np.asarray(positions, dtype=np.float64))
        store.append(
            self._pl.add_mesh(
                cloud,
                color=color,
                point_size=size,
                render_points_as_spheres=True,
                opacity=opacity,
                show_scalar_bar=False,
                render=False,
                reset_camera=False,
                **self._MARKER_SHADING,
            )
        )

    def _clear_member_indicators(self) -> None:
        for a in self._member_actors + self._selected_actors:
            self._pl.remove_actor(a, render=False)
        self._member_actors.clear()
        self._selected_actors.clear()

    def _add_group_central_indicator(self, position: np.ndarray) -> None:
        """Marker for the FOF central — appended to _member_actors so it
        clears with the rest of the group."""
        # Noticeably larger than a member (24) so a member sitting in front
        # of it still leaves a green rim showing.
        self._add_marker_points(
            [position], self._CENTRAL_COLOR, 42.0, 0.95, self._member_actors
        )

    def _add_member_indicators(
        self,
        positions: np.ndarray,
        regimes: np.ndarray | None = None,
    ) -> None:
        """Splats for FOF group members coloured by CGM/hot regime."""
        for a in self._member_actors:
            self._pl.remove_actor(a, render=False)
        self._member_actors.clear()
        if positions is None or len(positions) == 0:
            return
        positions = np.asarray(positions, dtype=np.float64)
        # Group by regime so each colour gets one draw call
        if regimes is not None and len(regimes) == len(positions):
            groups = {0: [], 1: [], -1: []}
            for i, r in enumerate(regimes):
                groups[r if r in (0, 1) else -1].append(i)
        else:
            groups = {-1: list(range(len(positions)))}
        for regime_key, idxs in groups.items():
            if not idxs:
                continue
            self._add_marker_points(
                positions[idxs],
                self._REGIME_COLORS[regime_key],
                24.0,
                0.80,
                self._member_actors,
            )

    def _add_selected_indicator(
        self,
        position: np.ndarray,
        regime: int | None = None,
        color: str | None = None,
    ) -> None:
        """Selected galaxy: white border sphere + coloured fill sphere.

        color: if given, overrides the regime-based fill colour.
        """
        for a in self._selected_actors:
            self._pl.remove_actor(a, render=False)
        self._selected_actors.clear()
        if position is None:
            return
        if color is None:
            color = self._REGIME_COLORS.get(
                regime if regime in (0, 1) else -1, "cyan"
            )
        self._add_marker_points(
            [position], color, 30.0, 0.95, self._selected_actors
        )

    @property
    def has_member_indicators(self) -> bool:
        return bool(self._member_actors or self._selected_actors)

    # ---- White central marker -----------------------------------------

    def _clear_central_indicator(self) -> None:
        if self._central_actor is None:
            return
        self._pl.remove_actor(self._central_actor, render=False)
        self._central_actor = None

    def _add_central_indicator(
        self, center: tuple[float, float, float]
    ) -> None:
        """Thin white screen-space ring marking the FOF central."""
        self._clear_central_indicator()
        cloud = pv.PolyData(np.array([center], dtype=np.float64))
        self._central_actor = self._pl.add_mesh(
            cloud,
            color="white",
            point_size=22.0,
            render_points_as_spheres=False,
            opacity=0.95,
            show_scalar_bar=False,
            render=False,
            reset_camera=False,
        )

    # ---- Group-sized red circle ---------------------------------------

    def _clear_group_ring(self) -> None:
        if self._group_ring_actor is None:
            return
        self._pl.remove_actor(self._group_ring_actor, render=False)
        self._group_ring_actor = None

    def _add_group_ring(
        self,
        center: tuple[float, float, float],
        radius: float,
    ) -> None:
        """Red ring face-on to the camera, sized to enclose the whole group."""
        self._clear_group_ring()
        if radius <= 0:
            return
        c = np.array(center, dtype=np.float64)
        cam = np.array(self._pl.camera.position, dtype=np.float64)
        view = cam - c
        norm = np.linalg.norm(view)
        if norm < 1e-10:
            return
        view /= norm
        up = np.array([0.0, 1.0, 0.0])
        if abs(np.dot(view, up)) > 0.99:
            up = np.array([1.0, 0.0, 0.0])
        right = np.cross(view, up)
        right /= np.linalg.norm(right)
        up_perp = np.cross(right, view)
        up_perp /= np.linalg.norm(up_perp)

        theta = np.linspace(0, 2 * np.pi, 128, endpoint=False)
        pts = c + radius * (
            np.outer(np.cos(theta), right) + np.outer(np.sin(theta), up_perp)
        )
        pts = np.vstack([pts, pts[0]])  # close the ring
        circle = pv.lines_from_points(pts)
        self._group_ring_actor = self._pl.add_mesh(
            circle,
            color="red",
            line_width=2.0,
            opacity=0.9,
            render=False,
            reset_camera=False,
        )

    def _add_box_indicator(
        self,
        xmin: float,
        xmax: float,
        ymin: float,
        ymax: float,
        zmin: float,
        zmax: float,
    ) -> None:
        self._clear_indicator()
        box = pv.Box(bounds=(xmin, xmax, ymin, ymax, zmin, zmax))
        self._indicator_actor = self._pl.add_mesh(
            box.extract_all_edges(),
            color=_INDICATOR_COLOR,
            opacity=_INDICATOR_OPACITY,
            line_width=_INDICATOR_WIDTH,
            render_lines_as_tubes=False,
        )

    def _add_sphere_indicator(
        self,
        center: tuple[float, float, float],
        radius: float,
    ) -> None:
        """Three orthogonal great-circle rings — minimal wireframe look,
        much sparser than a full sphere mesh."""
        self._clear_indicator()
        cx, cy, cz = center
        t = np.linspace(0.0, 2.0 * np.pi, 64, dtype=np.float64)
        c, s = np.cos(t) * radius, np.sin(t) * radius
        # 1 equator (XY) + 4 meridians around the polar (Z) axis at
        # azimuths 0°, 45°, 90°, 135°.
        rings = [np.column_stack([cx + c, cy + s, np.full_like(c, cz)])]
        for deg in (0.0, 45.0, 90.0, 135.0):
            ca, sa = np.cos(np.deg2rad(deg)), np.sin(np.deg2rad(deg))
            rings.append(np.column_stack([cx + c * ca, cy + c * sa, cz + s]))
        # Build one PolyData with 5 closed polylines.
        all_pts = np.vstack(rings)
        n = len(t)
        n_rings = len(rings)
        lines = []
        for i in range(n_rings):
            lines.append(n + 1)
            lines.extend([i * n + j for j in range(n)])
            lines.append(i * n)  # close the loop
        poly = pv.PolyData(all_pts)
        # pv.PolyData(points) auto-creates a verts cell per point, which
        # VTK renders as visible point markers. Strip them so we get
        # pure lines.
        poly.verts = np.empty(0, dtype=np.int64)
        poly.lines = np.array(lines, dtype=np.int64)
        self._indicator_actor = self._pl.add_mesh(
            poly,
            color=_INDICATOR_COLOR,
            opacity=_INDICATOR_OPACITY,
            line_width=_INDICATOR_WIDTH,
            style="wireframe",
            render_points_as_spheres=False,
            point_size=0,
        )
        # Belt-and-braces: explicitly turn off vertex rendering on the
        # actor's property so VTK never draws point markers at the ring
        # vertices.
        try:
            self._indicator_actor.GetProperty().SetRenderPointsAsSpheres(False)
            self._indicator_actor.GetProperty().SetVertexVisibility(False)
            self._indicator_actor.GetProperty().SetPointSize(0)
        except Exception:
            pass

    def _add_point_indicator(
        self,
        center: tuple[float, float, float],
        radius: float = 1.0,
    ) -> None:
        """Small sphere marking a specific halo or galaxy position."""
        self._clear_indicator()
        sphere = pv.Sphere(
            radius=radius,
            center=center,
            theta_resolution=12,
            phi_resolution=12,
        )
        self._indicator_actor = self._pl.add_mesh(
            sphere,
            color=_INDICATOR_COLOR,
            opacity=_INDICATOR_OPACITY,
            style="wireframe",
            line_width=_INDICATOR_WIDTH,
        )

    # ------------------------------------------------------------------
    # Navigation
    # ------------------------------------------------------------------

    def reset(self) -> None:
        """Fit the full simulation box in view, centred on the box midpoint.

        In lightcone mode, frame the actual point-cloud bounds (a cone is not
        a cube) with an isometric 3/4 view instead of the box math."""
        self._clear_indicator()
        if self._lc_bounds is not None:
            # Lay the cone out horizontally, centred in the viewport: look
            # straight down -Z so the long radial (X) axis is horizontal and
            # the observer→distance direction reads left→right, then fit.
            b = self._lc_bounds
            cx, cy, cz = (b[0] + b[1]) / 2.0
            span = float(np.max(b[1] - b[0])) or 1.0
            self._pl.camera_position = [
                (cx, cy, cz + span * 2.0),  # camera on +Z, looking down -Z
                (cx, cy, cz),  # focus at cone centre
                (0.0, 1.0, 0.0),  # +Y up  → X axis horizontal
            ]
            self._pl.reset_camera()
            # reset_camera fits the bounding SPHERE (radius ~ span/2) to the
            # viewport's SHORTER axis. A lightcone is long and thin, so that
            # leaves it a tiny sliver ringed by empty margin. Zoom so the
            # cone's widest extent fills the viewport WIDTH instead — properly
            # zoomed-in and centred, still showing the cone end to end.
            try:
                w, h = self._pl.window_size
                if w and h:
                    self._pl.camera.zoom(max(1.0, (w / h) * 0.9))
                else:
                    self._pl.camera.zoom(1.5)
            except Exception:
                self._pl.camera.zoom(1.5)
            self._pl.renderer.ResetCameraClippingRange()
            return
        half = self._box_size / 2.0
        self._pl.camera_position = [
            (half, half, self._box_size * 2.2),
            (half, half, half),
            (0.0, 1.0, 0.0),
        ]

    def focus_on_boxes(
        self,
        regions: list[tuple[float, float, float, float]],
    ) -> None:
        """Frame one or more simulation boxes.

        regions: list of (offset_x, offset_y, offset_z, box_size).
        A single box with zero offset reproduces the same view as reset().
        """
        if not regions:
            return
        xmin = min(ox for ox, oy, oz, bs in regions)
        xmax = max(ox + bs for ox, oy, oz, bs in regions)
        ymin = min(oy for ox, oy, oz, bs in regions)
        ymax = max(oy + bs for ox, oy, oz, bs in regions)
        zmin = min(oz for ox, oy, oz, bs in regions)
        zmax = max(oz + bs for ox, oy, oz, bs in regions)
        cx = (xmin + xmax) / 2
        cy = (ymin + ymax) / 2
        cz = (zmin + zmax) / 2
        span = max(xmax - xmin, ymax - ymin, zmax - zmin)
        self._pl.camera_position = [
            (cx, cy, cz + span * 1.7),
            (cx, cy, cz),
            (0.0, 1.0, 0.0),
        ]

    def go_to_lightcone_observer(self) -> None:
        """Place the camera AT the observer (coordinate origin) looking
        outward along the cone axis, with the sky spread horizontally — as if
        standing at the observer and gazing out into the lightcone."""
        self._clear_indicator()
        if self._lc_bounds is None:
            return
        b = self._lc_bounds
        center = (b[0] + b[1]) / 2.0
        n = float(np.linalg.norm(center))
        if n < 1e-9:
            return
        self._pl.camera.position = (0.0, 0.0, 0.0)
        # Focal point AT the cone centre (not a unit step out): this keeps the
        # focal distance ~half the cone depth, so the trackball orbits around
        # the cone and VTK's clipping-range heuristic (which scales with the
        # focal distance) stays sane. A focal point ~1 Mpc from the camera —
        # as before — made the near plane collapse and geometry pop in/out of
        # the far clip whenever the user rotated or panned.
        self._pl.camera.focal_point = tuple(center)  # look outward at the cone
        # +Z up puts the (wider) Y sky-extent left→right, so the cone opens
        # out horizontally in front of the viewer.
        self._pl.camera.up = (0.0, 0.0, 1.0)
        self._pl.renderer.ResetCameraClippingRange()

    def go_to_box_center(
        self,
        offset: tuple[float, float, float] = (0.0, 0.0, 0.0),
        box_size: float | None = None,
    ) -> None:
        """Place the camera AT the active box centre, looking along +z."""
        self._clear_indicator()
        bs = box_size if box_size is not None else self._box_size
        ox, oy, oz = offset
        half = bs / 2.0
        cx, cy, cz = ox + half, oy + half, oz + half
        self._pl.camera.focal_point = (cx, cy, cz + 1.0)
        self._pl.camera.position = (cx, cy, cz)
        self._pl.camera.up = (0.0, 1.0, 0.0)

    # ------------------------------------------------------------------
    # Keyboard fly movement
    # ------------------------------------------------------------------

    def fly(self, direction: str, step_frac: float = 0.012) -> None:
        """Translate the camera (and its focal point) one step in a view-
        relative direction. Because both move together this is a true fly —
        it carries the camera through the box centre and out the far side,
        unlike trackball dolly/orbit which stall at the focal point.

        direction: 'forward' | 'back' | 'left' | 'right' | 'up' | 'down'.
        """
        cam = self._pl.camera
        pos = np.array(cam.position, dtype=np.float64)
        focal = np.array(cam.focal_point, dtype=np.float64)

        view = focal - pos
        nv = np.linalg.norm(view)
        if nv < 1e-9:
            return
        view /= nv
        up = np.array(cam.up, dtype=np.float64)
        nu = np.linalg.norm(up)
        up = up / nu if nu > 1e-9 else np.array([0.0, 1.0, 0.0])
        right = np.cross(view, up)
        nr = np.linalg.norm(right)
        right = right / nr if nr > 1e-9 else np.array([1.0, 0.0, 0.0])
        up = np.cross(right, view)  # re-orthonormalise

        basis = {
            "forward": view,
            "back": -view,
            "right": right,
            "left": -right,
            "up": up,
            "down": -up,
        }
        d = basis.get(direction)
        if d is None:
            return

        delta = d * (self._box_size * float(step_frac))
        cam.position = tuple(pos + delta)
        cam.focal_point = tuple(focal + delta)

    def go_to_coords(
        self,
        x: float,
        y: float,
        z: float,
        distance: float = 5.0,
    ) -> None:
        """Point camera at (x, y, z) from a given standoff distance, and
        draw a wireframe sphere of radius `distance` at the target so the
        focus region is visible (matches the Box-mode wireframe convention)."""
        self._pl.camera.focal_point = (x, y, z)
        self._pl.camera.position = (x, y, z + distance)
        self._pl.camera.up = (0.0, 1.0, 0.0)
        self._add_sphere_indicator((x, y, z), distance)

    def go_to_halo(self, halo_idx: int, distance: float = 5.0) -> None:
        """Fly to the halo at halo_idx and mark it with a red circle."""
        pos = self._halo_index.position_of(halo_idx)
        cx, cy, cz = float(pos[0]), float(pos[1]), float(pos[2])
        self._pl.camera.focal_point = (cx, cy, cz)
        self._pl.camera.position = (cx, cy, cz + distance)
        self._pl.camera.up = (0.0, 1.0, 0.0)
        self._add_circle_indicator((cx, cy, cz), distance * 0.015)

    def go_to_nearest_halo(
        self,
        x: float,
        y: float,
        z: float,
        distance: float = 5.0,
    ) -> int:
        idx = self._halo_index.nearest((x, y, z))
        self.go_to_halo(idx, distance)
        return idx

    def go_to_galaxy(self, galaxy_idx: int, radius: float = 3.0) -> tuple:
        """Fly to galaxy_idx; camera sits on the sphere surface looking inward.

        A red circle marks the galaxy position. No sphere shown.
        Returns (cx, cy, cz) so callers can apply a focus sphere.
        """
        if self._galaxy_positions is None or len(self._galaxy_positions) == 0:
            return (0.0, 0.0, 0.0)
        pos = self._galaxy_positions[galaxy_idx]
        cx, cy, cz = float(pos[0]), float(pos[1]), float(pos[2])

        self._pl.camera.focal_point = (cx, cy, cz)
        self._pl.camera.position = (cx, cy, cz + radius)
        self._pl.camera.up = (0.0, 1.0, 0.0)

        self._add_circle_indicator((cx, cy, cz), radius * 0.015)
        return (cx, cy, cz)

    def _add_circle_indicator(
        self,
        center: tuple[float, float, float],
        radius: float,
    ) -> None:
        """The selection ring on the target (fixed size on screen)."""
        self.marker.set_target(center)

    def zoom_to_radius(
        self,
        center: tuple[float, float, float],
        radius: float,
    ) -> None:
        """Frame a sphere of the given radius and draw a wireframe sphere indicator."""
        cx, cy, cz = center
        fov_rad = np.deg2rad(self._pl.camera.view_angle)
        distance = radius / np.tan(fov_rad / 2.0) * 1.2
        self._pl.camera.focal_point = center
        self._pl.camera.position = (cx, cy, cz + distance)
        self._pl.camera.up = (0.0, 1.0, 0.0)
        self._add_sphere_indicator(center, radius)

    def zoom_to_box(
        self,
        xmin: float,
        xmax: float,
        ymin: float,
        ymax: float,
        zmin: float,
        zmax: float,
    ) -> None:
        """Frame an axis-aligned sub-box and draw a wireframe box indicator."""
        cx = (xmin + xmax) / 2.0
        cy = (ymin + ymax) / 2.0
        cz = (zmin + zmax) / 2.0
        radius = max(xmax - xmin, ymax - ymin, zmax - zmin) / 2.0
        fov_rad = np.deg2rad(self._pl.camera.view_angle)
        distance = radius / np.tan(fov_rad / 2.0) * 1.2
        self._pl.camera.focal_point = (cx, cy, cz)
        self._pl.camera.position = (cx, cy, cz + distance)
        self._pl.camera.up = (0.0, 1.0, 0.0)
        self._add_box_indicator(xmin, xmax, ymin, ymax, zmin, zmax)

add_sparkles(center, radius, n_stars=90000, seed=12345, point_size=2.0)

A star backdrop behind the isolated object.

n_stars covers the whole sphere around the camera, of which the field of view shows about 3-4%, so the count is far higher than the number on screen at any moment (~3,000 of 90,000).

The stars are anchored to the camera rather than scattered around the object: a shell centred on the object puts half of its stars between the viewer and the object, where they read as specks over the galaxies instead of sky behind them. Anchored to the camera they sit at a fixed distance in every direction — the parallax-free backdrop a real star field has — and they are repositioned as the camera moves. Round points of a fixed pixel size, so they look the same however the user zooms.

Source code in visage/scene/camera.py

def add_sparkles(
    self,
    center: tuple[float, float, float],
    radius: float,
    n_stars: int = 90_000,
    seed: int = 12345,
    point_size: float = 2.0,
) -> None:
    """A star backdrop behind the isolated object.

    `n_stars` covers the whole sphere around the camera, of which the
    field of view shows about 3-4%, so the count is far higher than the
    number on screen at any moment (~3,000 of 90,000).

    The stars are anchored to the camera rather than scattered around
    the object: a shell centred on the object puts half of its stars
    between the viewer and the object, where they read as specks over
    the galaxies instead of sky behind them.  Anchored to the camera
    they sit at a fixed distance in every direction — the parallax-free
    backdrop a real star field has — and they are repositioned as the
    camera moves.  Round points of a fixed pixel size, so they look the
    same however the user zooms.
    """
    self.clear_sparkles()
    if radius <= 0 or n_stars <= 0:
        return
    rng = np.random.default_rng(seed)

    dirs = rng.normal(size=(n_stars, 3))
    dirs /= np.maximum(np.linalg.norm(dirs, axis=1)[:, None], 1e-12)
    # A little depth variation so the field doesn't look like a decal.
    spread = 1.0 + 0.35 * rng.random(n_stars)

    # Colours down the spectral sequence, O/B through M, with the
    # white/yellow middle most common — real stars, not a rainbow.
    colors = self._STAR_COLORS[
        rng.choice(
            len(self._STAR_COLORS), size=n_stars, p=self._STAR_WEIGHTS
        )
    ]
    cloud = pv.PolyData(np.zeros((n_stars, 3)))

    self._sparkle_state = {
        "dirs": dirs,
        "spread": spread,
        "distance": max(float(radius), 1e-6) * self._STAR_DISTANCE,
        "center": np.asarray(center, dtype=float),
        "radius": max(float(radius), 1e-6),
        "rgba": np.concatenate(
            [colors, np.full((n_stars, 1), 255, np.uint8)], axis=1
        ),
        "cloud": cloud,
        "cam_key": None,
    }
    self._reposition_sparkles()

    self._sparkle_actors.append(
        self._pl.add_mesh(
            cloud,
            scalars="rgba",
            rgb=True,
            point_size=point_size,
            render_points_as_spheres=True,
            opacity=0.9,
            # Flat: a two-pixel dot has no shading to show, and a
            # white specular highlight would wash the reds out.
            lighting=False,
            show_scalar_bar=False,
            render=False,
            reset_camera=False,
        )
    )
    try:
        self._sparkle_obs = (
            self._pl.renderer.GetActiveCamera().AddObserver(
                "ModifiedEvent", self._on_camera_moved
            )
        )
    except Exception:
        self._sparkle_obs = None
    # The backdrop sits well beyond the isolated object, so the far
    # clipping plane has to be re-fitted or the stars are clipped away.
    self._pl.renderer.ResetCameraClippingRange()
    # Show them straight away — without this the stars only turn up on
    # the next render the view happens to do (e.g. a camera move).
    self._pl.render()

fly(direction, step_frac=0.012)

Translate the camera (and its focal point) one step in a view- relative direction. Because both move together this is a true fly — it carries the camera through the box centre and out the far side, unlike trackball dolly/orbit which stall at the focal point.

direction: 'forward' | 'back' | 'left' | 'right' | 'up' | 'down'.

Source code in visage/scene/camera.py

def fly(self, direction: str, step_frac: float = 0.012) -> None:
    """Translate the camera (and its focal point) one step in a view-
    relative direction. Because both move together this is a true fly —
    it carries the camera through the box centre and out the far side,
    unlike trackball dolly/orbit which stall at the focal point.

    direction: 'forward' | 'back' | 'left' | 'right' | 'up' | 'down'.
    """
    cam = self._pl.camera
    pos = np.array(cam.position, dtype=np.float64)
    focal = np.array(cam.focal_point, dtype=np.float64)

    view = focal - pos
    nv = np.linalg.norm(view)
    if nv < 1e-9:
        return
    view /= nv
    up = np.array(cam.up, dtype=np.float64)
    nu = np.linalg.norm(up)
    up = up / nu if nu > 1e-9 else np.array([0.0, 1.0, 0.0])
    right = np.cross(view, up)
    nr = np.linalg.norm(right)
    right = right / nr if nr > 1e-9 else np.array([1.0, 0.0, 0.0])
    up = np.cross(right, view)  # re-orthonormalise

    basis = {
        "forward": view,
        "back": -view,
        "right": right,
        "left": -right,
        "up": up,
        "down": -up,
    }
    d = basis.get(direction)
    if d is None:
        return

    delta = d * (self._box_size * float(step_frac))
    cam.position = tuple(pos + delta)
    cam.focal_point = tuple(focal + delta)

focus_on_boxes(regions)

Frame one or more simulation boxes.

regions: list of (offset_x, offset_y, offset_z, box_size). A single box with zero offset reproduces the same view as reset().

Source code in visage/scene/camera.py

def focus_on_boxes(
    self,
    regions: list[tuple[float, float, float, float]],
) -> None:
    """Frame one or more simulation boxes.

    regions: list of (offset_x, offset_y, offset_z, box_size).
    A single box with zero offset reproduces the same view as reset().
    """
    if not regions:
        return
    xmin = min(ox for ox, oy, oz, bs in regions)
    xmax = max(ox + bs for ox, oy, oz, bs in regions)
    ymin = min(oy for ox, oy, oz, bs in regions)
    ymax = max(oy + bs for ox, oy, oz, bs in regions)
    zmin = min(oz for ox, oy, oz, bs in regions)
    zmax = max(oz + bs for ox, oy, oz, bs in regions)
    cx = (xmin + xmax) / 2
    cy = (ymin + ymax) / 2
    cz = (zmin + zmax) / 2
    span = max(xmax - xmin, ymax - ymin, zmax - zmin)
    self._pl.camera_position = [
        (cx, cy, cz + span * 1.7),
        (cx, cy, cz),
        (0.0, 1.0, 0.0),
    ]

go_to_box_center(offset=(0.0, 0.0, 0.0), box_size=None)

Place the camera AT the active box centre, looking along +z.

Source code in visage/scene/camera.py

def go_to_box_center(
    self,
    offset: tuple[float, float, float] = (0.0, 0.0, 0.0),
    box_size: float | None = None,
) -> None:
    """Place the camera AT the active box centre, looking along +z."""
    self._clear_indicator()
    bs = box_size if box_size is not None else self._box_size
    ox, oy, oz = offset
    half = bs / 2.0
    cx, cy, cz = ox + half, oy + half, oz + half
    self._pl.camera.focal_point = (cx, cy, cz + 1.0)
    self._pl.camera.position = (cx, cy, cz)
    self._pl.camera.up = (0.0, 1.0, 0.0)

go_to_coords(x, y, z, distance=5.0)

Point camera at (x, y, z) from a given standoff distance, and draw a wireframe sphere of radius distance at the target so the focus region is visible (matches the Box-mode wireframe convention).

Source code in visage/scene/camera.py

def go_to_coords(
    self,
    x: float,
    y: float,
    z: float,
    distance: float = 5.0,
) -> None:
    """Point camera at (x, y, z) from a given standoff distance, and
    draw a wireframe sphere of radius `distance` at the target so the
    focus region is visible (matches the Box-mode wireframe convention)."""
    self._pl.camera.focal_point = (x, y, z)
    self._pl.camera.position = (x, y, z + distance)
    self._pl.camera.up = (0.0, 1.0, 0.0)
    self._add_sphere_indicator((x, y, z), distance)

go_to_galaxy(galaxy_idx, radius=3.0)

Fly to galaxy_idx; camera sits on the sphere surface looking inward.

A red circle marks the galaxy position. No sphere shown. Returns (cx, cy, cz) so callers can apply a focus sphere.

Source code in visage/scene/camera.py

def go_to_galaxy(self, galaxy_idx: int, radius: float = 3.0) -> tuple:
    """Fly to galaxy_idx; camera sits on the sphere surface looking inward.

    A red circle marks the galaxy position. No sphere shown.
    Returns (cx, cy, cz) so callers can apply a focus sphere.
    """
    if self._galaxy_positions is None or len(self._galaxy_positions) == 0:
        return (0.0, 0.0, 0.0)
    pos = self._galaxy_positions[galaxy_idx]
    cx, cy, cz = float(pos[0]), float(pos[1]), float(pos[2])

    self._pl.camera.focal_point = (cx, cy, cz)
    self._pl.camera.position = (cx, cy, cz + radius)
    self._pl.camera.up = (0.0, 1.0, 0.0)

    self._add_circle_indicator((cx, cy, cz), radius * 0.015)
    return (cx, cy, cz)

go_to_halo(halo_idx, distance=5.0)

Fly to the halo at halo_idx and mark it with a red circle.

Source code in visage/scene/camera.py

def go_to_halo(self, halo_idx: int, distance: float = 5.0) -> None:
    """Fly to the halo at halo_idx and mark it with a red circle."""
    pos = self._halo_index.position_of(halo_idx)
    cx, cy, cz = float(pos[0]), float(pos[1]), float(pos[2])
    self._pl.camera.focal_point = (cx, cy, cz)
    self._pl.camera.position = (cx, cy, cz + distance)
    self._pl.camera.up = (0.0, 1.0, 0.0)
    self._add_circle_indicator((cx, cy, cz), distance * 0.015)

go_to_lightcone_observer()

Place the camera AT the observer (coordinate origin) looking outward along the cone axis, with the sky spread horizontally — as if standing at the observer and gazing out into the lightcone.

Source code in visage/scene/camera.py

def go_to_lightcone_observer(self) -> None:
    """Place the camera AT the observer (coordinate origin) looking
    outward along the cone axis, with the sky spread horizontally — as if
    standing at the observer and gazing out into the lightcone."""
    self._clear_indicator()
    if self._lc_bounds is None:
        return
    b = self._lc_bounds
    center = (b[0] + b[1]) / 2.0
    n = float(np.linalg.norm(center))
    if n < 1e-9:
        return
    self._pl.camera.position = (0.0, 0.0, 0.0)
    # Focal point AT the cone centre (not a unit step out): this keeps the
    # focal distance ~half the cone depth, so the trackball orbits around
    # the cone and VTK's clipping-range heuristic (which scales with the
    # focal distance) stays sane. A focal point ~1 Mpc from the camera —
    # as before — made the near plane collapse and geometry pop in/out of
    # the far clip whenever the user rotated or panned.
    self._pl.camera.focal_point = tuple(center)  # look outward at the cone
    # +Z up puts the (wider) Y sky-extent left→right, so the cone opens
    # out horizontally in front of the viewer.
    self._pl.camera.up = (0.0, 0.0, 1.0)
    self._pl.renderer.ResetCameraClippingRange()

reset()

Fit the full simulation box in view, centred on the box midpoint.

In lightcone mode, frame the actual point-cloud bounds (a cone is not a cube) with an isometric 3/4 view instead of the box math.

Source code in visage/scene/camera.py

def reset(self) -> None:
    """Fit the full simulation box in view, centred on the box midpoint.

    In lightcone mode, frame the actual point-cloud bounds (a cone is not
    a cube) with an isometric 3/4 view instead of the box math."""
    self._clear_indicator()
    if self._lc_bounds is not None:
        # Lay the cone out horizontally, centred in the viewport: look
        # straight down -Z so the long radial (X) axis is horizontal and
        # the observer→distance direction reads left→right, then fit.
        b = self._lc_bounds
        cx, cy, cz = (b[0] + b[1]) / 2.0
        span = float(np.max(b[1] - b[0])) or 1.0
        self._pl.camera_position = [
            (cx, cy, cz + span * 2.0),  # camera on +Z, looking down -Z
            (cx, cy, cz),  # focus at cone centre
            (0.0, 1.0, 0.0),  # +Y up  → X axis horizontal
        ]
        self._pl.reset_camera()
        # reset_camera fits the bounding SPHERE (radius ~ span/2) to the
        # viewport's SHORTER axis. A lightcone is long and thin, so that
        # leaves it a tiny sliver ringed by empty margin. Zoom so the
        # cone's widest extent fills the viewport WIDTH instead — properly
        # zoomed-in and centred, still showing the cone end to end.
        try:
            w, h = self._pl.window_size
            if w and h:
                self._pl.camera.zoom(max(1.0, (w / h) * 0.9))
            else:
                self._pl.camera.zoom(1.5)
        except Exception:
            self._pl.camera.zoom(1.5)
        self._pl.renderer.ResetCameraClippingRange()
        return
    half = self._box_size / 2.0
    self._pl.camera_position = [
        (half, half, self._box_size * 2.2),
        (half, half, half),
        (0.0, 1.0, 0.0),
    ]

set_lightcone_bounds(bounds)

Enable lightcone framing: bounds = [[xmin,ymin,zmin],[xmax,ymax,zmax]].

Source code in visage/scene/camera.py

def set_lightcone_bounds(self, bounds: np.ndarray | None) -> None:
    """Enable lightcone framing: bounds = [[xmin,ymin,zmin],[xmax,ymax,zmax]]."""
    self._lc_bounds = (
        None if bounds is None else np.asarray(bounds, dtype=float)
    )

zoom_to_box(xmin, xmax, ymin, ymax, zmin, zmax)

Frame an axis-aligned sub-box and draw a wireframe box indicator.

Source code in visage/scene/camera.py

def zoom_to_box(
    self,
    xmin: float,
    xmax: float,
    ymin: float,
    ymax: float,
    zmin: float,
    zmax: float,
) -> None:
    """Frame an axis-aligned sub-box and draw a wireframe box indicator."""
    cx = (xmin + xmax) / 2.0
    cy = (ymin + ymax) / 2.0
    cz = (zmin + zmax) / 2.0
    radius = max(xmax - xmin, ymax - ymin, zmax - zmin) / 2.0
    fov_rad = np.deg2rad(self._pl.camera.view_angle)
    distance = radius / np.tan(fov_rad / 2.0) * 1.2
    self._pl.camera.focal_point = (cx, cy, cz)
    self._pl.camera.position = (cx, cy, cz + distance)
    self._pl.camera.up = (0.0, 1.0, 0.0)
    self._add_box_indicator(xmin, xmax, ymin, ymax, zmin, zmax)

zoom_to_radius(center, radius)

Frame a sphere of the given radius and draw a wireframe sphere indicator.

Source code in visage/scene/camera.py

def zoom_to_radius(
    self,
    center: tuple[float, float, float],
    radius: float,
) -> None:
    """Frame a sphere of the given radius and draw a wireframe sphere indicator."""
    cx, cy, cz = center
    fov_rad = np.deg2rad(self._pl.camera.view_angle)
    distance = radius / np.tan(fov_rad / 2.0) * 1.2
    self._pl.camera.focal_point = center
    self._pl.camera.position = (cx, cy, cz + distance)
    self._pl.camera.up = (0.0, 1.0, 0.0)
    self._add_sphere_indicator(center, radius)

visage.scene.halo_layer

visage.scene.halo_layer.HaloLayer

Manages the halo point-cloud actor(s) inside a PyVista Plotter.

Source code in visage/scene/halo_layer.py

class HaloLayer:
    """Manages the halo point-cloud actor(s) inside a PyVista Plotter."""

    def __init__(
        self,
        plotter: pv.Plotter,
        color_mode: ColorMode = "mvir",
        colormap: str = "viridis",
        opacity: float = 0.05,
        visible: bool = True,
    ) -> None:
        self._pl = plotter
        self._color_mode: ColorMode = color_mode
        self._colormap = colormap
        self._opacity = opacity
        self._visible = visible
        self._actors: list = []
        # (actor, floor, multiplier) so the opacity slider can move the
        # layered splats in place instead of rebuilding them.
        self._opacity_terms: list[tuple] = []
        self._snapshot: HaloSnapshot | None = None
        self._focus_mask: np.ndarray | None = None
        self._filter_mask: np.ndarray | None = None
        self._slice_mask: np.ndarray | None = None  # lightcone redshift cut
        self._offset: np.ndarray = np.zeros(3, dtype=np.float32)

    # ------------------------------------------------------------------
    # Public API
    # ------------------------------------------------------------------

    @property
    def visible(self) -> bool:
        return self._visible

    @visible.setter
    def visible(self, value: bool) -> None:
        self._visible = value
        for actor in self._actors:
            actor.SetVisibility(value)
        self._pl.render()

    @property
    def opacity(self) -> float:
        return self._opacity

    @opacity.setter
    def opacity(self, value: float) -> None:
        self._opacity = float(value)
        if self._snapshot is None:
            return
        if not self._opacity_terms:
            self._redraw()
            return
        for actor, floor, mul in self._opacity_terms:
            actor.prop.opacity = max(floor, self._opacity * mul)

    @property
    def color_mode(self) -> ColorMode:
        return self._color_mode

    @color_mode.setter
    def color_mode(self, value: ColorMode) -> None:
        self._color_mode = value
        if self._snapshot is not None:
            self._redraw()

    @property
    def colormap(self) -> str:
        return self._colormap

    @colormap.setter
    def colormap(self, value: str) -> None:
        self._colormap = value
        if self._snapshot is not None:
            self._redraw()

    def set_offset(self, offset: np.ndarray) -> None:
        self._offset = np.asarray(offset, dtype=np.float32)
        if self._snapshot is not None:
            self._redraw()

    def update(self, snapshot: HaloSnapshot) -> None:
        self._snapshot = snapshot
        self._redraw()

    def set_mask(self, mask: np.ndarray | None) -> None:
        self.set_focus_mask(mask)

    def set_focus_mask(self, mask: np.ndarray | None) -> None:
        self._focus_mask = mask
        if self._snapshot is not None:
            self._redraw()

    def set_filter_mask(self, mask: np.ndarray | None) -> None:
        self._filter_mask = mask
        if self._snapshot is not None:
            self._redraw()

    def set_slice_mask(self, mask: np.ndarray | None) -> None:
        """Lightcone redshift/time cut keep-mask (True = shown)."""
        self._slice_mask = mask
        if self._snapshot is not None:
            self._redraw()

    def _combined_mask(self) -> np.ndarray | None:
        masks = [
            m
            for m in (self._focus_mask, self._filter_mask, self._slice_mask)
            if m is not None
        ]
        if not masks:
            return None
        n = len(masks[0])
        if any(len(m) != n for m in masks):
            return None
        out = masks[0].copy()
        for m in masks[1:]:
            out = out & m
        return out

    # ------------------------------------------------------------------
    # Internal
    # ------------------------------------------------------------------

    def _clear_actors(self) -> None:
        for actor in self._actors:
            self._pl.remove_actor(actor, render=False)
        self._actors.clear()
        self._opacity_terms.clear()

    def _redraw(self) -> None:
        snap = self._snapshot
        if snap is None or snap.count == 0:
            self._clear_actors()
            return

        # Combined focus + filter mask
        mask = self._combined_mask()
        if mask is not None and len(mask) == snap.count:
            from visage.io.halo_reader import HaloSnapshot as _HS

            snap = _HS(
                positions=snap.positions[mask],
                masses=snap.masses[mask],
                vmax=snap.vmax[mask],
                rvir=snap.rvir[mask],
                vvir=snap.vvir[mask],
                snap_num=snap.snap_num,
            )
            if snap.count == 0:
                self._clear_actors()
                return

        colors = self._compute_colors(snap)
        radii = halo_world_radii(snap.masses)

        # Full rebuild every redraw — the layered NFW-style stack has
        # three actors per halo population so the in-place fast-path
        # doesn't apply (the same trade-off galaxies make for Structure).
        self._clear_actors()
        self._render_layered(
            snap.positions + self._offset,
            colors,
            radii,
            n_total=self._snapshot.count,
        )

    # ------------------------------------------------------------------
    # Layered NFW-style halo rendering — 3 stacked gaussian splats per
    # halo at decreasing radius and increasing opacity, giving a soft
    # density-profile look (bright core → faint Rvir boundary).
    # ------------------------------------------------------------------

    # (radius_scale, opacity_floor, opacity_multiplier).  The floors kept a
    # dragged opacity slider from fading the haloes to nothing; with the
    # slider gone and the value fixed low on purpose, a floor would just
    # override it — so the layers now scale straight from the setting and
    # only their relative weighting is kept.
    _LAYERS = (
        (1.00, 0.0, 0.35),  # outer envelope ~ Rvir boundary
        (0.45, 0.0, 0.60),  # inner halo
        (0.18, 0.0, 0.95),  # dense core
    )

    # Crowding: a big box puts so many overlapping splats along each line
    # of sight that they pile up into a solid wash at the opacity a sparse
    # box needs.  Fade with halo count, from full opacity at or below the
    # reference count down to a floor — 0.05 stays 0.05 for Millennium-like
    # boxes (~30k haloes) and lands at ~0.02 for microUchuu (~390k).
    _DENSITY_REFERENCE_COUNT = 150_000
    _MIN_DENSITY_SCALE = 0.4

    @classmethod
    def density_scale(cls, n_halos: int) -> float:
        """Opacity multiplier for a snapshot holding `n_halos` haloes."""
        if n_halos <= cls._DENSITY_REFERENCE_COUNT:
            return 1.0
        return max(
            cls._MIN_DENSITY_SCALE,
            cls._DENSITY_REFERENCE_COUNT / float(n_halos),
        )

    def _render_layered(
        self,
        positions: np.ndarray,
        colors: np.ndarray,
        radii: np.ndarray,
        n_total: int | None = None,
    ) -> None:
        if len(positions) == 0:
            return
        # Crowding is a property of the box, so it is measured on the whole
        # snapshot rather than on whatever survived the current masks.
        # Measuring the drawn subset made the haloes jump from 0.02 to 0.05
        # alpha the moment a focus, isolate or filter narrowed the view.
        crowding = self.density_scale(
            len(positions) if n_total is None else n_total
        )
        for r_scale, opa_floor, opa_mul_base in self._LAYERS:
            opa_mul = opa_mul_base * crowding
            cloud = pv.PolyData(positions)
            cloud["rgba"] = scalars_to_rgba(colors, self._colormap)
            cloud["radius"] = (radii * float(r_scale)).astype(np.float32)
            actor = self._pl.add_mesh(
                cloud,
                scalars="rgba",
                rgb=True,
                style="points_gaussian",
                emissive=False,
                opacity=max(opa_floor, self._opacity * opa_mul),
                show_scalar_bar=False,
                render=False,
                reset_camera=False,
            )
            mapper = actor.mapper
            mapper.SetScaleArray("radius")
            mapper.SetScaleFactor(1.0)
            if not self._visible:
                actor.SetVisibility(False)
            self._opacity_terms.append(
                (actor, float(opa_floor), float(opa_mul))
            )
            self._actors.append(actor)

    def _compute_colors(self, snap: HaloSnapshot) -> np.ndarray:
        vmin, vmax_r = _RANGES[self._color_mode]
        if self._color_mode == "mvir":
            return normalize_log(snap.masses, vmin, vmax_r)
        if self._color_mode == "rvir":
            return normalize_log(snap.rvir, vmin, vmax_r)
        if self._color_mode == "vvir":
            return normalize_log(snap.vvir, vmin, vmax_r)
        if self._color_mode == "vmax":
            return normalize_log(np.maximum(snap.vmax, 1e-3), vmin, vmax_r)
        return normalize_log(snap.masses, *_RANGES["mvir"])

density_scale(n_halos) classmethod

Opacity multiplier for a snapshot holding n_halos haloes.

Source code in visage/scene/halo_layer.py

@classmethod
def density_scale(cls, n_halos: int) -> float:
    """Opacity multiplier for a snapshot holding `n_halos` haloes."""
    if n_halos <= cls._DENSITY_REFERENCE_COUNT:
        return 1.0
    return max(
        cls._MIN_DENSITY_SCALE,
        cls._DENSITY_REFERENCE_COUNT / float(n_halos),
    )

set_slice_mask(mask)

Lightcone redshift/time cut keep-mask (True = shown).

Source code in visage/scene/halo_layer.py

def set_slice_mask(self, mask: np.ndarray | None) -> None:
    """Lightcone redshift/time cut keep-mask (True = shown)."""
    self._slice_mask = mask
    if self._snapshot is not None:
        self._redraw()

visage.scene.galaxy_layer

visage.scene.galaxy_layer.GalaxyLayer

Manages the galaxy point-cloud actor(s) inside a PyVista Plotter.

Source code in visage/scene/galaxy_layer.py

class GalaxyLayer:
    """Manages the galaxy point-cloud actor(s) inside a PyVista Plotter."""

    def __init__(
        self,
        plotter: pv.Plotter,
        color_mode: ColorMode = "structure",
        colormap: str = "plasma",
        opacity: float = 1.0,
        visible: bool = True,
    ) -> None:
        self._pl = plotter
        self._color_mode: ColorMode = color_mode
        self._colormap = colormap
        self._opacity = opacity
        self._visible = visible
        self._actors: list = []
        # (actor, floor, multiplier) for every actor whose opacity is an
        # actor property, so the opacity slider can change it in place
        # instead of rebuilding the geometry.
        self._opacity_terms: list[tuple] = []
        # The photometry stack bakes opacity into per-point alpha, so it
        # can only follow the slider through a redraw.
        self._opacity_baked = False
        # Emissive (additive) splats are the default look: overlapping
        # galaxies add light instead of alpha-blending toward grey.  The
        # haloes stay alpha-blended — they are a scaffold, not a source.
        self._emissive = True
        # Multiplies how much light the emissive layers put out.
        self._glow_strength = _DEFAULT_GLOW_STRENGTH
        self._cloud: pv.PolyData | None = None  # persistent geometry
        self._render_params: tuple = ()  # tracks need-to-rebuild
        # Per-galaxy glow dimming for the galaxies being drawn (see
        # `glow_crowding`), and the inputs it was computed from.
        self._dim: np.ndarray | None = None
        self._dim_snap: GalaxySnapshot | None = None
        self._dim_key: tuple | None = None
        self._dim_cached: np.ndarray | None = None
        self._snapshot: GalaxySnapshot | None = None
        self._focus_mask: np.ndarray | None = None
        self._filter_mask: np.ndarray | None = None
        self._offset: np.ndarray = np.zeros(3, dtype=np.float32)
        # ── Lightcone extensions (both no-ops in normal box Explore mode) ──
        # World-radius multiplier: galaxies are Rvir-sized splats, sub-pixel
        # across a ~575 Mpc/h lightcone, so the lightcone layer scales them up.
        self._radius_scale: float = 1.0
        # Lightcone redshift/time cut: a keep-mask (True = shown).  The slider
        # sets it to remove part of the cone (e.g. a near-side redshift cut).
        # Composed (AND) with the focus/filter masks; None = no cut.
        self._slice_mask: np.ndarray | None = None
        # Always render the disk/bulge inner stellar layers (normally focus-only)
        # — the lightcone enables these so galaxies look more defined.
        self._show_inner_layers: bool = False
        # Synthetic-photometry stack (lightcone SED): the ordered list of
        # sed_mags keys the "sedstack" colour mode composites into a false-
        # colour image (each band tinted its representative colour, summed).
        self._sed_bands: list[str] = []

    # False-colour composite tuning (see _compute_sed_stack_rgb):
    #   _SED_GAMMA      — power applied to each band's median-balanced flux to
    #                     exaggerate the small real SED-shape ratios into
    #                     visible hue differences.
    #   _SED_SATURATION — how far to push the per-galaxy colour away from grey.
    #   _SED_FLOOR      — minimum brightness for any galaxy that has light, so
    #                     the faint end stays visible instead of being crushed
    #                     away (every galaxy with flux renders).
    _SED_GAMMA = 2.5
    _SED_SATURATION = 2.0
    _SED_FLOOR = 0.18

    # ------------------------------------------------------------------
    # Public API
    # ------------------------------------------------------------------

    @property
    def visible(self) -> bool:
        return self._visible

    @visible.setter
    def visible(self, value: bool) -> None:
        self._visible = value
        for actor in self._actors:
            actor.SetVisibility(value)
        self._pl.render()

    @property
    def glow_strength(self) -> float:
        return self._glow_strength

    @glow_strength.setter
    def glow_strength(self, value: float) -> None:
        value = max(0.05, float(value))
        if value == self._glow_strength:
            return
        self._glow_strength = value
        # Emissive brightness lives in the colours, so this needs a rebuild.
        if self._snapshot is not None and self._emissive:
            self._redraw()

    @property
    def emissive(self) -> bool:
        return self._emissive

    @emissive.setter
    def emissive(self, value: bool) -> None:
        value = bool(value)
        if value == self._emissive:
            return
        self._emissive = value
        if self._snapshot is not None:
            self._redraw()

    @property
    def opacity(self) -> float:
        return self._opacity

    @opacity.setter
    def opacity(self, value: float) -> None:
        self._opacity = float(value)
        if self._snapshot is None:
            return
        if self._opacity_baked or self._emissive or not self._opacity_terms:
            self._redraw()
            return
        # Cheap path: opacity is just an actor property, so there is no
        # reason to rebuild point clouds while the slider is dragged.
        for actor, floor, mul in self._opacity_terms:
            actor.prop.opacity = max(floor, self._opacity * mul)

    @property
    def color_mode(self) -> ColorMode:
        return self._color_mode

    @color_mode.setter
    def color_mode(self, value: ColorMode) -> None:
        self._color_mode = value
        if self._snapshot is not None:
            self._redraw()

    @property
    def colormap(self) -> str:
        return self._colormap

    @colormap.setter
    def colormap(self, value: str) -> None:
        self._colormap = value
        if self._snapshot is not None:
            self._redraw()

    def set_offset(self, offset: np.ndarray) -> None:
        self._offset = np.asarray(offset, dtype=np.float32)
        if self._snapshot is not None:
            self._redraw()

    def update(self, snapshot: GalaxySnapshot) -> None:
        self._snapshot = snapshot
        self._redraw()

    def set_radius_scale(self, scale: float) -> None:
        """Multiply per-galaxy world radii (lightcone visibility). 1.0 = box."""
        self._radius_scale = max(float(scale), 1e-6)
        if self._snapshot is not None:
            self._redraw()

    def set_slice_mask(self, mask: np.ndarray | None) -> None:
        """Lightcone redshift/time cut: keep only galaxies where mask is True
        (removes the rest from the render).  None disables the cut."""
        self._slice_mask = mask
        if self._snapshot is not None:
            self._redraw()

    def set_sed_bands(self, bands: list[str]) -> None:
        """Set the ordered sed_mags keys composited by the 'sedstack' colour
        mode (a false-colour stack of one or more filters)."""
        self._sed_bands = list(bands)
        if self._snapshot is not None:
            self._redraw()

    def set_show_inner_layers(self, on: bool) -> None:
        """Always draw the disk/bulge inner stellar layers (otherwise only
        shown inside a focus region).  Used by the lightcone for more definition.
        """
        self._show_inner_layers = bool(on)
        if self._snapshot is not None:
            self._redraw()

    def set_mask(self, mask: np.ndarray | None) -> None:
        """Backwards-compatible: sets the focus mask."""
        self.set_focus_mask(mask)

    def set_focus_mask(self, mask: np.ndarray | None) -> None:
        """Spatial focus mask (from sphere/box zoom). None = no focus."""
        self._focus_mask = mask
        if self._snapshot is not None:
            self._redraw()

    def set_filter_mask(self, mask: np.ndarray | None) -> None:
        """Property filter mask (from Filters tab). None = no filtering."""
        self._filter_mask = mask
        if self._snapshot is not None:
            self._redraw()

    def _combined_mask(self) -> np.ndarray | None:
        # AND together focus + filter + slice (lightcone cut) keep-masks.
        masks = [
            m
            for m in (self._focus_mask, self._filter_mask, self._slice_mask)
            if m is not None
        ]
        if not masks:
            return None
        n = len(masks[0])
        if any(len(m) != n for m in masks):
            # Masks are from different snapshots mid-transition; can't safely
            # combine them.  Return None so _redraw() shows everything until
            # they're refreshed for the new snapshot.
            return None
        out = masks[0].copy()
        for m in masks[1:]:
            out = out & m
        return out

    # ------------------------------------------------------------------
    # Internal
    # ------------------------------------------------------------------

    def _layer_weight(self, floor: float, mul: float) -> float:
        """How strongly a layer contributes, from the opacity setting."""
        weight = max(floor, self._opacity * mul)
        if self._emissive:
            weight *= self._glow_strength
        return weight

    def _layer_alpha(self, floor: float, mul: float) -> float:
        """Actor opacity for a layer.

        Emissive splats are scaled by the actor's opacity before being
        added, so a layer at 0.15-0.5 alpha comes out clamped and dim.  In
        that mode the actor stays fully opaque and the layer's weighting
        is carried by the colour instead (see `_rgba`).
        """
        if self._emissive:
            return 1.0
        return self._layer_weight(floor, mul)

    def _rgba(self, scalars: np.ndarray, cmap: str, weight=1.0) -> np.ndarray:
        """Splat colours, asinh-shaped when the splats are emissive.

        The curve goes on the scalar, before the colormap lookup — not on
        the colour that comes out of it.  A sequential colormap already
        carries brightness, so scaling its output would count the value
        twice and darken everything instead of shaping the roll-off.

        ``weight`` is one number for the layer or one per point.
        """
        v = np.clip(np.asarray(scalars, dtype=np.float32), 0.0, 1.0)
        if not self._emissive:
            return scalars_to_rgba(v, cmap)
        v = np.arcsinh(_GLOW_KNEE * v) / np.arcsinh(_GLOW_KNEE)
        rgba = scalars_to_rgba(v, cmap).copy()
        w = np.asarray(weight, dtype=np.float32).reshape(-1, 1)
        rgba[:, :3] = np.clip(rgba[:, :3] * w, 0.0, 255.0).astype(np.uint8)
        return rgba

    def _hue_rgba(
        self,
        scalars: np.ndarray,
        rgb: tuple[int, int, int],
        weight=1.0,
        floor: float = 0.3,
    ) -> np.ndarray:
        """Single-hue splat colours whose brightness rises with the scalar.

        Brightness runs from `floor` (scalar 0) to full (scalar 1), with the
        same asinh shaping and layer weighting as `_rgba` when emissive.
        ``weight`` is one number for the layer or one per point.
        """
        v = np.clip(np.asarray(scalars, dtype=np.float32), 0.0, 1.0)
        if self._emissive:
            v = np.arcsinh(_GLOW_KNEE * v) / np.arcsinh(_GLOW_KNEE)
        else:
            weight = 1.0
        bright = (floor + (1.0 - floor) * v) * np.asarray(
            weight, dtype=np.float32
        )
        rgba = np.full((len(v), 4), 255, dtype=np.uint8)
        rgba[:, :3] = np.clip(
            bright[:, None] * np.asarray(rgb, dtype=np.float32), 0.0, 255.0
        ).astype(np.uint8)
        return rgba

    def _clear_actors(self) -> None:
        for actor in self._actors:
            self._pl.remove_actor(actor, render=False)
        self._actors.clear()
        self._opacity_terms.clear()
        self._opacity_baked = False

    def _track_opacity(self, actor, floor: float, mul: float) -> None:
        """Remember how this actor's opacity follows the slider."""
        self._opacity_terms.append((actor, float(floor), float(mul)))

    def _redraw(self) -> None:
        snap = self._snapshot
        if snap is None or snap.count == 0:
            self._clear_actors()
            self._cloud = None
            return

        # Combined focus + filter + lightcone-slice keep-mask.  A lightcone
        # redshift/time cut removes galaxies here just like a property filter.
        mask = self._combined_mask()
        if mask is not None and len(mask) == snap.count:
            snap = self._subset(snap, mask)
            if snap.count == 0:
                self._clear_actors()
                self._cloud = None
                return

        radii = galaxy_world_radii_rvir(
            snap.rvir, snap.mvir, scale=self._radius_scale
        )
        eff_pos = snap.positions + self._offset

        # Every mode shares the same Structure composition (BH core, cold-gas
        # envelope, stellar particles, CGM/Hot outer envelope).  When the
        # mode isn't 'structure', we add ONE more outermost layer whose
        # colour comes from the active Colour-by + the galaxy colormap.
        self._clear_actors()
        self._cloud = None
        self._render_params = ()
        self._render_composition(snap, radii, eff_pos)

    def _render_composition(
        self,
        snap: GalaxySnapshot,
        radii: np.ndarray,
        eff_pos: np.ndarray,
    ) -> None:
        """The complete Explore galaxy render: the structure composition plus,
        for non-structure Colour-by modes, the outer property halo.  Shared by
        box Explore and the lightcone (dim full cone + bright shell) so both
        paths draw identical splats, layers and colormaps."""
        if self._color_mode == "sedstack":
            # Photometry is its own thing: a SINGLE Rvir-sized splat layer
            # painted with the false-colour filter/M-L stack — NOT the nested
            # structure composition. This is the "build an image from the
            # filters" view, so it deliberately skips the CGM/Hot/cold/BH
            # inner layers the other modes draw.
            rgb = self._compute_sed_stack_rgb(snap)
            if rgb is not None:
                self._render_sed_stack(eff_pos, rgb, radii)
            return
        self._dim = self._crowd_dimming(snap, radii)
        self._render_structure(snap, radii, eff_pos)
        if self._color_mode == "type":
            mass_colors = normalize_log(
                snap.stellar_mass, *_RANGES["stellar_mass"]
            )
            for tmask, cmap in [
                (snap.gal_type == 0, _CENTRAL_CMAP),
                (snap.gal_type > 0, _SATELLITE_CMAP),
            ]:
                if not np.any(tmask):
                    continue
                self._render_outer_property(
                    eff_pos[tmask],
                    mass_colors[tmask],
                    radii[tmask],
                    cmap,
                    self._dim_of(tmask),
                )
        elif self._color_mode != "structure":
            colors = self._compute_colors(snap)
            self._render_outer_property(
                eff_pos, colors, radii, self._colormap, self._dim
            )

    def _crowd_dimming(
        self, snap: GalaxySnapshot, radii: np.ndarray
    ) -> np.ndarray | None:
        """Per-galaxy brightness factor that keeps overlapping glows from
        clipping to white (None when the splats are not emissive).

        The loader works out the crowding of every galaxy in a snapshot
        while reading it (a lightcone's is worked out here, once), so
        changing snapshot costs nothing here.  A focus region or lightcone
        cut removes whole regions, which leaves the crowding of the
        galaxies kept much as it was, so it is reused.  A property filter
        that leaves well under the whole snapshot thins out clusters
        instead, so the crowding is redone for the galaxies drawn and what
        remains brightens; that costs time in proportion to the galaxies
        drawn.  None of this depends on the camera or colours, so
        restyling reuses it.
        """
        if not self._emissive or snap.count == 0:
            return None
        full = self._snapshot
        mask = self._combined_mask()
        if mask is not None and mask.size != full.count:
            mask = None
        thinned = (
            mask is not None
            and self._filter_mask is not None
            and int(np.count_nonzero(self._filter_mask)) < 0.7 * full.count
        )
        reuse = self._radius_scale == 1.0 and not thinned
        key = (
            None if mask is None else hash(np.packbits(mask).tobytes()),
            self._radius_scale,
            reuse,
        )
        if full is self._dim_snap and key == self._dim_key:
            return self._dim_cached
        if not reuse:
            n = glow_crowding(snap.positions, radii)
        else:
            if full.crowding is None:
                full.crowding = glow_crowding(
                    full.positions,
                    galaxy_world_radii_rvir(full.rvir, full.mvir),
                )
            n = full.crowding if mask is None else full.crowding[mask]
        dim = np.minimum(1.0, (n / _CROWD_FREE) ** -_CROWD_POWER)
        self._dim_cached = dim.astype(np.float32)
        self._dim_snap, self._dim_key = full, key
        return self._dim_cached

    def _shape_splat(self, mapper) -> None:
        """Give an emissive layer's splats the core-and-glow profile."""
        if self._emissive:
            mapper.SetSplatShaderCode(_SPLAT_SHADER)

    def _dim_of(self, mask: np.ndarray) -> np.ndarray | float:
        """The crowding factor for a subset of the drawn galaxies."""
        return 1.0 if self._dim is None else self._dim[mask]

    @staticmethod
    def _subset(snap: GalaxySnapshot, mask: np.ndarray) -> GalaxySnapshot:
        """Field-agnostic slice of a GalaxySnapshot by a boolean mask."""
        from dataclasses import fields as _dc_fields

        from visage.io.galaxy_reader import GalaxySnapshot as _GS

        n = snap.count
        kwargs = {}
        for fld in _dc_fields(snap):
            v = getattr(snap, fld.name)
            if isinstance(v, np.ndarray) and len(v) == n:
                v = v[mask]
            elif isinstance(v, dict):
                # e.g. sed_mags: dict[str, np.ndarray] — slice every array
                # inside it too, so it stays aligned with the sliced snapshot.
                v = {
                    k: (
                        arr[mask]
                        if isinstance(arr, np.ndarray) and len(arr) == n
                        else arr
                    )
                    for k, arr in v.items()
                }
            kwargs[fld.name] = v
        return _GS(**kwargs)

    def _render_by_type(self, snap: GalaxySnapshot, radii: np.ndarray) -> None:
        mass_colors = normalize_log(
            snap.stellar_mass, *_RANGES["stellar_mass"]
        )
        for mask, cmap in [
            (snap.gal_type == 0, _CENTRAL_CMAP),
            (snap.gal_type > 0, _SATELLITE_CMAP),
        ]:
            if not np.any(mask):
                continue
            self._render_gaussian(
                snap.positions[mask], mass_colors[mask], radii[mask], cmap
            )

    def _render_structure(
        self,
        snap: GalaxySnapshot,
        radii: np.ndarray,
        positions: np.ndarray | None = None,
    ) -> None:
        """Multi-layer physically-suggestive galaxy rendering.

        Always rendered (all galaxies):
          • blue cold-gas envelope sized by ColdGas — the base every galaxy
            gets
          • a green CGM shell and a red hot-atmosphere shell, each drawn
            for the galaxies that actually carry that mass.  A galaxy with
            both gets both, at their own radii; a galaxy with one gets one.

        Only when a focus region is active:
          • pale-blue disk layer (StellarMass − BulgeMass)
          • amber bulge layer (BulgeMass)

        Gas shells and the bulge are single-hue: mass sets brightness and
        size, so a massive component is always the brighter one.

        All layers share the per-galaxy world-space `radii` envelope so the
        overall splat size stays consistent with the standard rendering.
        """
        if snap.count == 0:
            return

        pos = snap.positions if positions is None else positions
        # ---- Per-galaxy radii (Mpc/h) keyed off the subhalo Rvir --------
        r_outer = np.maximum(radii, 1e-4)
        r_cold = 0.45 * r_outer

        # Convenience clamped log10
        def _logn(x, vmin, vmax):
            log = np.log10(np.maximum(x, 1.0))
            return np.clip((log - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0)

        cold_scalar = _logn(snap.cold_gas, *_COLD_RANGE)
        # Crowding factor per galaxy (1.0 when the splats are not emissive).
        dim = 1.0 if self._dim is None else self._dim

        # ---- (1) Halo-gas shells --------------------------------------
        # CGM and hot atmosphere are separate components, not two sides of
        # a switch: a galaxy can carry both, and the regime flag only says
        # which one dominates.  Each shell is drawn for the galaxies whose
        # mass in that component is non-zero, so a mixed galaxy shows a
        # green CGM shell inside a red hot one and a single-component
        # galaxy shows just its own.  (Before, the regime flag picked one
        # and the other was never drawn — and a model with no Regime field
        # rendered every galaxy as hot.)
        shells = (
            ("cgm", snap.cgm_gas, _CGM_RGB, _CGM_RANGE, 0.85),
            ("hot", snap.hot_gas, _HOT_RGB, _HOT_RANGE, 1.00),
        )
        for _name, mass, rgb, (lo, hi), r_scale in shells:
            if mass.size != snap.count:
                continue
            mask = mass > 0.0
            if not np.any(mask):
                continue
            scalar = _logn(mass[mask], lo, hi)
            cloud = pv.PolyData(pos[mask])
            cloud["rgba"] = self._hue_rgba(
                scalar, rgb, self._layer_weight(0.2, 0.5) * self._dim_of(mask)
            )
            cloud["radius"] = (
                r_outer[mask] * r_scale * (0.5 + 0.5 * scalar)
            ).astype(np.float32)
            actor = self._pl.add_mesh(
                cloud,
                scalars="rgba",
                rgb=True,
                style="points_gaussian",
                emissive=self._emissive,
                opacity=self._layer_alpha(0.2, 0.5),
                show_scalar_bar=False,
                render=False,
                reset_camera=False,
            )
            mp = actor.mapper
            mp.SetScaleArray("radius")
            mp.SetScaleFactor(1.0)
            self._shape_splat(mp)
            if not self._visible:
                actor.SetVisibility(False)
            self._track_opacity(actor, 0.2, 0.5)
            self._actors.append(actor)

        # ---- (2) Cold-gas blue envelope -------------------------------
        cloud = pv.PolyData(pos)
        cloud["rgba"] = self._hue_rgba(
            cold_scalar, _COLD_RGB, self._layer_weight(0.2, 0.5) * dim
        )
        cloud["radius"] = (r_cold * (0.5 + 0.5 * cold_scalar)).astype(
            np.float32
        )
        actor = self._pl.add_mesh(
            cloud,
            scalars="rgba",
            rgb=True,
            style="points_gaussian",
            emissive=self._emissive,
            opacity=self._layer_alpha(0.2, 0.5),
            show_scalar_bar=False,
            render=False,
            reset_camera=False,
        )
        mp = actor.mapper
        mp.SetScaleArray("radius")
        mp.SetScaleFactor(1.0)
        self._shape_splat(mp)
        if not self._visible:
            actor.SetVisibility(False)
        self._track_opacity(actor, 0.2, 0.5)
        self._actors.append(actor)

        # (Per-galaxy star scatter and BH accretion-disk cores both
        # removed — invisible / negligible at typical zoom levels and
        # together they were the bulk of the per-frame splat cost.)

        # ---- (3) Inner stellar layers (disk/bulge) -------------------
        # Normally focus-only (a focus region active); at full-scene box scale
        # they'd be invisible.  The lightcone opts in via _show_inner_layers so
        # its galaxies get the extra disk/bulge definition everywhere.
        if self._focus_mask is not None or self._show_inner_layers:
            disk_mass = np.maximum(snap.stellar_mass - snap.bulge_mass, 0.0)
            disk_scalar = _logn(disk_mass, 7.0, 12.0)
            bulge_scalar = _logn(snap.bulge_mass, 6.0, 12.0)

            r_disk = (r_outer * 0.22 * (0.35 + 0.65 * disk_scalar)).astype(
                np.float32
            )
            r_bulge = (r_outer * 0.12 * (0.30 + 0.70 * bulge_scalar)).astype(
                np.float32
            )

            weight = self._layer_weight(0.5, 0.75) * dim
            for rgba, radii_arr in [
                (self._rgba(disk_scalar, "Blues_r", weight), r_disk),
                (self._hue_rgba(bulge_scalar, _BULGE_RGB, weight), r_bulge),
            ]:
                cloud = pv.PolyData(pos)
                cloud["rgba"] = rgba
                cloud["radius"] = radii_arr
                actor = self._pl.add_mesh(
                    cloud,
                    scalars="rgba",
                    rgb=True,
                    style="points_gaussian",
                    emissive=self._emissive,
                    opacity=self._layer_alpha(0.5, 0.75),
                    show_scalar_bar=False,
                    render=False,
                    reset_camera=False,
                )
                mp = actor.mapper
                mp.SetScaleArray("radius")
                mp.SetScaleFactor(1.0)
                self._shape_splat(mp)
                if not self._visible:
                    actor.SetVisibility(False)
                self._track_opacity(actor, 0.5, 0.75)
                self._actors.append(actor)

    def _render_outer_property(
        self,
        positions: np.ndarray,
        colors: np.ndarray,
        radii: np.ndarray,
        cmap: str,
        dim=1.0,
    ) -> None:
        """Outermost halo around the Structure composition, coloured by the
        active galaxy Colour-by mode + chosen colormap.  Slightly larger and
        more transparent than the CGM/Hot envelope so the inner Structure
        layers stay visible."""
        if len(positions) == 0:
            return
        cloud = pv.PolyData(positions)
        cloud["rgba"] = self._rgba(
            colors, cmap, self._layer_weight(0.12, 0.25) * dim
        )
        # Sit ~30% beyond the standard envelope.  This is the "Colour-by" halo.
        cloud["radius"] = (radii * 1.3).astype(np.float32)
        actor = self._pl.add_mesh(
            cloud,
            scalars="rgba",
            rgb=True,
            style="points_gaussian",
            emissive=self._emissive,
            # Subtle so the inner Structure detail isn't drowned
            opacity=self._layer_alpha(0.12, 0.25),
            show_scalar_bar=False,
            render=False,
            reset_camera=False,
        )
        mp = actor.mapper
        mp.SetScaleArray("radius")
        mp.SetScaleFactor(1.0)
        self._shape_splat(mp)
        if not self._visible:
            actor.SetVisibility(False)
        self._track_opacity(actor, 0.12, 0.25)
        self._actors.append(actor)

    def _render_gaussian(
        self,
        positions: np.ndarray,
        colors: np.ndarray,
        radii: np.ndarray,
        cmap: str,
    ) -> None:
        if len(positions) == 0:
            return
        cloud = pv.PolyData(positions)
        cloud["rgba"] = self._rgba(colors, cmap)
        cloud["radius"] = radii
        actor = self._pl.add_mesh(
            cloud,
            scalars="rgba",
            rgb=True,
            style="points_gaussian",
            emissive=self._emissive,
            opacity=self._layer_alpha(0.0, 1.0),
            show_scalar_bar=False,
            render=False,
            reset_camera=False,
        )
        # Make the gaussian splats sized in world coordinates (Mpc/h) via
        # the per-point "radius" array rather than fixed screen pixels.
        mapper = actor.mapper
        mapper.SetScaleArray("radius")
        mapper.SetScaleFactor(1.0)
        self._shape_splat(mapper)
        if not self._visible:
            actor.SetVisibility(False)
        self._cloud = cloud
        self._track_opacity(actor, 0.0, 1.0)
        self._actors.append(actor)

    def _compute_colors(self, snap: GalaxySnapshot) -> np.ndarray:
        # Photometry (sedstack) has its own RGB path and never reaches here.
        m = self._color_mode
        if m == "ssfr":
            return normalize_log(snap.ssfr, *_RANGES["ssfr"])
        if m == "sfr":
            return normalize_log(np.maximum(snap.sfr, 1e-6), *_RANGES["sfr"])
        if m == "cold_gas":
            return normalize_log(
                np.maximum(snap.cold_gas, 1.0), *_RANGES["cold_gas"]
            )
        if m == "bulge_mass":
            return normalize_log(
                np.maximum(snap.bulge_mass, 1.0), *_RANGES["bulge_mass"]
            )
        if m == "bh_mass":
            return normalize_log(
                np.maximum(snap.bh_mass, 1.0), *_RANGES["bh_mass"]
            )
        if m == "ics_mass":
            return normalize_log(
                np.maximum(snap.ics_mass, 1.0), *_RANGES["ics_mass"]
            )
        if m == "bt":
            bt = snap.bulge_mass / np.where(
                snap.stellar_mass > 0, snap.stellar_mass, np.inf
            )
            vmin, vmax = _RANGES["bt"]
            return np.clip(
                (bt - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0
            ).astype(np.float32)
        if m == "age":
            ages = snap.mean_age.astype(np.float32)
            vmin, vmax = _RANGES["age"]
            return np.clip((ages - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0)
        if m.startswith("extra:"):
            return self._extra_colors(snap, m[len("extra:") :])
        if m in _LOG_FIELDS:
            attr, floor = _LOG_FIELDS[m]
            return normalize_log(
                np.maximum(getattr(snap, attr), floor), *_RANGES[m]
            )
        return normalize_log(snap.stellar_mass, *_RANGES["stellar_mass"])

    def _extra_colors(self, snap: GalaxySnapshot, key: str) -> np.ndarray:
        """Normalised colours for an auto-discovered property.

        The field's meaning (and so its natural range) isn't known, so the
        scale comes from the data: log10 when it is non-negative and spans
        several decades, linear otherwise, clipped to the 1st–99th
        percentile so a handful of outliers can't flatten everything else.
        """
        v = snap.extra.get(key)
        if v is None or len(v) != snap.count:
            return np.zeros(snap.count, dtype=np.float32)
        v = np.asarray(v, dtype=np.float64)
        finite = np.isfinite(v)
        if not np.any(finite):
            return np.zeros(snap.count, dtype=np.float32)
        pos = v[finite & (v > 0)]
        if (
            v[finite].min() >= 0.0
            and pos.size > 0
            and pos.max() / max(pos.min(), 1e-300) >= 1.0e3
        ):
            v = np.log10(np.maximum(v, pos.min()))
            finite = np.isfinite(v)
        lo, hi = np.percentile(v[finite], [1.0, 99.0])
        if hi <= lo:
            hi = lo + 1e-10
        out = np.clip((v - lo) / (hi - lo), 0.0, 1.0)
        return np.where(finite, out, 0.0).astype(np.float32)

    def _sed_stack_amplitude(
        self, snap: GalaxySnapshot, item: str
    ) -> tuple[np.ndarray, str]:
        """(linear amplitude, filter_name) for one stack item.

        Filter band  -> relative flux f = 10**(-0.4*mag).
        'ml:' item   -> mass-to-light ratio M*/L in that band.

        The amplitude is then divided by its own median over the galaxies
        that have it, so the *typical* galaxy sits at ~1 in every channel
        (i.e. neutral/white) and it's the DEVIATIONS from typical — a galaxy
        brighter than average in r, fainter in u — that show up as colour.
        This is what makes combining filters produce reds / blues / yellows
        instead of a flat white: the composite keeps the relative amplitudes
        between bands (the real SED shape) rather than normalising each band
        independently to the same range."""
        key = item[len("ml:") :] if item.startswith("ml:") else item
        filt = key
        for pre in ("mag_rest_", "mag_obs_"):
            if filt.startswith(pre):
                filt = filt[len(pre) :]
                break

        mags = snap.sed_mags.get(key)
        if mags is None or len(mags) != snap.count:
            return np.zeros(snap.count), filt
        mags = np.asarray(mags, dtype=np.float64)
        finite = np.isfinite(mags)

        with np.errstate(over="ignore", invalid="ignore"):
            if item.startswith("ml:"):
                from visage.sed.filters import SOLAR_ABSMAG_AB

                m_sun = SOLAR_ABSMAG_AB.get(filt)
                if m_sun is None:
                    return np.zeros(snap.count), filt
                lum = 10.0 ** (-0.4 * (mags - m_sun))  # L / Lsun
                mass = np.maximum(snap.stellar_mass.astype(np.float64), 0.0)
                amp = np.where(lum > 0, mass / lum, 0.0)
            else:
                amp = 10.0 ** (-0.4 * mags)  # relative flux

        amp = np.where(finite & np.isfinite(amp), amp, 0.0)
        good = amp[amp > 0]
        if good.size:
            med = float(np.median(good))
            if med > 0:
                amp = amp / med  # typical galaxy -> ~1 in this channel
        return amp, filt

    def _compute_sed_stack_rgb(self, snap: GalaxySnapshot):
        """Per-galaxy RGB (N, 3) in [0, 1] for the 'sedstack' mode — a proper
        astronomical false-colour composite (à la Lupton et al. 2004).

        Hue and brightness are computed separately so the two goals don't
        fight each other:

        * HUE — each item's median-balanced flux/M*L is raised to a power
          (_SED_GAMMA) and tinted its representative colour; the per-galaxy
          channel ratios then carry the (contrast-boosted, saturated) SED
          shape, so galaxies read visibly blue / red / yellow.
        * BRIGHTNESS — an asinh stretch of the RAW (un-boosted) luminance with
          a floor, so EVERY galaxy that has light is visible. Brightness must
          not use the gamma-boosted amplitude or the faint end gets crushed to
          nothing and ~a quarter of the galaxies vanish.

        Returns None if there's nothing to show."""
        from visage.sed.filters import band_colour

        n = snap.count
        items = [
            it
            for it in self._sed_bands
            if it.replace("ml:", "") in snap.sed_mags
        ]
        if n == 0 or not items:
            return None

        rgb_hue = np.zeros((n, 3), dtype=np.float64)  # gamma-boosted -> hue
        lum = np.zeros(n, dtype=np.float64)  # raw flux sum -> brightness
        for item in items:
            amp, filt = self._sed_stack_amplitude(snap, item)
            cr, cg, cb = band_colour(filt)
            boosted = np.where(amp > 0, amp**self._SED_GAMMA, 0.0)
            rgb_hue[:, 0] += boosted * cr
            rgb_hue[:, 1] += boosted * cg
            rgb_hue[:, 2] += boosted * cb
            lum += amp

        lit = lum > 0
        if not lit.any():
            return np.zeros((n, 3), dtype=np.float32)

        # Hue: unit-ish colour, pushed away from grey so it reads vividly.
        hmax = rgb_hue.max(axis=1)
        hue = rgb_hue / np.maximum(hmax[:, None], 1e-12)
        gray = hue.mean(axis=1, keepdims=True)
        hue = np.clip(gray + self._SED_SATURATION * (hue - gray), 0.0, 1.0)

        # Brightness: asinh stretch of the raw luminance, softened by the
        # median so most galaxies sit in the mid-tones, then lifted off the
        # floor so even the faintest galaxies with light stay visible.
        soft = float(np.median(lum[lit])) or 1.0
        stretched = np.arcsinh(lum / soft)
        norm = float(np.percentile(stretched[lit], 99.0)) or 1.0
        s_raw = np.clip(stretched / norm, 0.0, 1.0)
        s = np.where(
            lit, self._SED_FLOOR + (1.0 - self._SED_FLOOR) * s_raw, 0.0
        )

        out = hue * s[:, None]
        return np.clip(out, 0.0, 1.0).astype(np.float32)

    # Nested gaussian shells for the photometry splat: (radius scale, alpha
    # fraction). A sharply peaked profile — a small, near-opaque bright core
    # inside a few smaller, progressively fainter haloes — gives each galaxy a
    # crisp, defined point rather than a soft blob. The outer halo is kept
    # well under the full Rvir so neighbouring splats stay distinct. The alpha
    # fraction scales the per-galaxy brightness (see _render_sed_stack).
    _SED_SHELLS = (
        (0.60, 0.12),
        (0.36, 0.25),
        (0.20, 0.45),
        (0.10, 0.75),
        (0.05, 1.0),
    )

    def _render_sed_stack(
        self,
        positions: np.ndarray,
        rgb: np.ndarray,
        radii: np.ndarray,
    ) -> None:
        """The whole photometry view: an Rvir-sized gaussian splat per galaxy
        painted with per-point RGBA directly (no colormap, no inner structure
        layers) — a clean false-colour image built from the filter stack.

        Each splat's ALPHA is its composite brightness (times the shell
        weight), so a faint galaxy is *transparent* rather than an opaque dark
        disk. That's what keeps dim galaxies in the foreground from painting
        dark spots over brighter ones behind them, while a bright galaxy still
        gets an opaque, defined core. The nested shells (small bright core +
        fainter haloes) give the splats their peaked, defined profile."""
        if len(positions) == 0:
            return
        rgb = np.clip(rgb, 0.0, 1.0)
        brightness = rgb.max(axis=1)
        keep = brightness > 0.01  # skip the truly-invisible; dim ones stay
        if not keep.any():
            return
        pos = positions[keep]
        rgb255 = (rgb[keep] * 255).astype(np.uint8)
        base_a = brightness[keep]  # 0..1 per-galaxy opacity
        rad = radii[keep]
        first = True
        for rscale, afrac in self._SED_SHELLS:
            cloud = pv.PolyData(pos)
            rgba = np.empty((len(pos), 4), dtype=np.uint8)
            rgba[:, :3] = rgb255
            rgba[:, 3] = np.clip(
                base_a * afrac * self._opacity * 255.0, 0, 255
            ).astype(np.uint8)
            # VTK point-gaussian reads 0..255 uint8 RGBA as direct scalars.
            cloud["sed_rgba"] = rgba
            cloud["radius"] = (rad * rscale).astype(np.float32)
            actor = self._pl.add_mesh(
                cloud,
                scalars="sed_rgba",
                rgb=True,  # direct scalars; 4 components => RGBA
                style="points_gaussian",
                # Never emissive: the photometry composite is already a
                # calibrated brightness (median-balanced flux, asinh
                # stretch) carried in per-point alpha.  Adding it instead
                # of blending it would throw that calibration away.
                emissive=False,
                opacity=1.0,  # per-point alpha carries the opacity
                show_scalar_bar=False,
                render=False,
                reset_camera=False,
            )
            mp = actor.mapper
            mp.SetScaleArray("radius")
            mp.SetScaleFactor(1.0)
            if not self._visible:
                actor.SetVisibility(False)
            if first:
                self._cloud = cloud
                first = False
            self._opacity_baked = True
            self._actors.append(actor)

set_filter_mask(mask)

Property filter mask (from Filters tab). None = no filtering.

Source code in visage/scene/galaxy_layer.py

def set_filter_mask(self, mask: np.ndarray | None) -> None:
    """Property filter mask (from Filters tab). None = no filtering."""
    self._filter_mask = mask
    if self._snapshot is not None:
        self._redraw()

set_focus_mask(mask)

Spatial focus mask (from sphere/box zoom). None = no focus.

Source code in visage/scene/galaxy_layer.py

def set_focus_mask(self, mask: np.ndarray | None) -> None:
    """Spatial focus mask (from sphere/box zoom). None = no focus."""
    self._focus_mask = mask
    if self._snapshot is not None:
        self._redraw()

set_mask(mask)

Backwards-compatible: sets the focus mask.

Source code in visage/scene/galaxy_layer.py

def set_mask(self, mask: np.ndarray | None) -> None:
    """Backwards-compatible: sets the focus mask."""
    self.set_focus_mask(mask)

set_radius_scale(scale)

Multiply per-galaxy world radii (lightcone visibility). 1.0 = box.

Source code in visage/scene/galaxy_layer.py

def set_radius_scale(self, scale: float) -> None:
    """Multiply per-galaxy world radii (lightcone visibility). 1.0 = box."""
    self._radius_scale = max(float(scale), 1e-6)
    if self._snapshot is not None:
        self._redraw()

set_sed_bands(bands)

Set the ordered sed_mags keys composited by the 'sedstack' colour mode (a false-colour stack of one or more filters).

Source code in visage/scene/galaxy_layer.py

def set_sed_bands(self, bands: list[str]) -> None:
    """Set the ordered sed_mags keys composited by the 'sedstack' colour
    mode (a false-colour stack of one or more filters)."""
    self._sed_bands = list(bands)
    if self._snapshot is not None:
        self._redraw()

set_show_inner_layers(on)

Always draw the disk/bulge inner stellar layers (otherwise only shown inside a focus region). Used by the lightcone for more definition.

Source code in visage/scene/galaxy_layer.py

def set_show_inner_layers(self, on: bool) -> None:
    """Always draw the disk/bulge inner stellar layers (otherwise only
    shown inside a focus region).  Used by the lightcone for more definition.
    """
    self._show_inner_layers = bool(on)
    if self._snapshot is not None:
        self._redraw()

set_slice_mask(mask)

Lightcone redshift/time cut: keep only galaxies where mask is True (removes the rest from the render). None disables the cut.

Source code in visage/scene/galaxy_layer.py

def set_slice_mask(self, mask: np.ndarray | None) -> None:
    """Lightcone redshift/time cut: keep only galaxies where mask is True
    (removes the rest from the render).  None disables the cut."""
    self._slice_mask = mask
    if self._snapshot is not None:
        self._redraw()

Parallel loader

visage.parallel.loader.SnapshotLoader

Prefetch and cache haloes + galaxies around the current snapshot.

Uses a thread pool so the main render thread never blocks on disk I/O. Keeps an LRU in-memory cache of recently loaded snapshots.

Parameters:

Name	Type	Description	Default
config	SimConfig		required
snap_table	SnapshotTable		required
n_jobs	int | None	the CPUs this process may use, minus one)	None
prefetch_radius	int		2
cache_size	int		8
min_halo_mass	float		10000000000.0
min_stellar_mass	float		100000000.0
cache_bytes	int | None	physical RAM, or $VISAGE_CACHE_GB)	None
max_halos	int | None		None
max_galaxies	int | None		None

Source code in visage/parallel/loader.py

class SnapshotLoader:
    """Prefetch and cache haloes + galaxies around the current snapshot.

    Uses a thread pool so the main render thread never blocks on disk I/O.
    Keeps an LRU in-memory cache of recently loaded snapshots.

    Parameters
    ----------
    config:          parsed SimConfig
    snap_table:      SnapshotTable for the simulation
    n_jobs:          worker threads for parallel halo file reads (default:
                     the CPUs this process may use, minus one)
    prefetch_radius: number of snapshots ahead/behind to prefetch
    cache_size:      max snapshots kept in memory
    min_halo_mass:   Msun halo mass floor
    min_stellar_mass: Msun stellar mass floor
    cache_bytes:     memory budget for cached snapshots (None = half of
                     physical RAM, or $VISAGE_CACHE_GB)
    max_halos:       optional downsample ceiling (None = every halo)
    max_galaxies:    optional downsample ceiling (None = every galaxy)
    """

    def __init__(
        self,
        config: SimConfig,
        snap_table: SnapshotTable,
        n_jobs: int | None = None,
        prefetch_radius: int = 2,
        cache_size: int = 8,
        min_halo_mass: float = 1.0e10,
        min_stellar_mass: float = 1.0e8,
        max_halos: int | None = None,
        max_galaxies: int | None = None,
        cache_bytes: int | None = None,
        min_particles: int | None = None,
    ) -> None:
        self._cfg = config
        self._snap_table = snap_table
        self._n_jobs = (
            n_jobs if n_jobs is not None else system.default_workers()
        )
        self._prefetch_radius = prefetch_radius
        self._min_halo_mass = min_halo_mass
        self._min_stellar_mass = min_stellar_mass
        self._max_halos = max_halos
        self._max_galaxies = max_galaxies
        # Resolution (visage.utils.resolution): haloes of at least this
        # many particles, galaxies from the stellar mass that resolves
        # them in each snapshot; None keeps the fixed mass cuts.
        self._min_particles = min_particles
        if min_particles:
            self._min_halo_mass = (
                min_particles * config.part_mass * 1.0e10 / config.hubble_h
            )

        self._executor = ThreadPoolExecutor(
            max_workers=max(2, prefetch_radius * 2)
        )
        self._futures: dict[int, Future] = {}
        self._lock = Lock()
        self._tree_cache: dict[int, KDTree] = {}

        # Repeated requests for a snapshot skip disk entirely.  The cache
        # is bounded by memory rather than by a snapshot count, so a big
        # box fills it and then rolls, instead of holding every snapshot
        # of the run resident until the machine runs out of RAM.
        self._cache_bytes = (
            int(cache_bytes)
            if cache_bytes is not None
            else default_cache_bytes()
        )
        self._cached_load = _SnapshotCache(
            self._cache_bytes, self._load, on_evict=self._on_evict
        )

    def _load(self, snap_num: int) -> tuple[HaloSnapshot, GalaxySnapshot]:
        halos = load_halo_snapshot(
            tree_dir=self._cfg.tree_dir,
            tree_name=self._cfg.tree_name,
            tree_type=getattr(self._cfg, "tree_type", ""),
            snap_num=snap_num,
            first_file=self._cfg.first_file,
            last_file=self._cfg.last_file,
            mass_cut=self._min_halo_mass,
            max_halos=self._max_halos,
            hubble_h=self._cfg.hubble_h,
            n_jobs=self._n_jobs,
            box_size=self._cfg.box_size,
        )
        galaxies = load_galaxy_snapshot(
            hdf5_path=self._cfg.hdf5_path,
            snap_num=snap_num,
            min_stellar_mass=self._min_stellar_mass,
            max_galaxies=self._max_galaxies,
            hubble_h=self._cfg.hubble_h,
            scale_factors=self._snap_table.scale_factors,
            omega_m=self._cfg.omega,
            omega_l=self._cfg.omega_lambda,
            min_particles=self._min_particles,
            particle_mass=self._cfg.part_mass,
        )
        # Build the spatial index while still on the background thread so
        # snap navigation never blocks on KDTree construction (~50 ms / snap).
        if len(halos.positions) > 0:
            self._tree_cache[snap_num] = KDTree(halos.positions)
        # Likewise the glow crowding the galaxy layer dims by (~150 ms for
        # 100k galaxies).
        galaxies.crowding = glow_crowding(
            galaxies.positions,
            galaxy_world_radii_rvir(galaxies.rvir, galaxies.mvir),
        )
        return halos, galaxies

    def _on_evict(self, snap_num: int) -> None:
        """Drop everything else tied to an evicted snapshot."""
        self._tree_cache.pop(snap_num, None)

    @property
    def cache_nbytes(self) -> int:
        """Bytes currently held by cached snapshots."""
        return self._cached_load.nbytes

    def get_tree(self, snap_num: int) -> KDTree | None:
        """Return the pre-built KDTree for snap_num, or None if not ready."""
        return self._tree_cache.get(snap_num)

    def get(self, snap_num: int) -> tuple[HaloSnapshot, GalaxySnapshot]:
        """Return (HaloSnapshot, GalaxySnapshot) for snap_num.

        Blocks only on a cold-cache miss; otherwise returns from memory.
        Triggers background prefetch of neighbouring snapshots as a side-effect.
        """
        result = self._cached_load(snap_num)
        self._prefetch_neighbours(snap_num)
        return result

    def preload_all(self) -> list[Future]:
        """Kick off background loads of every snapshot. Returns the futures
        so a caller can track progress. Already-loaded / in-flight snapshots
        are not resubmitted."""
        # Silence per-snapshot load chatter so it doesn't bury the startup
        # browser URL in the terminal.
        from visage.io import halo_reader, galaxy_reader

        halo_reader.VERBOSE = False
        galaxy_reader.VERBOSE = False
        n = self._snap_table.count
        futures: list[Future] = []
        with self._lock:
            for snap in range(n):
                if snap not in self._futures:
                    self._futures[snap] = self._submit(snap)
                futures.append(self._futures[snap])
        return futures

    def _submit(self, snap_num: int) -> Future:
        """Background load whose Future is released once it completes.

        A Future keeps its result alive, so holding finished ones would pin
        every snapshot in memory and make the cache budget meaningless.
        """
        future = self._executor.submit(self._background_load, snap_num)
        future.add_done_callback(
            lambda _f, s=snap_num: self._futures.pop(s, None)
        )
        return future

    def _background_load(self, snap_num: int):
        _wait_until_quiet()
        return self._cached_load(snap_num)

    def _prefetch_neighbours(self, current: int) -> None:
        n = self._snap_table.count
        with self._lock:
            for offset in range(1, self._prefetch_radius + 1):
                for snap in (current - offset, current + offset):
                    if 0 <= snap < n and snap not in self._futures:
                        self._futures[snap] = self._submit(snap)

    def shutdown(self) -> None:
        self._executor.shutdown(wait=False, cancel_futures=True)

cache_nbytes property

Bytes currently held by cached snapshots.

get(snap_num)

Return (HaloSnapshot, GalaxySnapshot) for snap_num.

Blocks only on a cold-cache miss; otherwise returns from memory. Triggers background prefetch of neighbouring snapshots as a side-effect.

Source code in visage/parallel/loader.py

def get(self, snap_num: int) -> tuple[HaloSnapshot, GalaxySnapshot]:
    """Return (HaloSnapshot, GalaxySnapshot) for snap_num.

    Blocks only on a cold-cache miss; otherwise returns from memory.
    Triggers background prefetch of neighbouring snapshots as a side-effect.
    """
    result = self._cached_load(snap_num)
    self._prefetch_neighbours(snap_num)
    return result

get_tree(snap_num)

Return the pre-built KDTree for snap_num, or None if not ready.

Source code in visage/parallel/loader.py

def get_tree(self, snap_num: int) -> KDTree | None:
    """Return the pre-built KDTree for snap_num, or None if not ready."""
    return self._tree_cache.get(snap_num)

preload_all()

Kick off background loads of every snapshot. Returns the futures so a caller can track progress. Already-loaded / in-flight snapshots are not resubmitted.

Source code in visage/parallel/loader.py

def preload_all(self) -> list[Future]:
    """Kick off background loads of every snapshot. Returns the futures
    so a caller can track progress. Already-loaded / in-flight snapshots
    are not resubmitted."""
    # Silence per-snapshot load chatter so it doesn't bury the startup
    # browser URL in the terminal.
    from visage.io import halo_reader, galaxy_reader

    halo_reader.VERBOSE = False
    galaxy_reader.VERBOSE = False
    n = self._snap_table.count
    futures: list[Future] = []
    with self._lock:
        for snap in range(n):
            if snap not in self._futures:
                self._futures[snap] = self._submit(snap)
            futures.append(self._futures[snap])
    return futures

Utilities

visage.utils.colormap

cmap_css_gradient(name, n=12)

CSS linear-gradient string for a matplotlib colormap.

Source code in visage/utils/colormap.py

def cmap_css_gradient(name: str, n: int = 12) -> str:
    """CSS linear-gradient string for a matplotlib colormap."""
    import matplotlib.pyplot as plt

    cmap = plt.get_cmap(name)
    stops = []
    for i in range(n):
        t = i / (n - 1)
        r, g, b, a = cmap(float(t))
        stops.append(
            f"rgba({int(r*255)},{int(g*255)},{int(b*255)},{a:.2f}) {int(t*100)}%"
        )
    return "linear-gradient(to right, " + ", ".join(stops) + ")"

colormap_lut(name, n=256) cached

(n, 4) uint8 RGBA lookup table for a matplotlib colormap.

Bin centres, matching how VTK samples its own lookup table, so the colours come out identical to the ones VTK would have produced.

Source code in visage/utils/colormap.py

@lru_cache(maxsize=64)
def colormap_lut(name: str, n: int = 256) -> np.ndarray:
    """(n, 4) uint8 RGBA lookup table for a matplotlib colormap.

    Bin centres, matching how VTK samples its own lookup table, so the
    colours come out identical to the ones VTK would have produced.
    """
    import matplotlib

    cmap = matplotlib.colormaps[name]
    return (cmap((np.arange(n) + 0.5) / n) * 255.0).astype(np.uint8)

normalize_log(values, vmin, vmax)

Generic log10 normalisation to [0, 1].

Source code in visage/utils/colormap.py

def normalize_log(
    values: np.ndarray,
    vmin: float,
    vmax: float,
) -> np.ndarray:
    """Generic log10 normalisation to [0, 1]."""
    log_v = np.log10(np.maximum(values, 1e-30))
    return np.clip((log_v - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0).astype(
        np.float32
    )

normalize_log_halo_mass(mass, vmin=HALO_MASS_RANGE[0], vmax=HALO_MASS_RANGE[1])

Map halo masses (Msun) to [0, 1] via log10.

Source code in visage/utils/colormap.py

def normalize_log_halo_mass(
    mass: np.ndarray,
    vmin: float = HALO_MASS_RANGE[0],
    vmax: float = HALO_MASS_RANGE[1],
) -> np.ndarray:
    """Map halo masses (Msun) to [0, 1] via log10."""
    log_m = np.log10(np.maximum(mass, 1.0))
    return np.clip((log_m - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0)

normalize_log_mass(mass, vmin=STELLAR_MASS_RANGE[0], vmax=STELLAR_MASS_RANGE[1])

Map stellar or halo masses (Msun) to [0, 1] via log10.

Source code in visage/utils/colormap.py

def normalize_log_mass(
    mass: np.ndarray,
    vmin: float = STELLAR_MASS_RANGE[0],
    vmax: float = STELLAR_MASS_RANGE[1],
) -> np.ndarray:
    """Map stellar or halo masses (Msun) to [0, 1] via log10."""
    log_m = np.log10(np.maximum(mass, 1.0))
    return np.clip((log_m - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0)

normalize_log_ssfr(ssfr, vmin=SSFR_RANGE[0], vmax=SSFR_RANGE[1])

Map specific SFR (yr^-1) to [0, 1] via log10.

Source code in visage/utils/colormap.py

def normalize_log_ssfr(
    ssfr: np.ndarray,
    vmin: float = SSFR_RANGE[0],
    vmax: float = SSFR_RANGE[1],
) -> np.ndarray:
    """Map specific SFR (yr^-1) to [0, 1] via log10."""
    ssfr_safe = np.maximum(ssfr, 1e-14)
    log_ssfr = np.log10(ssfr_safe)
    return np.clip((log_ssfr - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0)

scalars_to_rgba(values, cmap, n=256)

Map already-normalised (0-1) scalars to RGBA, vectorised.

Handing VTK a scalar array plus a colormap makes PyVista convert the values to colours one at a time in Python — hundreds of milliseconds for a large point cloud, and it happens on every redraw. A table lookup here is the same 256-entry quantisation VTK's own lookup table uses, done in numpy.

Alpha is left at full: layer opacity stays an actor property, so the opacity slider can still move it without rebuilding anything.

Source code in visage/utils/colormap.py

def scalars_to_rgba(values: np.ndarray, cmap: str, n: int = 256) -> np.ndarray:
    """Map already-normalised (0-1) scalars to RGBA, vectorised.

    Handing VTK a scalar array plus a colormap makes PyVista convert the
    values to colours one at a time in Python — hundreds of milliseconds
    for a large point cloud, and it happens on every redraw.  A table
    lookup here is the same 256-entry quantisation VTK's own lookup table
    uses, done in numpy.

    Alpha is left at full: layer opacity stays an actor property, so the
    opacity slider can still move it without rebuilding anything.
    """
    v = np.asarray(values, dtype=np.float32)
    idx = np.clip((v * n).astype(np.int32), 0, n - 1)
    return colormap_lut(cmap, n)[idx]

visage.utils.sizing

galaxy_point_sizes(stellar_mass, size_min=HALO_SIZE_BINS[0] * GALAXY_SIZE_SCALE, size_max=HALO_SIZE_BINS[-1] * GALAXY_SIZE_SCALE, mass_range=STELLAR_MASS_RANGE)

Scale galaxy point sizes by stellar mass (Msun) using a fixed log10 range.

Source code in visage/utils/sizing.py

def galaxy_point_sizes(
    stellar_mass: np.ndarray,
    size_min: float = HALO_SIZE_BINS[0] * GALAXY_SIZE_SCALE,
    size_max: float = HALO_SIZE_BINS[-1] * GALAXY_SIZE_SCALE,
    mass_range: tuple[float, float] = STELLAR_MASS_RANGE,
) -> np.ndarray:
    """Scale galaxy point sizes by stellar mass (Msun) using a fixed log10 range."""
    if len(stellar_mass) == 0:
        return np.array([], dtype=np.float32)
    log_m = np.log10(np.maximum(stellar_mass, 1.0))
    vmin, vmax = mass_range
    norm = np.clip((log_m - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0)
    return (size_min + norm * (size_max - size_min)).astype(np.float32)

galaxy_world_radii(stellar_mass, r_min=0.025, r_max=0.35, mass_range=STELLAR_MASS_RANGE)

Per-galaxy world-space gaussian radius (Mpc/h), scaling with stellar mass.

Source code in visage/utils/sizing.py

def galaxy_world_radii(
    stellar_mass: np.ndarray,
    r_min: float = 0.025,  # Mpc/h  (low-mass end)
    r_max: float = 0.35,  # Mpc/h  (high-mass end)
    mass_range: tuple[float, float] = STELLAR_MASS_RANGE,
) -> np.ndarray:
    """Per-galaxy world-space gaussian radius (Mpc/h), scaling with stellar mass."""
    if len(stellar_mass) == 0:
        return np.array([], dtype=np.float32)
    log_m = np.log10(np.maximum(stellar_mass, 1.0))
    vmin, vmax = mass_range
    norm = np.clip((log_m - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0)
    return (r_min + norm * (r_max - r_min)).astype(np.float32)

galaxy_world_radii_rvir(rvir, mvir_raw=None, scale=1.0, r_min=0.02)

Per-galaxy world-space gaussian radius (Mpc/h): the subhalo Rvir.

The whole layered galaxy composition (outer CGM/Hot envelope, cold-gas envelope, focus-only disk/bulge) is keyed off this radius, so a galaxy's rendered extent is its subhalo's virial radius rather than a fixed stellar-mass mapping.

Where Rvir is missing (<= 0), falls back to the analytic virial radius of mvir_raw (raw SAGE units, 10^10 Msun/h; Δ=200 × ρ_crit), then to the r_min visibility floor.

Source code in visage/utils/sizing.py

def galaxy_world_radii_rvir(
    rvir: np.ndarray,
    mvir_raw: np.ndarray | None = None,
    scale: float = 1.0,
    r_min: float = 0.02,  # Mpc/h visibility floor
) -> np.ndarray:
    """Per-galaxy world-space gaussian radius (Mpc/h): the subhalo Rvir.

    The whole layered galaxy composition (outer CGM/Hot envelope, cold-gas
    envelope, focus-only disk/bulge) is keyed off this radius, so a galaxy's
    rendered extent is its subhalo's virial radius rather than a fixed
    stellar-mass mapping.

    Where Rvir is missing (<= 0), falls back to the analytic virial radius
    of ``mvir_raw`` (raw SAGE units, 10^10 Msun/h; Δ=200 × ρ_crit), then to
    the ``r_min`` visibility floor.
    """
    if len(rvir) == 0:
        return np.array([], dtype=np.float32)
    r = np.asarray(rvir, dtype=np.float64).copy()
    bad = ~(r > 0)
    if mvir_raw is not None and np.any(bad):
        m = np.maximum(  # Msun/h
            np.asarray(mvir_raw, dtype=np.float64)[bad] * 1.0e10, 1.0
        )
        rho_crit = 2.775e11  # h^2 Msun / Mpc^3
        r[bad] = (3.0 * m / (4.0 * np.pi * 200.0 * rho_crit)) ** (1.0 / 3.0)
    return np.maximum(r * scale, r_min).astype(np.float32)

glow_crowding(positions, sigma)

How many galaxy glows overlap at each galaxy, on average over viewing directions (1 = only its own).

A splat is a gaussian of width sigma. Seen from a random direction, a neighbour at 3D distance d adds, on average, D(x) / x of its peak brightness at this galaxy's position on screen, where x = d / (sigma * sqrt(2)) and D is Dawson's integral (the isotropic average of the projected gaussian). Averaging over directions keeps the estimate the same however the camera turns, so it is computed once per snapshot.

Source code in visage/utils/sizing.py

def glow_crowding(positions: np.ndarray, sigma: np.ndarray) -> np.ndarray:
    """How many galaxy glows overlap at each galaxy, on average over
    viewing directions (1 = only its own).

    A splat is a gaussian of width ``sigma``.  Seen from a random direction,
    a neighbour at 3D distance d adds, on average, D(x) / x of its peak
    brightness at this galaxy's position on screen, where x = d / (sigma *
    sqrt(2)) and D is Dawson's integral (the isotropic average of the
    projected gaussian).  Averaging over directions keeps the estimate the
    same however the camera turns, so it is computed once per snapshot.
    """
    from scipy.spatial import cKDTree
    from scipy.special import dawsn

    pos = np.asarray(positions, dtype=np.float64)
    sig = np.maximum(np.asarray(sigma, dtype=np.float64), 1e-6)
    n = len(pos)
    out = np.ones(n)
    if n < 2:
        return out
    tree = cKDTree(pos)
    # Neighbours are found per band of splat size, each band searched out
    # to the reach of its largest splat: one search at the reach of the
    # largest splat overall (a cluster's) would test every pair in the box.
    lo, hi = np.log(sig.min()), np.log(sig.max())
    edges = np.exp(np.linspace(lo, hi, max(2, int((hi - lo) / 0.4) + 2)))
    edges[-1] = np.inf
    band = np.searchsorted(edges, sig, side="right") - 1
    for b in np.unique(band):
        src = np.flatnonzero(band == b)
        reach = _CROWD_REACH * sig[src].max()
        pairs = cKDTree(pos[src]).sparse_distance_matrix(
            tree, reach, output_type="ndarray"
        )
        j, i, d = src[pairs["i"]], pairs["j"], pairs["v"]
        keep = (i != j) & (d < _CROWD_REACH * sig[j])
        j, i, d = j[keep], i[keep], d[keep]
        x = np.maximum(d / (np.sqrt(2.0) * sig[j]), 1e-6)
        out += np.bincount(i, dawsn(x) / x, minlength=n)
    return out

halo_point_sizes(masses, size_min=HALO_SIZE_BINS[0], size_max=HALO_SIZE_BINS[-1], mass_range=HALO_MASS_RANGE)

Scale halo point sizes by Mvir (Msun) using a fixed log10 range.

Fixed range prevents flickering when the mass distribution changes between snapshots.

Source code in visage/utils/sizing.py

def halo_point_sizes(
    masses: np.ndarray,
    size_min: float = HALO_SIZE_BINS[0],
    size_max: float = HALO_SIZE_BINS[-1],
    mass_range: tuple[float, float] = HALO_MASS_RANGE,
) -> np.ndarray:
    """Scale halo point sizes by Mvir (Msun) using a fixed log10 range.

    Fixed range prevents flickering when the mass distribution changes
    between snapshots.
    """
    if len(masses) == 0:
        return np.array([], dtype=np.float32)
    log_m = np.log10(np.maximum(masses, 1.0))
    vmin, vmax = mass_range
    norm = np.clip((log_m - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0)
    return (size_min + norm * (size_max - size_min)).astype(np.float32)

halo_world_radii(masses, r_min=0.15, r_max=1.5, mass_range=HALO_MASS_RANGE)

Per-halo world-space gaussian radius (Mpc/h), scaling with Mvir.

Source code in visage/utils/sizing.py

def halo_world_radii(
    masses: np.ndarray,
    r_min: float = 0.15,  # Mpc/h  (low-mass end)
    r_max: float = 1.5,  # Mpc/h  (high-mass end)
    mass_range: tuple[float, float] = HALO_MASS_RANGE,
) -> np.ndarray:
    """Per-halo world-space gaussian radius (Mpc/h), scaling with Mvir."""
    if len(masses) == 0:
        return np.array([], dtype=np.float32)
    log_m = np.log10(np.maximum(masses, 1.0))
    vmin, vmax = mass_range
    norm = np.clip((log_m - vmin) / (vmax - vmin + 1e-10), 0.0, 1.0)
    return (r_min + norm * (r_max - r_min)).astype(np.float32)

size_bin_mask(sizes, bin_edges)

Return a boolean mask per size bin for split-bin rendering.

Source code in visage/utils/sizing.py

def size_bin_mask(
    sizes: np.ndarray, bin_edges: list[float]
) -> list[np.ndarray]:
    """Return a boolean mask per size bin for split-bin rendering."""
    masks = []
    for i, edge in enumerate(bin_edges):
        lo = bin_edges[i - 1] if i > 0 else 0.0
        if i < len(bin_edges) - 1:
            masks.append((sizes >= lo) & (sizes < edge))
        else:
            masks.append(sizes >= lo)
    return masks

visage.utils.kdtree.NearestHaloIndex

Thin wrapper around scipy KDTree for nearest-halo spatial queries.

Rebuilt automatically when positions change (i.e. on snapshot change).

Source code in visage/utils/kdtree.py

class NearestHaloIndex:
    """Thin wrapper around scipy KDTree for nearest-halo spatial queries.

    Rebuilt automatically when positions change (i.e. on snapshot change).
    """

    def __init__(self, positions: np.ndarray | None = None) -> None:
        self._tree: KDTree | None = None
        self._positions: np.ndarray | None = None
        if positions is not None and len(positions) > 0:
            self.update(positions)

    def update(
        self, positions: np.ndarray, tree: KDTree | None = None
    ) -> None:
        self._positions = positions
        self._tree = tree if tree is not None else KDTree(positions)

    def nearest(self, point: tuple[float, float, float]) -> int:
        """Return index of the closest halo/point to the given (x, y, z)."""
        if self._tree is None:
            raise RuntimeError("Index is empty — call update() first.")
        _, idx = self._tree.query(point)
        return int(idx)

    def within_radius(
        self, center: tuple[float, float, float], radius: float
    ) -> np.ndarray:
        """Return indices of all points within radius of center."""
        if self._tree is None:
            return np.array([], dtype=np.int64)
        return np.array(
            self._tree.query_ball_point(center, radius), dtype=np.int64
        )

    def position_of(self, idx: int) -> np.ndarray:
        if self._positions is None:
            raise RuntimeError("Index is empty — call update() first.")
        return self._positions[idx]

nearest(point)

Return index of the closest halo/point to the given (x, y, z).

Source code in visage/utils/kdtree.py

def nearest(self, point: tuple[float, float, float]) -> int:
    """Return index of the closest halo/point to the given (x, y, z)."""
    if self._tree is None:
        raise RuntimeError("Index is empty — call update() first.")
    _, idx = self._tree.query(point)
    return int(idx)

within_radius(center, radius)

Return indices of all points within radius of center.

Source code in visage/utils/kdtree.py

def within_radius(
    self, center: tuple[float, float, float], radius: float
) -> np.ndarray:
    """Return indices of all points within radius of center."""
    if self._tree is None:
        return np.array([], dtype=np.int64)
    return np.array(
        self._tree.query_ball_point(center, radius), dtype=np.int64
    )

App

visage.app.create_app(par_path, par_dir=None, initial_snap=None, n_jobs=-1, min_halo_mass=10000000000.0, min_stellar_mass=100000000.0, max_halos=None, max_galaxies=None, cache_gb=None, port=8080, lightcone_path=None, perf=False, min_particles=None)

Source code in visage/app.py

def create_app(
    par_path: str | Path,
    par_dir: str | Path | None = None,
    initial_snap: int | None = None,
    n_jobs: int = -1,
    min_halo_mass: float = 1.0e10,
    min_stellar_mass: float = 1.0e8,
    max_halos: int | None = None,
    max_galaxies: int | None = None,
    cache_gb: float | None = None,
    port: int = 8080,
    lightcone_path: str | Path | None = None,
    perf: bool = False,
    min_particles: int | None = None,
):
    scene = Scene(
        primary_par_path=par_path,
        off_screen=system.prefer_offscreen_rendering(),
        initial_snap=initial_snap,
        n_jobs=n_jobs,
        min_halo_mass=min_halo_mass,
        min_stellar_mass=min_stellar_mass,
        max_halos=max_halos,
        max_galaxies=max_galaxies,
        cache_bytes=(
            int(cache_gb * 1024**3) if cache_gb is not None else None
        ),
        lightcone_path=lightcone_path,
        min_particles=min_particles,
    )

    server = get_server(client_type="vue3")
    server.enable_module(has_capabilities)

    # Say which GPU actually renders, and warn when OpenGL has fallen back
    # to the CPU (Remote Desktop, a VM, missing drivers) — that, not
    # ViSAGE, is what makes the view sluggish.
    def _schedule_hd_warm_up() -> None:
        """A few seconds after start: spawn HD's build workers and compile
        its shaders, so the first HD view doesn't wait for either."""
        try:
            loop = asyncio.get_running_loop()
        except RuntimeError:
            return

        async def _warm() -> None:
            from visage.hd.group import warm_pool
            from visage.hd.layer import warm_up_prepare, warm_up_shaders

            try:
                warm_pool()
                prep = await loop.run_in_executor(None, warm_up_prepare)
                # Compiling holds the main thread (0.4 s here, over a
                # second on some PCs): wait until no frames have been sent
                # for a couple of seconds, so it doesn't stall a drag.
                while time.monotonic() - frame_delivery.last_frame[0] < 2.0:
                    await asyncio.sleep(0.5)
                st = server.state
                if not (st.hd_active or st.hd_building):
                    warm_up_shaders(scene.plotter, prep)
            except Exception as exc:  # only a head start; HD works without
                print(f"HD warm-up skipped: {exc}", flush=True)

        loop.call_later(3.0, lambda: loop.create_task(_warm()))

    def _on_renderer(info: system.GPUInfo) -> None:
        print(f"GPU      : {info.describe()}", flush=True)
        hd_reason = system.hd_unavailable_reason(info)
        if hd_reason:
            # No GPU: HD mode is off (its button greyed out, saying why).
            server.state.hd_available = False
            server.state.hd_unavailable = hd_reason
            print(f"HD mode  : off -- {hd_reason}", flush=True)
        else:
            _schedule_hd_warm_up()
        if not info.software:
            return
        msg = (
            "Rendering on the CPU (software OpenGL) — the view will be slow. "
            "Update the GPU driver, or run ViSAGE locally rather than over "
            "Remote Desktop."
        )
        print(f"WARNING  : {msg}", flush=True)
        try:
            loop = asyncio.get_running_loop()
        except RuntimeError:
            return

        def _notify() -> None:
            server.state.notice_text = msg
            server.state.notice_color = "warning"
            server.state.notice_timeout = 15000
            server.state.notice_show = True
            server.state.flush()

        loop.call_soon(_notify)

    system.report_renderer_once(scene.plotter.ren_win, _on_renderer)
    frame_delivery.install()

    async def _frame_sent() -> None:
        # Two frames may still be on their way while the next is drawn; wait
        # only when more are queued, so the backlog stays at two frames.
        monitor = server.context.network_monitor
        while monitor.pending > 2 or frame_delivery.browser_behind():
            await asyncio.sleep(0.002)

    camera_motion.set_frame_gate(_frame_sent)
    # Mouse drags too: frames beyond two in flight are skipped, not queued.
    frame_delivery.backlog[0] = lambda: server.context.network_monitor.pending

    # ... and the browser's own backlog: a ping after each frame that the
    # page answers once it has caught up (visage.js).
    server.state.visage_ping = 0

    def _send_ping(ping_id: int) -> None:
        def _go() -> None:
            server.state.visage_ping = ping_id
            server.state.flush()

        try:
            asyncio.get_running_loop().call_soon(_go)
        except RuntimeError:
            frame_delivery._ping["sent"] = None

    frame_delivery.ping_sender[0] = _send_ping

    @server.controller.trigger("visage_pong")
    def _on_pong(ping_id):
        frame_delivery.pong(ping_id)

    # While frames flow, snapshot preloading waits rather than competing
    # with the renderer for the CPU.
    frame_delivery.on_frame.append(lambda: hold_background(0.3))

    # A lightcone is framed to the window's shape, which is only known once
    # the browser sends its first frame size: frame it again then (once).
    if scene.is_lightcone:
        _fit_size = [tuple(scene.plotter.ren_win.GetSize()[:2])]

        def _refit_lightcone() -> None:
            size = tuple(scene.plotter.ren_win.GetSize()[:2])
            if _fit_size[0] is None or size == _fit_size[0]:
                return
            _fit_size[0] = None

            def _go() -> None:
                scene.camera.reset()
                if hasattr(server.controller, "view_update"):
                    server.controller.view_update()

            try:
                asyncio.get_running_loop().call_soon(_go)
            except RuntimeError:
                pass

        frame_delivery.on_frame.append(_refit_lightcone)
    if perf:
        from visage.utils import perf as _perf

        _perf.install(
            scene.plotter.ren_win, out=lambda s: print(s, flush=True)
        )

    # Serve ViSAGE's client-side helpers (pop-out drag, Enter-to-
    # click). Vue 3 silently drops <script> tags from templates so we
    # have to inject these via the module/static-asset system.
    import os as _os_static

    _sage_static_dir = _os_static.path.join(
        _os_static.path.dirname(__file__), "static"
    )

    # Cache-bust the served JS/CSS with each file's mtime so browsers
    # pick up edits on restart instead of replaying a stale cached copy
    # (which silently breaks client-side features like the pop-out
    # fullscreen toggle until a manual hard refresh).
    def _bust(rel):
        # Map "sage_static/<subpath>" to the file under the static dir so
        # nested assets (e.g. katex/katex.min.js) are cache-busted too.
        sub = rel.split("?", 1)[0]
        if sub.startswith("sage_static/"):
            sub = sub[len("sage_static/") :]
        try:
            mtime = int(
                _os_static.path.getmtime(
                    _os_static.path.join(_sage_static_dir, sub)
                )
            )
        except OSError:
            return rel
        return f"{rel}?v={mtime}"

    server.enable_module(
        {
            "serve": {"sage_static": _sage_static_dir},
            # KaTeX (vendored, woff2-only) must load before visage.js,
            # which calls katex.render() for Story Mode equation overlays.
            "scripts": [
                _bust("sage_static/katex/katex.min.js"),
                _bust("sage_static/visage.js"),
            ],
            "styles": [
                _bust("sage_static/katex/katex.min.css"),
                _bust("sage_static/sage_theme.css"),
            ],
        }
    )

    # xterm.js — browser-side terminal emulator for the console PTY.
    server.enable_module(
        {
            "styles": [
                "https://cdn.jsdelivr.net/npm/xterm@5.3.0/css/xterm.css",
            ],
            "scripts": [
                "https://cdn.jsdelivr.net/npm/xterm@5.3.0/lib/xterm.js",
                "https://cdn.jsdelivr.net/npm/@xterm/addon-fit@0.10.0/lib/addon-fit.js",
                # Canvas renderer: block characters (SAGE26's banner) fill
                # their cells, as in a normal terminal.
                "https://cdn.jsdelivr.net/npm/xterm-addon-canvas@0.5.0/lib/xterm-addon-canvas.js",
            ],
        }
    )

    # Single-theme config — DOS Blue is now the only palette.
    _vuetify_config = {
        "theme": {
            "defaultTheme": "dos_blue",
            "themes": {
                "dos_blue": {
                    "dark": True,
                    "colors": {
                        "primary": "#ffff55",  # DOS yellow
                        "secondary": "#ffffff",
                        "background": "#000000",
                        "surface": "#000000",
                        "on-surface": "#ffffff",
                        "on-background": "#ffffff",
                    },
                },
            },
        }
    }
    # Story Mode player + reactive state / controllers.
    init_story_mode(server, scene)

    _NAV_TABS = [
        ("Structure", "layers"),
        ("Filters", "filters"),
        ("Record", "record"),
        ("Coords", "coords"),
        ("Box", "box"),
        ("Console", "console"),
        ("Library", "library"),
    ]

    # ---- Model discovery ------------------------------------------------
    # We list models by scanning the SAGE OUTPUT directory (each subfolder
    # with a model_0.hdf5 is one model, named after the folder).  The .par
    # files are still used internally for tree paths but they are no longer
    # surfaced in the UI.
    # A lightcone is a single standalone cloud — there are no sibling SAGE
    # model folders to discover, so skip the scan (and the par-based paths,
    # which don't apply).
    if lightcone_path is not None:
        discovered = []
    else:
        primary_hdf5 = Path(scene.primary.cfg.hdf5_path)
        output_dir = primary_hdf5.parent.parent
        par_d = Path(par_dir) if par_dir else Path(par_path).parent
        discovered = find_models(output_dir, par_dir=par_d)
    # Index by model name for fast lookup
    discovered_by_name: dict[str, dict] = {m["name"]: m for m in discovered}

    # Per-box profile storage (Python side only; keys = model names).
    _profiles: dict = {}

    def _build_box_strip_items() -> list[dict]:
        """Build the reactive list for the viewport box strip."""
        strip = []
        pm = scene.primary
        s_lbl = pm.snap_table.label(max(0, pm.current_snap))
        strip.append(
            {
                "name": pm.name,
                "label": f"{pm.name}  {s_lbl}",
                "active": scene.active_box_name == pm.name,
                "primary": True,
            }
        )
        for adj_name in scene._adjacent_order:
            m = scene._models.get(adj_name)
            if m is None:
                continue
            al = m.snap_table.label(max(0, m.current_snap))
            strip.append(
                {
                    "name": adj_name,
                    "label": f"{adj_name}  {al}",
                    "active": scene.active_box_name == adj_name,
                    "primary": False,
                }
            )
        return strip

    def _build_models_list() -> list[dict]:
        """Build the reactive list of model entries for the menu."""
        loaded = {m.name: m for m in scene.list_models()}
        out = []
        for entry in discovered:
            name = entry["name"]
            is_loaded = name in loaded
            is_primary = name == scene.primary_name
            is_adjacent = is_loaded and scene.is_adjacent(name)
            compatible = is_loaded and scene.is_compatible_for_overlay(name)
            overlay_on = (
                is_loaded
                and not is_primary
                and not is_adjacent
                and loaded[name].visible
            )
            out.append(
                {
                    "name": name,
                    "path": str(entry["par"]),
                    "loaded": is_loaded,
                    "primary": is_primary,
                    "compatible": compatible,
                    "overlay": overlay_on,
                    "adjacent": is_adjacent,
                }
            )
        # Loaded-but-not-on-disk (e.g. primary model whose output dir wasn't scanned)
        for name, m in loaded.items():
            if not any(e["name"] == name for e in out):
                is_adjacent = scene.is_adjacent(name)
                out.append(
                    {
                        "name": name,
                        "path": str(m.path),
                        "loaded": True,
                        "primary": name == scene.primary_name,
                        "compatible": scene.is_compatible_for_overlay(name),
                        "overlay": m.visible
                        and name != scene.primary_name
                        and not is_adjacent,
                        "adjacent": is_adjacent,
                    }
                )
        return out

    server.state.models_list = _build_models_list()

    # ── Session models registry ─────────────────────────────────────────
    # Remember every model (box or lightcone) opened this session so the user
    # can jump back to any of them from the Launch dropdown.  Re-opening one
    # relaunches ViSAGE on it (like the wizard's own Explore/Visualize launch).
    from visage.utils import recent_models as _recent

    if lightcone_path is not None:
        _cur_kind = "lightcone"
        _cur_path = str(Path(lightcone_path).expanduser().resolve())
        _cur_name = Path(lightcone_path).stem
    else:
        _cur_kind = "box"
        _cur_path = (
            str(Path(par_path).expanduser().resolve()) if par_path else ""
        )
        _cur_name = scene.primary_name
    _sess_entries = (
        _recent.record(_cur_name, _cur_kind, _cur_path)
        if _cur_path
        else _recent.load()
    )
    server.state.session_models = [
        {
            "name": e.get("name", ""),
            "kind": e.get("kind", "box"),
            "path": e.get("path", ""),
            # Shown as the visible subtitle (full path is the hover tooltip)
            # so two entries with the same name/kind — e.g. two lightcones —
            # are actually distinguishable in the dropdown.
            "parent": str(Path(e.get("path", "")).parent),
            "active": e.get("kind") == _cur_kind
            and e.get("path") == _cur_path,
        }
        for e in _sess_entries
    ]

    server.state.relaunching = ""  # name of the model being opened

    @server.controller.set("open_recent_model")
    async def _open_recent_model(path, kind):
        import asyncio as _asyncio

        argv = (
            ["--lightcone", str(path)]
            if kind == "lightcone"
            else ["--par", str(path)]
        ) + ["--port", str(port)]
        # Show the switching card and have the page reload itself once
        # the new ViSAGE is up (visage.js watches `relaunching`).
        server.state.model_loading = True
        server.state.relaunching = Path(str(path)).stem or "model"
        server.state.flush()
        await _asyncio.sleep(0.8)
        system.relaunch_visage(argv)

    server.state.active_box_name = scene.primary_name
    server.state.box_strip_items = []
    server.state.model_loading = False
    # Rotating quip shown on the "switching models" overlay. Updated by
    # an asyncio task that runs while model_loading is True.
    server.state.model_quip = "Switching models, please hold..."
    _MODEL_QUIPS: list[str] = [
        "Reticulating splines...",
        "Herding electrons into formation...",
        "Asking the universe nicely for the haloes...",
        "Spinning up galaxies — please don't shake the box.",
        "Negotiating with dark matter (it drives a hard bargain).",
        "Counting black holes... lost count.",
        "Convincing photons to travel faster. No luck.",
        "Reading SAGE bedtime stories to the trees...",
        "Stirring the cold gas. Gently.",
        "Polishing the CGM. Won't be long.",
        "Aligning angular momenta — close your eyes.",
        "Subhalo abundance matching the vibes...",
        "Hubble tension intensifies...",
        "Letting the H2 form on dust grains. Patience.",
        "Renormalizing the friend-of-friend friendships.",
        "Tracing merger trees — they're surprisingly deep.",
        "Reminding satellites who's central.",
        "Quenching star formation (sorry).",
        "Calibrating feedback — too much, less, more, less...",
        "Recomputing 1/H(z). Again.",
        "Reading the par file out loud, slowly.",
    ]

    async def _quip_loop():
        import asyncio as _asyncio
        import random as _random

        try:
            while bool(server.state.model_loading):
                server.state.model_quip = _random.choice(_MODEL_QUIPS)
                server.state.flush()
                await _asyncio.sleep(2.5)
        finally:
            server.state.model_quip = "Switching models, please hold..."
            server.state.flush()

    _quip_task: list = [None]

    @server.state.change("model_loading")
    def on_model_loading_change(model_loading, **_):
        import asyncio as _asyncio

        if model_loading:
            if _quip_task[0] is None or _quip_task[0].done():
                _quip_task[0] = _asyncio.ensure_future(_quip_loop())
        else:
            if _quip_task[0] is not None and not _quip_task[0].done():
                _quip_task[0].cancel()

    server.state.model_fields = scene.primary.fields_available
    server.state.primary_model = scene.primary.name
    # Snackbar for overlay-compatibility errors etc.
    server.state.notice_show = False
    server.state.notice_text = ""
    server.state.notice_color = "warning"
    server.state.notice_timeout = 4500

    # Update checker
    server.state.update_checking = False

    # Per-model flags used by static menu items (dict keyed by name)
    def _model_flags() -> dict:
        loaded = {m.name: m for m in scene.list_models()}
        out = {}
        for entry in discovered:
            n = entry["name"]
            is_adjacent = scene.is_adjacent(n) if n in loaded else False
            out[n] = {
                "primary": n == scene.primary.name,
                "loaded": n in loaded,
                "overlay": n in loaded
                and loaded[n].visible
                and n != scene.primary.name
                and not is_adjacent,
                "compatible": (
                    scene.is_compatible_for_overlay(n) if n in loaded else True
                ),
                "adjacent": is_adjacent,
            }
        return out

    server.state.model_flags = _model_flags()

    def _refresh_models_state() -> None:
        server.state.models_list = _build_models_list()
        server.state.model_fields = scene.primary.fields_available
        server.state.primary_model = scene.primary.name
        server.state.model_flags = _model_flags()
        server.state.box_strip_items = _build_box_strip_items()
        server.state.flush()

    scene.register_model_change_callback(_refresh_models_state)

    @server.controller.set("switch_model")
    async def on_switch_model(name: str):
        import asyncio

        # Capture the outgoing model's full UI state (filters, structure,
        # colormaps, opacity, visibility, snapshot) so it carries over to
        # the model we switch to.  Redshift — not the raw snap index — is
        # what we preserve, so models with different snapshot lists land on
        # the matching cosmic time rather than the matching slider position.
        carried = save_profile(server.state)
        old_table = scene.primary.snap_table
        cur_snap = int(carried.get("snap_num") or 0)
        cur_snap = max(0, min(cur_snap, old_table.count - 1))
        carried_z = old_table.snap_to_z(cur_snap)

        server.state.model_loading = True
        server.state.flush()
        await asyncio.sleep(
            0
        )  # yield so the browser receives model_loading=True
        try:
            if not scene.has_model(name):
                entry = discovered_by_name.get(name)
                if entry is None:
                    return
                scene.add_model(entry["par"])
            scene.switch_primary(name)
            # Preload all snapshots of the new primary in the background.
            server.controller.preload_snapshots(scene.primary.loader)
        finally:
            server.state.model_loading = False
            # Map the carried-over redshift onto the new model's snapshot
            # list (closest match), clamped to its valid range.
            new_max = scene.primary.snap_count - 1
            target_snap = scene.primary.snap_table.z_to_snap(carried_z)
            target_snap = max(0, min(target_snap, new_max))
            # Re-apply the carried settings to the new primary, with the
            # snapshot keys adjusted to the new model's range / matched z.
            carried["snap_num"] = target_snap
            carried["snap_max"] = new_max
            load_profile(server.state, carried)
            # dirty() forces the @state.change handlers to fire even when a
            # value is numerically unchanged, so every setting is re-applied
            # to the now-active (new primary) model's layers and filters.
            server.state.dirty(*BOX_PROFILE_KEYS)
            # Keep this model's stored profile in sync so a later
            # side-by-side activation restores the same carried-over state.
            _profiles[name] = save_profile(server.state)
            server.state.snap_label = scene.snap_label
            _refresh_models_state()
            server.state.flush()

    @server.controller.set("toggle_overlay")
    async def on_toggle_overlay(name: str):
        import asyncio

        server.state.model_loading = True
        server.state.flush()
        await asyncio.sleep(
            0
        )  # yield so the browser receives model_loading=True
        try:
            if not scene.has_model(name):
                entry = discovered_by_name.get(name)
                if entry is None:
                    return
                scene.add_model(entry["par"])
            model = scene._models[name]
            # Try the toggle; capture any compatibility error
            err = scene.set_overlay_visible(name, not model.visible)
            if err is not None:
                server.state.notice_text = err
                server.state.notice_color = "warning"
                server.state.notice_show = True
            elif model.visible:
                server.controller.preload_snapshots(model.loader)
        finally:
            server.state.model_loading = False
            _refresh_models_state()

    @server.controller.set("toggle_adjacent")
    async def on_toggle_adjacent(name: str):
        import asyncio

        entry = discovered_by_name.get(name)
        if entry is not None:
            par_path = entry["par"]
        elif scene.has_model(name):
            par_path = scene._models[name].path
        else:
            return
        server.state.model_loading = True
        server.state.flush()
        await asyncio.sleep(
            0
        )  # yield so the browser receives model_loading=True
        try:
            is_now_adj, err = scene.toggle_adjacent(par_path)
            if err:
                server.state.notice_text = err
                server.state.notice_color = "warning"
                server.state.notice_show = True
                return
            if is_now_adj:
                adj_m = scene._models.get(name)
                if adj_m and name not in _profiles:
                    _profiles[name] = default_profile(adj_m.snap_count)
                if adj_m:
                    server.controller.preload_snapshots(adj_m.loader)
                # Stop any active rotation — it can't be per-box with a shared camera
                if getattr(server.state, "rotate_mode", "off") != "off":
                    server.state.rotate_mode = "off"
                # Reframe to show all loaded boxes
                regions = [(0.0, 0.0, 0.0, scene.primary.box_size)]
                for adj_name in scene._adjacent_order:
                    m = scene._models.get(adj_name)
                    if m is not None:
                        off = m.offset
                        regions.append(
                            (
                                float(off[0]),
                                float(off[1]),
                                float(off[2]),
                                m.box_size,
                            )
                        )
                scene.camera.focus_on_boxes(regions)
            else:
                _profiles.pop(name, None)
                server.state.active_box_name = scene.active_box_name
        finally:
            server.state.model_loading = False
            server.state.box_strip_items = _build_box_strip_items()
            _refresh_models_state()
            if hasattr(server.controller, "view_update"):
                server.controller.view_update()

    @server.controller.set("set_active_box")
    def on_set_active_box(name: str):
        if not scene.has_model(name):
            return
        old_name = scene.active_box_name
        if old_name == name:
            return
        _profiles[old_name] = save_profile(server.state)
        # HD and the selection are of the old box's catalogue.
        for hook in ("hd_drop", "sel_reset"):
            try:
                getattr(server.controller, hook)()
            except Exception:
                pass
        scene.set_active_box(name)
        server.state.active_box_name = name
        incoming = _profiles.get(name)
        if incoming is None:
            m = scene._models.get(name)
            incoming = default_profile(m.snap_count) if m else {}
            _profiles[name] = incoming
        load_profile(server.state, incoming)
        server.state.dirty(*BOX_PROFILE_KEYS)
        server.state.box_strip_items = _build_box_strip_items()
        server.state.flush()
        if hasattr(server.controller, "sync_active_snap_count"):
            server.controller.sync_active_snap_count()
        if hasattr(server.controller, "view_update"):
            server.controller.view_update()

    @server.controller.set("clear_box")
    def on_clear_box(name: str):
        m = scene._models.get(name)
        if m is None:
            return
        fresh = default_profile(m.snap_count)
        _profiles[name] = fresh
        m.set_snapshot(m.snap_count - 1)
        m.halo_layer.set_filter_mask(None)
        m.galaxy_layer.set_filter_mask(None)
        if name == scene.active_box_name:
            load_profile(server.state, fresh)
            server.state.dirty(*BOX_PROFILE_KEYS)
            server.state.box_strip_items = _build_box_strip_items()
            server.state.flush()
        else:
            server.state.box_strip_items = _build_box_strip_items()
            server.state.flush()

    # ── Launch Mode wizard (embedded overlay) ────────────────────────────────
    from visage.wizard.controller import WizardController
    from visage.wizard.ui import build_wizard_ui
    from visage.wizard.pso_gallery import build_pso_gallery

    server.state.wiz_active = False

    _wiz_ctrl = WizardController(
        server,
        port=port,
        scene=scene,
        auto_start=False,
    )

    @server.controller.set("open_wizard")
    def _open_wizard():
        _wiz_ctrl.reset_and_start()
        server.state.wiz_active = True
        server.state.flush()

    # Each tool opens on its own first screen.
    @server.controller.set("open_sage26_wizard")
    def _open_sage26_wizard():
        _wiz_ctrl.reset_and_start(flow="sage26")
        server.state.wiz_active = True
        server.state.flush()

    @server.controller.set("open_sageswarm_wizard")
    def _open_sageswarm_wizard():
        _wiz_ctrl.reset_and_start(flow="sageswarm")
        server.state.wiz_active = True
        server.state.flush()

    @server.controller.set("open_sagelightcone_wizard")
    def _open_sagelightcone_wizard():
        _wiz_ctrl.reset_and_start(flow="sagelightcone")
        server.state.wiz_active = True
        server.state.flush()

    @server.controller.set("check_for_updates")
    async def _on_check_for_updates():
        import asyncio
        import json
        import sys
        import urllib.request

        from visage._version import __version__ as _current

        server.state.update_checking = True
        server.state.flush()
        await asyncio.sleep(0)

        loop = asyncio.get_event_loop()

        def _fetch_latest() -> str:
            import ssl

            try:
                import certifi

                ctx = ssl.create_default_context(cafile=certifi.where())
            except ImportError:
                ctx = ssl.create_default_context()
            try:
                with urllib.request.urlopen(
                    "https://pypi.org/pypi/sage-viewer/json",
                    timeout=8,
                    context=ctx,
                ) as r:
                    return json.loads(r.read())["info"]["version"]
            except ssl.SSLError:
                # macOS Python may lack system CA certs; retry unverified
                ctx = ssl._create_unverified_context()
                with urllib.request.urlopen(
                    "https://pypi.org/pypi/sage-viewer/json",
                    timeout=8,
                    context=ctx,
                ) as r:
                    return json.loads(r.read())["info"]["version"]

        try:
            latest = await loop.run_in_executor(None, _fetch_latest)
        except Exception as exc:
            server.state.notice_text = f"Update check failed: {exc}"
            server.state.notice_color = "error"
            server.state.notice_timeout = 5000
            server.state.notice_show = True
            server.state.update_checking = False
            server.state.flush()
            return

        if latest == _current:
            server.state.notice_text = f"ViSAGE {_current} is up to date."
            server.state.notice_color = "success"
            server.state.notice_timeout = 4500
            server.state.notice_show = True
            server.state.update_checking = False
            server.state.flush()
            return

        server.state.notice_text = f"Updating {_current} → {latest}…"
        server.state.notice_color = "info"
        server.state.notice_timeout = -1
        server.state.notice_show = True
        server.state.flush()

        try:
            proc = await asyncio.create_subprocess_exec(
                sys.executable,
                "-m",
                "pip",
                "install",
                "--upgrade",
                "sage-viewer",
                stdout=asyncio.subprocess.PIPE,
                stderr=asyncio.subprocess.PIPE,
                start_new_session=os.name != "nt",
            )
            from visage import procs

            procs.track(proc)  # stopped if ViSAGE closes mid-install
            try:
                _, stderr = await proc.communicate()
            finally:
                procs.untrack(proc)
            if proc.returncode == 0:
                server.state.notice_text = (
                    f"Updated to {latest}. Restart ViSAGE to apply."
                )
                server.state.notice_color = "success"
                server.state.notice_timeout = -1
            else:
                server.state.notice_text = (
                    f"pip failed: {stderr.decode(errors='replace')[:140]}"
                )
                server.state.notice_color = "error"
                server.state.notice_timeout = 7000
        except Exception as exc:
            server.state.notice_text = f"Install failed: {exc}"
            server.state.notice_color = "error"
            server.state.notice_timeout = 7000

        server.state.notice_show = True
        server.state.update_checking = False
        server.state.flush()

    # `theme=("ui_theme",)` reactively binds the active Vuetify theme to
    # our state variable — Vuetify swaps the entire palette plus the root
    # `v-theme--<name>` class instantly when ui_theme changes.
    server.state.ui_theme = "dos_blue"
    with SinglePageLayout(
        server,
        full_height=True,
        vuetify_config=_vuetify_config,
        theme=("ui_theme",),
        style="background:#000000;",
    ) as layout:
        # Hide the SinglePageLayout's auto-built title — we render our own
        # later inside the toolbar so we can control its position relative
        # to the hamburger menu.
        layout.title.style = "display:none;"
        # Hide the default VAppBarNavIcon — replaced by our custom menu button below
        layout.icon.style = "display:none;"

        with layout.toolbar as tb:
            tb.density = "compact"
            tb.color = "#000000"
            tb.elevation = 0

            # ── Launch Mode button — wizard + models ───────────────────────
            with v3.VMenu(close_on_content_click=True):
                with v3.Template(v_slot_activator="{ props }"):
                    with v3.VBtn(
                        variant="text",
                        density="compact",
                        v_bind="props",
                        title="Launch Mode",
                        style="padding:2px 4px;min-width:36px;",
                    ):
                        html.Img(
                            src="/sage_static/SAGElogo.jpg",
                            style=(
                                "height:30px;width:30px;"
                                "object-fit:cover;border-radius:50%;"
                            ),
                        )
                with v3.VList(density="compact", bg_color="transparent"):
                    # ── Launch Mode ────────────────────────────────────────
                    v3.VListSubheader(
                        "LAUNCH MODE",
                        style="color:#06b6d4;font-size:0.65rem;",
                    )
                    v3.VListItem(
                        title="SAGE26",
                        prepend_icon="mdi-cube-outline",
                        click=server.controller.open_sage26_wizard,
                        color="cyan",
                        density="compact",
                    )
                    v3.VListItem(
                        title="LightSAGE",
                        prepend_icon="mdi-telescope",
                        click=server.controller.open_sagelightcone_wizard,
                        color="cyan",
                        density="compact",
                    )
                    v3.VListItem(
                        title="SAGEswarm",
                        prepend_icon="mdi-chart-scatter-plot",
                        click=server.controller.open_sageswarm_wizard,
                        color="cyan",
                        density="compact",
                    )
                    v3.VListItem(
                        title="Setup Wizard",
                        prepend_icon="mdi-rocket-launch",
                        click=server.controller.open_wizard,
                        color="cyan",
                        density="compact",
                    )
                    # ── Session models (boxes + lightcones loaded this
                    #    session; click to re-open one) ──────────────────────
                    v3.VDivider(
                        v_show=(
                            "session_models && session_models.length > 0",
                        ),
                        style="margin:4px 0;",
                    )
                    v3.VListSubheader(
                        "SESSION MODELS",
                        style="color:#06b6d4;font-size:0.65rem;",
                        v_show=(
                            "session_models && session_models.length > 0",
                        ),
                    )
                    with html.Div(
                        v_for=("m in session_models",),
                        key=("'sess-' + m.path",),
                        title=("m.path",),  # full path on hover
                    ):
                        v3.VListItem(
                            title=("m.name",),
                            subtitle=("m.parent",),
                            prepend_icon=(
                                "m.kind === 'lightcone' "
                                "? 'mdi-telescope' : 'mdi-cube-outline'",
                            ),
                            click=(
                                server.controller.open_recent_model,
                                "[m.path, m.kind]",
                            ),
                            active=("m.active",),
                            color="cyan",
                            density="compact",
                        )
                    v3.VDivider(style="margin:4px 0;")
                    # ── Models (switch rows) ───────────────────────────────
                    v3.VListSubheader(
                        "MODELS",
                        style="color:#06b6d4;font-size:0.65rem;",
                        v_show=("models_list && models_list.length > 0",),
                    )
                    with html.Div(
                        v_for=(
                            "m in [...models_list].sort("
                            "(a,b) => (b.primary ? 1 : 0) - (a.primary ? 1 : 0))",
                        ),
                        key=("'sw-' + m.name",),
                    ):
                        v3.VListItem(
                            title=("m.name",),
                            prepend_icon=(
                                "m.primary ? 'mdi-check-circle' : 'mdi-circle-outline'",
                            ),
                            click=(server.controller.switch_model, "[m.name]"),
                            active=("m.primary",),
                            color="cyan",
                            density="compact",
                        )
                    # ── Divider + overlay rows ─────────────────────────────
                    v3.VDivider(
                        v_show=("models_list && models_list.length > 1",),
                        style="margin:4px 0;",
                    )
                    v3.VListSubheader(
                        "OVERLAYS",
                        style="color:#9ca3af;font-size:0.65rem;",
                        v_show=("models_list && models_list.length > 1",),
                    )
                    with html.Div(
                        v_for=("m in models_list",),
                        key=("'ov-' + m.name",),
                        v_show=("!m.primary",),
                    ):
                        v3.VListItem(
                            title=(
                                "m.overlay ? '✓ ' + m.name : '+ ' + m.name",
                            ),
                            prepend_icon=(
                                "m.overlay ? 'mdi-layers' : 'mdi-layers-plus'",
                            ),
                            click=(
                                server.controller.toggle_overlay,
                                "[m.name]",
                            ),
                            active=("m.overlay",),
                            color="cyan",
                            density="compact",
                        )
                    # ── Side by side rows ─────────────────────────────────
                    v3.VDivider(
                        v_show=("models_list && models_list.length > 1",),
                        style="margin:4px 0;",
                    )
                    v3.VListSubheader(
                        "SIDE BY SIDE",
                        style="color:#06b6d4;font-size:0.65rem;",
                        v_show=("models_list && models_list.length > 1",),
                    )
                    with html.Div(
                        v_for=("m in models_list",),
                        key=("'sb-' + m.name",),
                        v_show=("!m.primary",),
                    ):
                        v3.VListItem(
                            title=(
                                "m.adjacent ? '✓ ' + m.name : '⊞ ' + m.name",
                            ),
                            prepend_icon=(
                                "m.adjacent "
                                "? 'mdi-check-circle-outline' "
                                ": 'mdi-view-split-vertical'",
                            ),
                            click=(
                                server.controller.toggle_adjacent,
                                "[m.name]",
                            ),
                            active=("m.adjacent",),
                            color="cyan",
                            density="compact",
                        )
                    # ── Close application ──────────────────────────────────
                    v3.VDivider(style="margin:4px 0;")
                    v3.VListItem(
                        title="Close Everything",
                        prepend_icon="mdi-close-box-outline",
                        click=server.controller.close_app,
                        color="#ef4444",
                        density="compact",
                    )

            # ── Explore Mode menu (hamburger) — tabs.  Sits right next to the
            #    Launch Mode button. ─────────────────────────────────────────
            with v3.VMenu(close_on_content_click=True):
                with v3.Template(v_slot_activator="{ props }"):
                    v3.VBtn(
                        icon="mdi-menu",
                        variant="text",
                        density="compact",
                        v_bind="props",
                        title="Tabs",
                        style="margin-left:2px;",
                    )
                with v3.VList(density="compact", bg_color="transparent"):
                    for label, value in _NAV_TABS:
                        v3.VListItem(
                            title=label,
                            value=value,
                            click=f"nav_active_tab = '{value}'",
                            active=(f"nav_active_tab === '{value}'",),
                            color="cyan",
                        )

            # ── SAGE26 setup button — the box (opens the SAGE26 setup
            #    wizard). ─────────────────────────────────────────────────
            v3.VBtn(
                icon="mdi-cube-outline",
                variant="text",
                density="compact",
                color="white",
                title="SAGE26 Setup",
                click=server.controller.open_sage26_wizard,
                style="margin-left:2px;",
            )

            # ── SAGEswarm button — opens the wizard straight into the PSO
            #    calibration flow. ───────────────────────────────────────────
            v3.VBtn(
                icon="mdi-chart-scatter-plot",
                variant="text",
                density="compact",
                color="white",
                title="SAGEswarm — PSO Calibration",
                click=server.controller.open_sageswarm_wizard,
                style="margin-left:2px;",
            )

            # ── LightSAGE button — opens the wizard straight into the
            #    lightcone-extraction flow. ────────────────────────────────
            v3.VBtn(
                icon="mdi-telescope",
                variant="text",
                density="compact",
                color="white",
                title="LightSAGE — Lightcone Extraction",
                click=server.controller.open_sagelightcone_wizard,
                style="margin-left:2px;",
            )

            # ── Export catalogue button ────────────────────────────────────
            v3.VBtn(
                icon="mdi-database-export-outline",
                variant="text",
                density="compact",
                color="white",
                title="Export galaxy catalogue",
                click="export_dialog_show = true",
                style="margin-left:4px;",
            )

            # ── Fly-through button ─────────────────────────────────────────
            v3.VBtn(
                icon="mdi-rotate-orbit",
                variant="text",
                density="compact",
                color=("flythrough_active ? 'cyan' : 'white'",),
                title="Fly-through: approach box then orbit (click to stop)",
                click=server.controller.toggle_flythrough,
                style="margin-left:2px;",
            )

            # ── Story Mode button + dropdown ───────────────────────────────
            build_story_button(server)

            # ── Check / install updates ────────────────────────────────────
            v3.VBtn(
                icon="mdi-update",
                variant="text",
                density="compact",
                color="white",
                title="Check for updates",
                loading=("update_checking",),
                click=server.controller.check_for_updates,
                style="margin-left:2px;",
            )

            # Title
            v3.VToolbarTitle(
                "ViSAGE",
                style="padding-left:4px;",
            )

            build_toolbar(server, scene)

        with layout.content:
            # Pixel-style monospace for the retro palettes — loaded via a
            # real <link> tag so the browser treats it as a normal external
            # stylesheet (works even when injected inside a Vue template).
            html.Link(
                rel="stylesheet",
                href=(
                    "https://fonts.googleapis.com/css2?"
                    "family=VT323&display=swap"
                ),
            )
            # Vuetify DOS-blue theme overrides are loaded via enable_module
            # (sage_static/sage_theme.css) — that's the only path that reliably
            # reaches html/body from inside a Trame/Vue template.
            html.Style(_THEME_CSS)

            # ── Story Mode playback overlay (hidden unless a story is active)
            build_story_hud(server)

            # ── Export catalogue dialog ────────────────────────────────────
            # Scope options are built in navigation_panel (state var) so the
            # lightcone-only "Whole Lightcone" entry appears only in that mode.
            _FMT_ITEMS = [
                {"title": "CSV", "value": "csv"},
                {"title": "HDF5", "value": "hdf5"},
                {"title": "FITS", "value": "fits"},
                {"title": "TXT", "value": "txt"},
            ]
            with v3.VDialog(
                v_model=("export_dialog_show",),
                max_width=500,
                persistent=False,
            ):
                with v3.VCard(
                    style="background:transparent !important;border:1px solid #06b6d4;color:#e2e8f0;",
                    color="transparent",
                    elevation=0,
                    rounded=False,
                ):
                    with html.Div(
                        style=(
                            "display:flex;align-items:center;gap:8px;"
                            "padding:14px 16px 10px;"
                            "border-bottom:1px solid #374151;"
                        ),
                    ):
                        v3.VIcon(
                            "mdi-database-export-outline",
                            color="cyan",
                            size="small",
                        )
                        html.Span(
                            "Export Galaxy Catalogue",
                            style=(
                                "color:#06b6d4;font-weight:700;"
                                "letter-spacing:0.06em;font-size:0.95rem;"
                            ),
                        )
                        v3.VSpacer()
                        v3.VBtn(
                            icon="mdi-close",
                            size="x-small",
                            variant="text",
                            color="#9ca3af",
                            click="export_dialog_show = false",
                        )
                    with v3.VCardText(
                        style="padding:16px;display:flex;flex-direction:column;gap:14px;"
                    ):
                        # Scope
                        v3.VSelect(
                            v_model=("export_scope",),
                            items=("export_scope_items",),
                            label="Selection scope",
                            variant="outlined",
                            density="compact",
                            color="cyan",
                            bg_color="#0d0d1a",
                            hide_details=True,
                        )
                        # Format toggle
                        with html.Div(
                            style="display:flex;flex-direction:column;gap:4px;"
                        ):
                            html.Span(
                                "Format",
                                style="font-size:0.75rem;color:#9ca3af;",
                            )
                            with v3.VBtnToggle(
                                v_model=("export_format",),
                                mandatory=True,
                                variant="outlined",
                                density="compact",
                                color="cyan",
                                style="width:100%;",
                            ):
                                for _fi in _FMT_ITEMS:
                                    v3.VBtn(
                                        _fi["title"],
                                        value=_fi["value"],
                                        style="flex:1;font-family:monospace;background:#000000;",
                                    )
                        # Optional filename
                        v3.VTextField(
                            v_model=("export_filename",),
                            label="Filename (optional, no extension)",
                            variant="outlined",
                            density="compact",
                            color="cyan",
                            bg_color="#0d0d1a",
                            hide_details=True,
                            placeholder="auto-generated if blank",
                            style="font-family:monospace;",
                        )
                        # Status
                        with html.Div(
                            v_show=("export_status",),
                            style=(
                                "font-size:0.72rem;font-family:monospace;"
                                "background:#0d0d1a;padding:8px 10px;"
                                "border:1px solid #374151;word-break:break-all;"
                                "color:#9ca3af;"
                            ),
                        ):
                            html.Span("{{ export_status }}")
                    with v3.VCardActions(
                        style="padding:8px 16px 14px;gap:8px;justify-content:flex-end;"
                    ):
                        v3.VBtn(
                            "Close",
                            variant="text",
                            color="#9ca3af",
                            click="export_dialog_show = false",
                        )
                        v3.VBtn(
                            "Export",
                            variant="outlined",
                            color="cyan",
                            prepend_icon="mdi-export",
                            loading=("export_busy",),
                            disabled=("export_busy",),
                            click=server.controller.do_export,
                            style="font-family:monospace;",
                        )

            with v3.VSheet(
                classes="sage-content",
                style=(
                    "position:fixed;"
                    "top:var(--v-layout-top,48px);"
                    "left:var(--v-layout-left,0px);"
                    "right:var(--v-layout-right,0px);"
                    "bottom:var(--v-layout-bottom,36px);"
                    "display:flex;flex-direction:row;"
                    "overflow:hidden;"
                ),
                rounded=False,
                elevation=0,
                color="#0a0a0f",
            ):
                # ── Launch Mode wizard overlay ─────────────────────────────
                with html.Div(
                    v_show=("wiz_active",),
                    style=(
                        "position:absolute;inset:0;z-index:50;"
                        "background:#0a0a1a;overflow:hidden;"
                    ),
                ):
                    build_wizard_ui(server, _wiz_ctrl)

                # Live PSO plot gallery (right-docked; shown during a run)
                build_pso_gallery(server)

                # Render window + loading overlay
                with v3.VSheet(
                    style="position:relative;flex:1;height:100%;display:flex;min-width:0;",
                    color="transparent",
                    rounded=False,
                    elevation=0,
                ):
                    view = VtkRemoteView(
                        scene.plotter.ren_win,
                        style="flex:1;height:100%;display:block;min-width:0;",
                        # Full quality at all times, no interactive level-of-
                        # detail. ViSAGE shows real data at full fidelity in every
                        # situation: all points, full render resolution, full JPEG
                        # quality — both while the camera moves and when it
                        # settles, so there is no interactive→still "flash".
                        # Reducing resolution/quality during motion was tried and
                        # rejected: the transition is too striking, and degrading
                        # the view is not acceptable for real data. Interactivity
                        # on large clouds is a rendering-horsepower problem (use a
                        # real GPU / native rendering), NOT something to solve by
                        # dropping quality or subsampling. Any thinning of the
                        # cloud must stay a deliberate, user-driven choice.
                        interactive_ratio=1.0,
                        # (95 on Windows, 100 elsewhere; VISAGE_JPEG_QUALITY sets it.)
                        interactive_quality=_JPEG_QUALITY,
                        still_quality=_JPEG_QUALITY,
                    )
                    server.controller.view_update = view.update

                    # Story Mode text / equation overlays (over the render).
                    build_story_overlays(server)

                    # The selection card (selected galaxy, group, actions).
                    selection_card(server)

                    # Components: labels on the HD galaxy's parts, drawn by
                    # visage.js from hd_labels (an SVG over the view).
                    html.Div(
                        id="hd-labels",
                        style=(
                            "position:absolute;inset:0;z-index:5;"
                            "pointer-events:none;"
                        ),
                    )

                    # HD build message: what is being built and how far.
                    with html.Div(
                        v_show=("hd_building",),
                        style=(
                            "position:absolute;inset:0;z-index:7;"
                            "display:flex;align-items:center;"
                            "justify-content:center;pointer-events:none;"
                        ),
                    ):
                        with v3.VCard(
                            style=(
                                "padding:22px 30px;min-width:420px;"
                                "background:rgba(10,10,26,0.88);"
                                "border:1px solid #06b6d4;"
                                "pointer-events:auto;"
                            ),
                        ):
                            html.Div(
                                "{{ hd_progress }}",
                                style=(
                                    "font-size:1.5rem;font-weight:600;"
                                    "color:#e2e8f0;margin-bottom:12px;"
                                ),
                            )
                            v3.VProgressLinear(
                                model_value=("hd_percent",),
                                indeterminate=("hd_percent < 0",),
                                color="cyan",
                                height=10,
                                rounded=True,
                            )
                            with html.Div(
                                style=(
                                    "display:flex;align-items:center;"
                                    "margin-top:8px;"
                                ),
                            ):
                                v3.VBtn(
                                    "Cancel",
                                    size="small",
                                    variant="outlined",
                                    color="#9ca3af",
                                    click="trigger('hd_cancel')",
                                    title="Stop building (Esc)",
                                )
                                html.Div(
                                    "{{ hd_percent >= 0 ? hd_percent + '%' : '' }}",
                                    style=(
                                        "flex:1;font-size:1.1rem;"
                                        "color:#67e8f9;text-align:right;"
                                    ),
                                )

                    # Spectrum build message (Photometry tab's Make
                    # spectrum): what is happening and how far.
                    with html.Div(
                        v_show=("spec_busy",),
                        style=(
                            "position:absolute;inset:0;z-index:7;"
                            "display:flex;align-items:center;"
                            "justify-content:center;pointer-events:none;"
                        ),
                    ):
                        with v3.VCard(
                            style=(
                                "padding:22px 30px;min-width:420px;"
                                "background:rgba(10,10,26,0.88);"
                                "border:1px solid #06b6d4;"
                                "pointer-events:auto;"
                            ),
                        ):
                            html.Div(
                                "{{ spec_title }}",
                                style=(
                                    "font-size:1.5rem;font-weight:600;"
                                    "color:#e2e8f0;margin-bottom:4px;"
                                ),
                            )
                            html.Div(
                                "{{ spec_progress }}",
                                style=(
                                    "font-size:0.95rem;color:#a3adbb;"
                                    "margin-bottom:12px;"
                                ),
                            )
                            v3.VProgressLinear(
                                model_value=("spec_percent",),
                                color="cyan",
                                height=10,
                                rounded=True,
                            )
                            with html.Div(
                                style=(
                                    "display:flex;align-items:center;"
                                    "margin-top:8px;"
                                ),
                            ):
                                v3.VBtn(
                                    "Cancel",
                                    size="small",
                                    variant="outlined",
                                    color="#9ca3af",
                                    click="trigger('spectrum_cancel')",
                                    title="Stop making the spectrum",
                                )
                                html.Div(
                                    "{{ spec_percent }}%",
                                    style=(
                                        "flex:1;font-size:1.1rem;"
                                        "color:#67e8f9;text-align:right;"
                                    ),
                                )

                    # Loading card: every snapshot loading into memory (or
                    # playback frames being prepared), over a dimmed view.
                    with html.Div(
                        v_show=("load_card_show",),
                        style=(
                            "position:absolute;inset:0;z-index:8;"
                            "display:flex;align-items:center;"
                            "justify-content:center;"
                            "background:rgba(0,0,0,0.35);"
                        ),
                    ):
                        with v3.VCard(
                            style=(
                                "padding:22px 30px;min-width:420px;"
                                "background:rgba(10,10,26,0.88);"
                                "border:1px solid #06b6d4;"
                            ),
                        ):
                            html.Div(
                                "{{ load_card_text }}",
                                style=(
                                    "font-size:1.5rem;font-weight:600;"
                                    "color:#e2e8f0;margin-bottom:4px;"
                                ),
                            )
                            html.Div(
                                "Loaded once, so moving through time is "
                                "instant afterwards.",
                                v_show=("load_card_skip",),
                                style=(
                                    "font-size:0.8rem;color:#9ca3af;"
                                    "margin-bottom:10px;"
                                ),
                            )
                            v3.VProgressLinear(
                                model_value=("load_card_percent",),
                                color="cyan",
                                height=10,
                                rounded=True,
                            )
                            with html.Div(
                                style=(
                                    "display:flex;align-items:center;"
                                    "margin-top:8px;"
                                ),
                            ):
                                v3.VBtn(
                                    "Continue now",
                                    v_show=("load_card_skip",),
                                    size="small",
                                    variant="outlined",
                                    color="#9ca3af",
                                    click=server.controller.load_card_continue,
                                    title="Use ViSAGE now; the rest loads "
                                    "in the background",
                                )
                                html.Div(
                                    "{{ load_card_percent }}%",
                                    style=(
                                        "flex:1;font-size:1.1rem;"
                                        "color:#67e8f9;text-align:right;"
                                    ),
                                )

                    # Pre-rendered playback overlay — shown only during
                    # playback, where it flips through the cached frames. While
                    # frames are being rendered the live view stays put (the
                    # progress shows in the toolbar chip instead).
                    # Opacity-driven (not v_show) so dropping it CROSSFADES to
                    # the live view underneath — the server always renders the
                    # staged scene before clearing playback_active, so the fade
                    # lands on matching content instead of a hard cut.
                    with html.Div(
                        style=(
                            "'position:absolute;inset:0;z-index:6;"
                            "background:#000;display:flex;"
                            "align-items:center;justify-content:center;"
                            "transition:opacity 0.35s ease;' + "
                            "(playback_active ? 'opacity:1;'"
                            " : 'opacity:0;pointer-events:none;')",
                        ),
                    ):
                        html.Img(
                            src=("playback_frame",),
                            style=(
                                "width:100%;height:100%;object-fit:contain;"
                                "display:block;"
                            ),
                        )

                    # Pop-out console — floats over the viewport, mirrors
                    # the active session's history. Toggle from the
                    # Console tab's "Pop-out" button. Drag-free for now;
                    # pinned to the bottom-left corner of the viewport.
                    with v3.VCard(
                        v_show=("console_popout_show",),
                        classes="sage-popout",
                        style=(
                            "position:absolute;left:24px;bottom:24px;"
                            "width:560px;min-width:240px;"
                            "height:360px;min-height:160px;"
                            "background:rgba(13,13,26,0.92);"
                            "border:1px solid #06b6d4;"
                            "box-shadow:0 0 18px rgba(6,182,212,0.30);"
                            "display:flex;flex-direction:column;"
                            "z-index:10;color:#e2e8f0;"
                            "resize:both;overflow:hidden;"
                        ),
                        elevation=0,
                        rounded=False,
                    ):
                        # Title bar — also the drag handle (cursor:move +
                        # ".sage-popout-handle" picked up by the global
                        # drag script).
                        with html.Div(
                            classes="sage-popout-handle",
                            style=(
                                "display:flex;align-items:center;"
                                "padding:4px 8px;gap:8px;"
                                "border-bottom:1px solid #1f2937;"
                                "flex-shrink:0;cursor:move;"
                                "user-select:none;"
                            ),
                        ):
                            html.Span(
                                "CONSOLE  (Console {{ console_active_id }})",
                                style=(
                                    "font-size:0.75rem;font-weight:700;"
                                    "letter-spacing:0.08em;color:#06b6d4;"
                                ),
                            )
                            v3.VSpacer()
                            v3.VBtn(
                                icon="mdi-fullscreen",
                                size="x-small",
                                variant="text",
                                color="#9ca3af",
                                classes="sage-popout-max-btn",
                            )
                            v3.VBtn(
                                icon="mdi-close",
                                size="x-small",
                                variant="text",
                                color="#9ca3af",
                                click=server.controller.console_toggle_popout,
                            )
                        # Terminal mode: xterm.js instance in the pop-out.
                        html.Div(
                            id=("'sage-pty-popout-' + console_active_id",),
                            v_show=("console_mode === 'terminal'",),
                            style=("flex:1 1 0;min-height:0;overflow:hidden;"),
                        )
                        # Command mode: mirrored history + input.
                        with v3.VSheet(
                            color="#0a0a0f",
                            classes="sage-console-scroll",
                            v_show=("console_mode === 'command'",),
                            style=(
                                "flex:1 1 0;min-height:0;overflow-y:auto;"
                                "padding:6px 10px;font-family:monospace;"
                                "font-size:0.72rem;line-height:1.4;"
                            ),
                        ):
                            with html.Div(
                                v_for=("entry in console_history",),
                                key=("'po-' + entry.id",),
                                style=(
                                    "padding:2px 0 4px;"
                                    "border-bottom:1px solid #1f2937;"
                                ),
                            ):
                                html.Div(
                                    "{{ entry.cmd }}",
                                    style=(
                                        "color:cyan;white-space:pre-wrap;"
                                        "font-family:monospace;"
                                    ),
                                )
                                html.Div(
                                    "{{ entry.out }}",
                                    style=(
                                        "color:#9ca3af;white-space:pre-wrap;"
                                    ),
                                )
                        with html.Div(
                            v_show=("console_mode === 'command'",),
                            style="padding:6px 8px;flex-shrink:0;",
                        ):
                            with html.Div(
                                raw_attrs=[
                                    'data-enter-click="btn-console-run"'
                                ],
                            ):
                                v3.VTextField(
                                    v_model=("console_input",),
                                    label="SAGE Commands",
                                    hide_details=True,
                                    variant="outlined",
                                    bg_color="#1a1a2e",
                                    density="compact",
                                    style="font-family:monospace;",
                                    keydown_enter=server.controller.console_submit,
                                )

                    # Snackbar for compatibility / status messages
                    v3.VSnackbar(
                        "{{ notice_text }}",
                        v_model=("notice_show",),
                        timeout=("notice_timeout",),
                        color=("notice_color",),
                        location="top",
                        contained=True,
                        style="margin-top:16px;",
                    )

                    # Library media pop-outs — one draggable card per open item,
                    # floating over the viewport just like the other sage-popout cards.
                    with v3.VCard(
                        v_for=("item in library_items",),
                        key=("item.id",),
                        classes="sage-popout",
                        style=(
                            "`position:absolute;"
                            "top:${item.top_px}px;"
                            "right:${item.right_px}px;"
                            "width:480px;min-width:240px;"
                            "height:380px;min-height:200px;"
                            "display:flex;flex-direction:column;"
                            "resize:both;overflow:hidden;"
                            "background:rgba(17,24,39,0.92);"
                            "backdrop-filter:blur(6px);"
                            "border:1px solid #374151;"
                            "color:#e2e8f0;"
                            # In Story Mode, lift library cards above the story
                            # overlay layer (z40) + HUD (z50) so opened items sit
                            # on top of the scene; normal-use layering unchanged.
                            "z-index:${story_active ? 60 : 5};`",
                        ),
                        elevation=8,
                    ):
                        with v3.VCardTitle(
                            classes="sage-popout-handle",
                            style=(
                                "display:flex;align-items:center;"
                                "font-size:0.85rem;letter-spacing:0.06em;"
                                "padding:10px 12px 6px;cursor:move;"
                                "flex-shrink:0;"
                            ),
                        ):
                            v3.VIcon(
                                "mdi-folder-multimedia-outline",
                                size="small",
                                color="cyan",
                                style="margin-right:6px;",
                            )
                            html.Span("{{ item.name }}")
                            v3.VSpacer()
                            v3.VBtn(
                                icon="mdi-fullscreen",
                                size="x-small",
                                variant="text",
                                classes="sage-popout-max-btn",
                            )
                            v3.VBtn(
                                icon="mdi-close",
                                size="x-small",
                                variant="text",
                                click=(
                                    server.controller.library_close_item,
                                    "[item.id]",
                                ),
                            )
                        v3.VDivider()
                        with v3.VCardText(
                            style=(
                                "padding:8px 12px;flex:1 1 0;min-height:0;"
                                "overflow:auto;display:flex;"
                                "align-items:center;justify-content:center;"
                            ),
                        ):
                            html.Img(
                                src=("item.data_url",),
                                v_if=("item.kind === 'image'",),
                                style=(
                                    "max-width:100%;max-height:100%;"
                                    "width:auto;height:auto;"
                                    "object-fit:contain;display:block;"
                                ),
                            )
                            html.Video(
                                src=("item.data_url",),
                                v_if=("item.kind === 'video'",),
                                controls=True,
                                autoplay=True,
                                muted=True,
                                loop=True,
                                style=(
                                    "max-width:100%;max-height:100%;"
                                    "width:auto;height:auto;"
                                    "object-fit:contain;display:block;"
                                ),
                            )

                    # Loading overlay shown while a model loads
                    with v3.VOverlay(
                        v_model=("model_loading",),
                        contained=True,
                        persistent=True,
                        classes="d-flex align-center justify-center",
                        scrim="rgba(0,0,0,0.78)",
                    ):
                        with v3.VSheet(
                            color="#1a1a2e",
                            style=(
                                "padding:40px 56px;border-radius:6px;"
                                "border:1px solid #06b6d4;"
                                "display:flex;flex-direction:column;align-items:center;"
                                "gap:18px;min-width:480px;max-width:640px;"
                                "box-shadow:0 0 24px rgba(6,182,212,0.35);"
                            ),
                        ):
                            v3.VLabel(
                                "SWITCHING MODELS, PLEASE HOLD…",
                                style=(
                                    "font-size:1.35rem;font-weight:700;"
                                    "letter-spacing:0.12em;color:#06b6d4;"
                                    "text-align:center;"
                                ),
                            )
                            v3.VProgressLinear(
                                indeterminate=True,
                                color="cyan",
                                height=4,
                                style="width:100%;",
                            )
                            v3.VLabel(
                                "{{ model_quip }}",
                                style=(
                                    "font-size:0.95rem;color:#e2e8f0;"
                                    "text-align:center;font-style:italic;"
                                    "min-height:1.4em;"
                                ),
                            )

                    # Box strip — shown when at least one adjacent box is loaded
                    with html.Div(
                        v_show=(
                            "box_strip_items && box_strip_items.length > 1",
                        ),
                        style=(
                            "position:absolute;bottom:0;left:0;right:0;"
                            "z-index:10;height:44px;"
                            "background:rgba(0,0,0,0.85);"
                            "display:flex;align-items:center;gap:6px;padding:0 10px;"
                            "border-top:1px solid rgba(255,255,255,0.08);"
                        ),
                    ):
                        with html.Div(
                            v_for=("b in box_strip_items",),
                            key=("'bs-' + b.name",),
                            click=(
                                server.controller.set_active_box,
                                "[b.name]",
                            ),
                            style=(
                                "b.active "
                                "? 'display:flex;align-items:center;gap:6px;"
                                "padding:4px 10px;border-radius:4px;cursor:pointer;"
                                "border:1px solid #00ffff;"
                                "background:rgba(0,255,255,0.08)' "
                                ": 'display:flex;align-items:center;gap:6px;"
                                "padding:4px 10px;border-radius:4px;cursor:pointer;"
                                "border:1px solid #374151'",
                            ),
                        ):
                            html.Span(
                                "{{ b.label }}",
                                style=(
                                    "b.active "
                                    "? 'color:#00ffff;font-weight:700;"
                                    "font-size:0.7rem;font-family:monospace' "
                                    ": 'color:#9ca3af;"
                                    "font-size:0.7rem;font-family:monospace'",
                                ),
                            )
                            v3.VBtn(
                                "CLR",
                                v_if=("!b.primary",),
                                size="x-small",
                                variant="text",
                                style="font-size:0.6rem;min-width:28px;color:#6b7280;",
                                click=(
                                    server.controller.clear_box,
                                    "[b.name]",
                                ),
                            )

                # Right panel — layers + navigation tabs. Fixed 300px flex
                # child. v_show lets a Story Mode scene hide it
                # (chrome.hide_panel); there is no general-use toggle.
                with v3.VSheet(
                    v_show=("!panels_hidden",),
                    style=(
                        "width:300px;min-width:300px;max-width:300px;"
                        "flex:0 0 300px;"
                        "height:100%;box-sizing:border-box;overflow:hidden;"
                    ),
                    color="#000000",
                    rounded=False,
                    elevation=0,
                ):
                    build_navigation_panel(server, scene)
                    build_selection(server, scene)

        with layout.footer as footer:
            footer.height = 0
            footer.style = "display:none;"
            build_info_panel(server, scene)

    return server, scene

Example: headless snapshot render

from visage.io import parse_par, SnapshotTable
from visage.io.halo_reader import load_halo_snapshot
from visage.io.galaxy_reader import load_galaxy_snapshot
import pyvista as pv

cfg = parse_par("input/millennium.par")
snap_table = SnapshotTable(cfg.snap_list_path)

halos = load_halo_snapshot(cfg.tree_dir, cfg.tree_name, snap_num=63)
galaxies = load_galaxy_snapshot(cfg.hdf5_path, snap_num=63)

pl = pv.Plotter(off_screen=True)
pl.add_points(halos.positions, color="cyan", opacity=0.02, point_size=3)
pl.add_points(galaxies.positions, color="white", point_size=2)
pl.screenshot("snapshot_63.png")