VegasAfterglow 2.0.6

High-performance C++ framework with Python interface for modeling Gamma-Ray Burst (GRB) afterglows.

VegasAfterglow

Latest Release Notes | Interactive Web Tool | Validation Report | Install Now

VegasAfterglow is a high-performance C++ framework with a user-friendly Python interface designed for the comprehensive modeling of Gamma-Ray Burst (GRB) afterglows. It achieves exceptional computational efficiency, enabling the generation of multi-wavelength light curves in milliseconds and facilitating robust Markov Chain Monte Carlo (MCMC) parameter inference. The framework incorporates advanced models for shock dynamics (both forward and reverse shocks), diverse radiation mechanisms (synchrotron with self-absorption, and inverse Compton scattering with Klein-Nishina corrections), and complex structured jet configurations.

Ecosystem & Integration

• Advanced Inference: While VegasAfterglow includes a standalone fitter, its models are also integrated into Redback for more advanced Bayesian inference, model comparison, and multi-messenger analysis.
• Community Tools: We also recommend checking out PyFRS, Jetsimpy, ASGARD as tool for afterglow modeling.

---

Quick Start

pip install VegasAfterglow

import numpy as np
import matplotlib.pyplot as plt
from VegasAfterglow import ISM, TophatJet, Observer, Radiation, Model

model = Model(
    jet=TophatJet(theta_c=0.1, E_iso=1e52, Gamma0=300),
    medium=ISM(n_ism=1),
    observer=Observer(lumi_dist=1e26, z=0.1, theta_obs=0),
    fwd_rad=Radiation(eps_e=0.1, eps_B=1e-3, p=2.3),
)

times = np.logspace(2, 8, 100)            # 100 s to ~3 years
bands = [1e9, 1e14, 1e17]                 # radio, near-IR, X-ray
flux = model.flux_density_grid(times, bands).total

for f, nu in zip(flux, bands):
    plt.loglog(times, f, label=f"{nu:.0e} Hz")
plt.xlabel("Time (s)"); plt.ylabel("Flux density (erg/cm²/s/Hz)"); plt.legend()

Prefer no install? Open the interactive web tool — parameter sliders, observational data upload, mobile-friendly.

---

Table of Contents

• VegasAfterglow
  • Quick Start
  • Table of Contents
  • Features
  • Installation
    • Python Installation
    • C++ Installation
  • Usage
    • Quick Start
      • Command-Line Tool (vegasgen)
      • Python API
    • Light Curve & Spectrum Calculation
    • Internal Quantities Evolution
    • Sky Image
    • MCMC Parameter Fitting
  • Performance Highlights
  • Architecture
  • Validation & Testing
  • Documentation
  • Contributing
  • License
  • Acknowledgments & Citation

---

Features

Shock Dynamics

• Forward and Reverse Shock Modeling: Simulates both shocks via shock crossing dynamics with arbitrary magnetization levels and shell thicknesses.
• Relativistic and Non-Relativistic Regimes: Accurately models shock evolution across all velocity regimes.
• Adiabatic and Radiative Blast Waves: Supports smooth transition between adiabatic and radiative blast waves.
• Ambient Medium: Supports uniform Interstellar Medium (ISM), stellar wind environments, and user-defined density profiles.
• Energy and Mass Injection: Supports user-defined profiles for continuous energy and/or mass injection into the blast wave.

Jet Structure & Geometry

• Structured Jet Profiles: Allows user-defined angular profiles for energy distribution, initial Lorentz factor, and magnetization.
• Arbitrary Viewing Angles: Supports off-axis observers at any viewing angle relative to the jet axis.
• Jet Spreading: Includes lateral expansion dynamics for realistic jet evolution (experimental).
• Non-Axisymmetric Jets: Capable of modeling complex, non-axisymmetric jet structures.

Radiation Mechanisms

• Synchrotron Radiation: Calculates synchrotron emission from shocked electrons.
• Synchrotron Self-Absorption (SSA): Includes SSA effects, crucial at low frequencies.
• Inverse Compton (IC) Scattering: Models IC processes, including:
  • Synchrotron Self-Compton (SSC) from both forward and reverse shocks.
  • Pairwise IC between forward and reverse shock electron and photon populations (experimental).
  • Includes Klein-Nishina corrections for accurate synchrotron and IC emission.

---

Installation

VegasAfterglow is available as a Python package with C++ source code also provided for direct use.

Python Installation

To install VegasAfterglow using pip:

pip install VegasAfterglow

This installs the core physics engine. To also install MCMC fitting support:

pip install VegasAfterglow[mcmc]

Performance note: PyPI wheels are compiled for broad CPU compatibility and run on any machine. For the fastest runs on your own hardware, build from source (see Install from Source below) — the build auto-tunes to your CPU's instruction set.

The interactive web tool is deployed across multiple global regions for low latency. If none of the servers provide a satisfactory experience from your location, you can run the interactive web tool locally for the best responsiveness (FastAPI + Next.js):

# terminal 1: backend
cd webtool/backend
python3 -m venv .venv
source .venv/bin/activate
pip install -e '../..[webtool]'
uvicorn app.main:app --reload --host 0.0.0.0 --port 8000

# terminal 2: frontend
cd webtool/frontend
cp .env.local.example .env.local
npm install
npm run dev

Open http://localhost:3000. Detailed local/deploy guide: webtool/README.md and webtool/DEPLOY.md.

The hosted Interactive Web Tool provides real-time parameter-tunable afterglow modeling, mobile-phone compatible layout, and observational data input/upload support.

VegasAfterglow requires Python 3.9 or higher.

Alternative: Install from Source (click to expand/collapse)

Build from source when pip installation is not viable, when the development version is required, or for maximum performance — a source build auto-tunes to your CPU's instruction set, which is typically several times faster than the portable PyPI wheel.

The quickest route compiles only VegasAfterglow from source while keeping binary dependencies:

pip install --no-binary VegasAfterglow VegasAfterglow

Or clone and install:

1. Clone this repository:

git clone https://github.com/YihanWangAstro/VegasAfterglow.git

2. Navigate to the directory and install the Python package:

cd VegasAfterglow
pip install .

Standard development environments typically include the necessary prerequisites (C++20 compatible compiler). For build-related issues, refer to the prerequisites section in C++ Installation.

C++ Installation

For advanced users who need to compile and use the C++ library directly:

Instructions for C++ Installation (click to expand/collapse)

1. Clone the repository (if not previously done):

git clone https://github.com/YihanWangAstro/VegasAfterglow.git
cd VegasAfterglow

2. Compile and run tests:

make tests

Upon successful compilation, you can create custom C++ problem generators using the VegasAfterglow interfaces. For implementation details, refer to the Creating Custom Problem Generators with C++ section or examine the example problem generators in tests/demo/.

Build Prerequisites (click to expand for dependency information)

The following development tools are required:

• C++20 compatible compiler:

  • Linux: GCC 10+ or Clang 13+
  • macOS: Apple Clang 13+ (with Xcode 13+) or GCC 10+ (via Homebrew)
  • Windows: MSVC 19.29+ (Visual Studio 2019 16.10+) or MinGW-w64 with GCC 10+

• Build tools:

  • Make (GNU Make 4.0+ recommended)

---

Usage

Quick Start

VegasAfterglow can be used in three ways:

• Interactive web tool — real-time, parameter-tunable afterglow modeling for light curves, spectra, and sky images, with mobile-phone compatible UI and observational data input/upload support (CSV/TXT/XLS/XLSX). No installation required.
• Command line — generate light curves instantly with the vegasgen command, no code needed.
• Python API — full programmatic control for custom analysis, MCMC fitting, and access to internal quantities.

Command-Line Tool (vegasgen)

# Default tophat jet, ISM, on-axis — just run it
vegasgen

# Gaussian jet at z=0.5 with plot output
vegasgen --jet gaussian --z 0.5 --plot

# Use filter names for frequencies, save to file
vegasgen --nu R J 1e18 -o lightcurve.csv

# Wind medium with SSC, save plot as PNG
vegasgen --medium wind --ssc --plot -o lightcurve.png

All parameters have sensible defaults. Run vegasgen --help for the full option list, or see the CLI documentation for detailed usage.

Python API

from VegasAfterglow import ISM, TophatJet, Observer, Radiation, Model
import numpy as np

model = Model(TophatJet(0.1, 1e52, 300), ISM(1), Observer(1e26, 0.1, 0), Radiation(0.1, 1e-3, 2.3))
result = model.flux_density_grid(np.logspace(2, 8, 100), [1e9, 1e14, 1e17])

See the detailed sections below and the example notebooks (script/quick.ipynb, script/details.ipynb, script/vegas-mcmc.ipynb) for Python usage including light curves, spectra, internal quantities, and MCMC fitting.

Light Curve & Spectrum Calculation

The example below walks through the main components needed to model a GRB afterglow, from setting up the physical parameters to producing light curves and spectra via script/quick.ipynb.

Model Setup (click to expand/collapse)

import numpy as np
import matplotlib.pyplot as plt
from VegasAfterglow import ISM, TophatJet, Observer, Radiation, Model

medium = ISM(n_ism=1)                                    # constant density ISM [cm^-3]
jet = TophatJet(theta_c=0.1, E_iso=1e52, Gamma0=300)    # top-hat jet [rad, erg, dimensionless]
obs = Observer(lumi_dist=1e26, z=0.1, theta_obs=0)      # observer [cm, dimensionless, rad]
rad = Radiation(eps_e=1e-1, eps_B=1e-3, p=2.3)          # microphysics

model = Model(jet=jet, medium=medium, observer=obs, fwd_rad=rad)

Light Curve Calculation (click to expand/collapse)

Compute and plot multi-wavelength light curves:

times = np.logspace(2, 8, 100)
bands = np.array([1e9, 1e14, 1e17])  # radio, near-IR, X-ray

# times must be in ascending order; frequencies can be in any order
results = model.flux_density_grid(times, bands)

plt.figure(figsize=(4.8, 3.6), dpi=200)
for i, nu in enumerate(bands):
    exp = int(np.floor(np.log10(nu)))
    base = nu / 10**exp
    plt.loglog(times, results.total[i,:], label=fr'${base:.1f} \times 10^{{{exp}}}$ Hz')

plt.xlabel('Time (s)')
plt.ylabel('Flux Density (erg/cm²/s/Hz)')
plt.legend()
plt.title('Light Curves')

Spectrum Analysis (click to expand/collapse)

Broadband spectra at different epochs:

frequencies = np.logspace(5, 22, 100)
epochs = np.array([1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8])

# times must be in ascending order
results = model.flux_density_grid(epochs, frequencies)

plt.figure(figsize=(4.8, 3.6), dpi=200)
colors = plt.cm.viridis(np.linspace(0, 1, len(epochs)))
for i, t in enumerate(epochs):
    exp = int(np.floor(np.log10(t)))
    base = t / 10**exp
    plt.loglog(frequencies, results.total[:,i], color=colors[i], label=fr'${base:.1f} \times 10^{{{exp}}}$ s')

plt.xlabel('Frequency (Hz)')
plt.ylabel('Flux Density (erg/cm²/s/Hz)')
plt.legend(ncol=2)
plt.title('Synchrotron Spectra')

Time-Frequency Pairs Calculation (click to expand/collapse)

For paired time-frequency points (t_i, nu_i) instead of a grid:

times = np.logspace(2, 8, 200)
frequencies = np.logspace(9, 17, 200)  # same length as times

results = model.flux_density(times, frequencies)  # times must be ascending

Three flux methods:

• flux_density_grid(times, freqs): NxM grid output
• flux_density(times, freqs): paired points (N output from N pairs)
• flux_density_exposures(times, freqs, expo): paired with exposure-time averaging

All return a FluxDict with .total, .fwd (.sync, .ssc), and .rvs (.sync, .ssc).

Internal Quantities Evolution

The example below walks through how you can check the evolution of internal quantities under various reference frames via script/details.ipynb.

Model Setup (click to expand/collapse)

import numpy as np
import matplotlib.pyplot as plt
from VegasAfterglow import ISM, TophatJet, Observer, Radiation, Model

z = 0.1
medium = ISM(n_ism=1)
jet = TophatJet(theta_c=0.3, E_iso=1e52, Gamma0=100)
obs = Observer(lumi_dist=1e26, z=z, theta_obs=0.)
rad = Radiation(eps_e=1e-1, eps_B=1e-3, p=2.3)

model = Model(jet=jet, medium=medium, observer=obs, fwd_rad=rad, resolutions=(0.06, 0.15, 6))

Get the simulation quantities (click to expand/collapse)

details = model.details(t_min=1e0, t_max=1e8)

Returns a SimulationDetails object. All quantities are 3D numpy arrays with shape (n_phi, n_theta, n_time), where axis 0 is the azimuthal angle, axis 1 is the polar angle, and axis 2 is the time step. For example, details.fwd.Gamma[0, 0, :] gives the Lorentz factor time evolution along the jet symmetry axis. Access forward shock via details.fwd (reverse shock via details.rvs when enabled).

Category	Attributes
Grid	phi, theta, t_src (source frame), fwd.t_comv (comoving), fwd.t_obs (observer)
Dynamics	Gamma, Gamma_th, r, theta, Doppler
Microphysics	B_comv [G], N_p, N_e, I_nu_max [erg/cm²/s/Hz]
Electron Lorentz factors	gamma_a, gamma_m, gamma_c, gamma_M
Characteristic frequencies	nu_a, nu_m, nu_c, nu_M [Hz]
Spectra (callable)	sync_spectrum, ssc_spectrum, Y_spectrum

Frame: gamma_*, nu_*, I_nu_max, and B_comv are all in the comoving (rest) frame of the shocked fluid. Convert a comoving frequency to the observer frame with nu_obs = nu_comv * Doppler / (1 + z) (the inverse applies for observer → comoving). See the internal-quantities docs for details and the inline example below for usage.

Checking the evolution of various parameters (click to expand/collapse)

Plot shock quantities across source, comoving, and observer frames:

attrs =['Gamma', 'B_comv', 'N_p','r','N_e','I_nu_max']
ylabels = [r'$\Gamma$', r'$B^\prime$ [G]', r'$N_p$', r'$r$ [cm]', r'$N_e$', r'$I_{\nu, \rm max}^\prime$ [erg/s/Hz]']

frames = ['t_src', 't_comv', 't_obs']
titles = ['source frame', 'comoving frame', 'observer frame']
colors = ['C0', 'C1', 'C2']
xlabels = [r'$t_{\rm src}$ [s]', r'$t^\prime$ [s]', r'$t_{\rm obs}$ [s]']
plt.figure(figsize= (4.2*len(frames), 3*len(attrs)))

# All 3D arrays have shape (n_phi, n_theta, n_time).
# Using [0, 0, :] extracts the time evolution at the first azimuthal and
# polar bin, i.e. along the jet symmetry axis.
for i, frame in enumerate(frames):
    for j, attr in enumerate(attrs):
        plt.subplot(len(attrs), len(frames) , j * len(frames) + i + 1)
        if j == 0:
            plt.title(titles[i])
        value = getattr(details.fwd, attr)
        if frame == 't_src':
            t = getattr(details, frame)
        else:
            t = getattr(details.fwd, frame)
        plt.loglog(t[0, 0, :], value[0, 0, :], color='k',lw=2.5)
        plt.loglog(t[0, 0, :], value[0, 0, :], color=colors[i])

        plt.xlabel(xlabels[i])
        plt.ylabel(ylabels[j])

plt.tight_layout()
plt.savefig('shock_quantities.png', dpi=300,bbox_inches='tight')

See script/details.ipynb for additional plots including characteristic electron energies, synchrotron frequencies, Doppler factor maps, and equal arrival time surfaces.

Per-Cell Spectrum Evaluation (click to expand/collapse)

In addition to scalar quantities, details() provides callable spectrum accessors that let you evaluate the comoving-frame synchrotron, SSC, and Compton-Y spectra at arbitrary frequencies for each grid cell. To use SSC and Y spectrum, enable SSC in the radiation model:

rad = Radiation(eps_e=1e-1, eps_B=1e-3, p=2.3, ssc=True)
model = Model(jet=jet, medium=medium, observer=obs, fwd_rad=rad, resolutions=(0.06, 0.15, 6))
details = model.details(t_min=1e0, t_max=1e8)

nu_comv = np.logspace(8, 20, 200)  # comoving frame frequency [Hz]

# Synchrotron spectrum at cell (phi=0, theta=0, t=5)
I_syn = details.fwd.sync_spectrum[0, 0, 5](nu_comv)   # erg/s/Hz/cm²/sr

# SSC spectrum at the same cell (requires ssc=True)
I_ssc = details.fwd.ssc_spectrum[0, 0, 5](nu_comv)    # erg/s/Hz/cm²/sr

# Compton-Y parameter as a function of electron Lorentz factor
gamma = np.logspace(1, 8, 200)
Y = details.fwd.Y_spectrum[0, 0, 5](gamma)            # dimensionless

These callable accessors are also available on details.rvs when a reverse shock is configured. The sync_spectrum and Y_spectrum are always available; ssc_spectrum is None unless ssc=True.

Sky Image

Generate spatially resolved images of the afterglow at any observer time and frequency. See script/sky-image.ipynb for the full example.

Single Frame (click to expand/collapse)

from matplotlib.colors import LogNorm
from VegasAfterglow.units import uas

img = model.sky_image([1e6], nu_obs=1e9, fov=500 * uas, npixel=256)

fig, ax = plt.subplots(dpi=100)
im = ax.imshow(
    img.image[0].T,
    origin="lower",
    extent=img.extent / uas,
    cmap="inferno",
    norm=LogNorm(),
)
ax.set_xlabel(r"$\Delta x$ ($\mu$as)")
ax.set_ylabel(r"$\Delta y$ ($\mu$as)")
fig.colorbar(im, label=r"Surface brightness (erg/cm$^2$/s/Hz/sr)")

The return value is a SkyImage object with:

• img.image: 3D array (n_frames, npixel, npixel) — surface brightness in erg/cm²/s/Hz/sr
• img.extent: angular extent in radians (pass to imshow(extent=...))
• img.pixel_solid_angle: pixel solid angle in sr

Multi-Frame Movie (click to expand/collapse)

Pass an array of observer times to generate an image sequence efficiently — dynamics are solved once, each frame only re-renders the sky projection:

times = np.logspace(4, 8, 60)
imgs = model.sky_image(times, nu_obs=1e9, fov=2000 * uas, npixel=128)
# imgs.image.shape == (60, 128, 128)

Off-Axis Observer (click to expand/collapse)

For off-axis observers, the image centroid drifts across the sky (superluminal apparent motion):

model_offaxis = Model(
    jet=TophatJet(theta_c=0.1, E_iso=1e52, Gamma0=200),
    medium=ISM(n_ism=1),
    observer=Observer(lumi_dist=1e26, z=0.1, theta_obs=0.4),
    fwd_rad=Radiation(eps_e=1e-1, eps_B=1e-3, p=2.3),
)

times_oa = np.logspace(5, 8, 30)
imgs_oa = model_offaxis.sky_image(times_oa, nu_obs=1e9, fov=5000 * uas, npixel=128)

MCMC Parameter Fitting

We provide some example data files in the data folder. Remember to keep your copy in the same directory as the original to ensure all data paths work correctly.

1. Configuring the Model (click to expand/collapse)

# Requires: pip install VegasAfterglow[mcmc]
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import corner
from VegasAfterglow import Fitter, ParamDef, Scale

# Create the fitter with model configuration and add data
fitter = Fitter(
    z=1.58,
    lumi_dist=3.364e28,
    jet="powerlaw",
    medium="wind",
)

All model configuration is passed directly to the Fitter constructor as keyword arguments. Check the documentation for all available options.

2. Defining Parameters and Running MCMC (click to expand/collapse)

Define parameters with ParamDef(name, min, max, scale). Use Scale.log for parameters spanning orders of magnitude, Scale.linear for narrow ranges, and Scale.fixed to hold a parameter constant:

mc_params = [
    ParamDef("E_iso",      1e50,  1e54,  Scale.log),       # Isotropic energy [erg]
    ParamDef("Gamma0",        5,  1000,  Scale.log),       # Lorentz factor at the core
    ParamDef("theta_c",     0.0,   0.5,  Scale.linear),    # Core half-opening angle [rad]
    ParamDef("k_e",           2,     2,  Scale.fixed),     # Energy power law index
    ParamDef("k_g",           2,     2,  Scale.fixed),     # Lorentz factor power law index
    ParamDef("theta_v",     0.0,   0.0,  Scale.fixed),     # Viewing angle [rad]
    ParamDef("p",             2,     3,  Scale.linear),    # Shocked electron power law index
    ParamDef("eps_e",      1e-2,   0.5,  Scale.log),       # Electron energy fraction
    ParamDef("eps_B",      1e-4,   0.5,  Scale.log),       # Magnetic field energy fraction
    ParamDef("A_star",     1e-3,     1,  Scale.log),       # Wind parameter
    ParamDef("xi_e",       1e-3,     1,  Scale.log),       # Electron acceleration fraction
]

Add observational data and run the fit:

# Add light curves at specific frequencies (all quantities in CGS units)
t_data = [1e3, 2e3, 5e3, 1e4, 2e4]  # Time in seconds
flux_data = [1e-26, 8e-27, 5e-27, 3e-27, 2e-27]  # erg/cm²/s/Hz
flux_err = [1e-28, 8e-28, 5e-28, 3e-28, 2e-28]    # erg/cm²/s/Hz
fitter.add_flux_density(nu=4.84e14, t=t_data, f_nu=flux_data, err=flux_err)

# Load from CSV files
import pandas as pd
bands = [2.4e17, 4.84e14, 1.4e14]
lc_files = ["data/ep.csv", "data/r.csv", "data/vt-r.csv"]
for nu, fname in zip(bands, lc_files):
    df = pd.read_csv(fname)
    fitter.add_flux_density(nu=nu, t=df["t"], f_nu=df["Fv_obs"], err=df["Fv_err"])

# Add spectra at specific times
fitter.add_spectrum(t=3000, nu=nu_data, f_nu=spectrum_data, err=spectrum_err)

# Option 1: Nested sampling with dynesty (computes evidence, robust for multimodal posteriors)
result = fitter.fit(
    mc_params,
    resolution=(0.06, 0.15, 6),    # Grid resolution (phi, theta, t)
    sampler="dynesty",             # Nested sampling algorithm
    nlive=1000,                    # Number of live points
    walks=100,                     # Number of random walks per live point
    dlogz=0.5,                     # Stopping criterion (evidence tolerance)
    npool=8,                       # Number of parallel threads
    top_k=10,                      # Number of best-fit parameters to return
)

# Option 2: MCMC with emcee (faster, good for unimodal posteriors, not optimal for multimodal posteriors)
result = fitter.fit(
    mc_params,
    resolution=(0.06, 0.15, 6),    # Grid resolution (phi, theta, t)
    sampler="emcee",               # MCMC sampler
    nsteps=25000,                  # Number of steps per walker
    nburn=5000,                    # Burn-in steps to discard
    thin=1,                        # Save every nth sample
    npool=8,                       # Number of parallel threads
    top_k=10,                      # Number of best-fit parameters to return
)

The result object contains samples, labels, latex_labels, top_k_params, top_k_log_probs, and bilby_result. It also exposes summary(), save(), and load() for inspecting and persisting fits — see below.

3. Analyzing Results and Generating Predictions (click to expand/collapse)

Inspect the top-K best fits:

print(result.summary())

Generate predictions with the best-fit parameters:

best_params = result.top_k_params[0]

t_out = np.logspace(2, 9, 150)
nu_out = np.logspace(16, 20, 150)

lc = fitter.flux_density_grid(best_params, t_out, bands)
spec = fitter.flux_density_grid(best_params, times, nu_out)

Corner plot:

flat_chain = result.samples.reshape(-1, result.samples.shape[-1])
corner.corner(flat_chain, labels=result.latex_labels, quantiles=[0.16, 0.5, 0.84], show_titles=True)

Persist the fit to disk and reload later (bilby-native HDF5 / JSON — interoperable with bilby.read_in_result(path)):

# After fitting -- saves the full Fitter state (constructor args, observation
# data, parameter definitions) alongside the samples.
fitter.save("my_fit.h5")

# Later session: one classmethod call reloads the configured Fitter and the
# FitResult. No manual reconfiguration, no MCMC re-run.
from VegasAfterglow import Fitter
fitter = Fitter.load("my_fit.h5")
result = fitter.result
print(result.summary())
lc = fitter.flux_density_grid(result.top_k_params[0], t_out, bands)

See script/vegas-mcmc.ipynb for a complete working example with data loading, fitting, and visualization.

---

Performance Highlights

VegasAfterglow is designed for speed. A tophat synchrotron light curve completes in under 1 ms; even the most demanding case — a structured jet with full SSC and Klein-Nishina cooling off-axis — finishes in ~140 ms. This makes MCMC parameter estimation lightning fast on a laptop:

• Tophat jet: 10,000 MCMC steps, 8 parameters, 15 data points for ~0.3 M samples ~15 s (Apple M2, 8 cores).
• Structured jet: same setup ~1 minute.

The charts below benchmark single-core wall-clock time by stage across four jet profiles (Tophat, Gaussian, Power-law, Two-component), two media (ISM/Wind), and $\theta_v/\theta_c = 0, 1, 2, 4$, all at the shipped default resolutions — calibrated per emission component to hold the validation gates across jet families and viewing angles (reverse-shock runs automatically use a denser grid). Cost scales with two factors: off-axis observers require a full 3D equal-arrival-time grid vs. 2D on-axis, and structured jets compute radiation independently at each $\theta$ vs. reusing one solution for tophat. Reverse shock adds a second shocked shell, roughly doubling the dynamics cost.

---

Architecture

---

Validation & Testing

VegasAfterglow ships a unified test framework (see TESTING.md) with three tiers under tests/:

• Code correctness — C++ unit tests (Boost.Test) and Python API tests (pytest), run on every push across Linux/macOS/Windows.
• Physics correctness — closure relations against standard afterglow theory, exact invariants (distance/redshift scalings, component additivity), golden-baseline regression, and parameter-space corner sweeps.
• Full validation (tests/validation/) — benchmark tests (resolution convergence, performance timing) and regression tests (shock dynamics scaling laws, characteristic frequency evolution, spectral power-law indices across Blandford-McKee and Sedov-Taylor phases). The suite's pinned fiducial resolution (0.1, 0.25, 10) achieves < 5% mean error for most configurations; the shipped defaults are gated separately per emission component (tests/python/test_default_convergence.py).

make test               # everything: C++ + Python suites + full validation + HTML report (slow)
make test-quick         # C++ + Python suites only (seconds)
make test-validation    # the full validation suite only (slow)

Every run writes test-report.html — a single self-contained page with all test outcomes, physics validation figures, and performance results.

See the latest Validation Report for full results (refreshed on every release and by a weekly scheduled run); per-release snapshots are archived in the reports index. To run the full validation locally:

pip install -e ".[test]" --config-settings=cmake.define.AFTERGLOW_PROFILE=ON
python tests/validation/run_validation.py --all        # full suite (timing + convergence + regression)
python tests/validation/run_validation.py --all --fast  # fast mode: skip convergence scans

Powered by Claude Code

---

Documentation

Full documentation — installation, examples, MCMC fitting guide, parameter reference, Python/C++ API, and troubleshooting — is available at Documentation.

For a history of changes and new features, see the Changelog.

---

Contributing

If you encounter any issues, have questions about the code, or want to request new features:

1. GitHub Issues - The most straightforward and fastest way to get help:

   • Open an issue at Issues
   • You can report bugs, suggest features, or ask questions
   • This allows other users to see the problem/solution as well
   • Can be done anonymously if preferred

2. Pull Requests - If you've implemented a fix or feature:

   • Fork the repository
   • Create a branch for your changes
   • Submit a pull request with your changes

3. Email - For private questions or discussions:

   • Contact the maintainers directly via email
   • Yihan Wang: [e-mail]
   • Bing Zhang: [e-mail]
   • Connery Chen: [e-mail]

We value all contributions and aim to respond to issues promptly. All communications are extremely welcome!

---

License

VegasAfterglow is released under the BSD-3-Clause License.

---

Acknowledgments & Citation

We would like to thank the contributors who helped improve VegasAfterglow. Special thanks to Weihua Lei, Shaoyu Fu, Liang-Jun Chen, Iris Yin, Cuiyuan Dai, Binbin Zhang and Nikhil Sarin for their invaluable work as beta testers, providing feedback and helping with bug fixes during development. We also thank the broader community for their suggestions and support.

If you find VegasAfterglow useful, please consider giving us a star — it helps others discover the project!

If you use VegasAfterglow in your research, we would be grateful if you could cite:

Wang, Y., Chen, C., & Zhang, B. (2026). VegasAfterglow: A high-performance framework for gamma-ray burst afterglows. Journal of High Energy Astrophysics, 50, 100490. ADS

@ARTICLE{2026JHEAp..5000490W,
       author = {{Wang}, Yihan and {Chen}, Connery and {Zhang}, Bing},
        title = "{VegasAfterglow: A high-performance framework for gamma-ray burst afterglows}",
      journal = {Journal of High Energy Astrophysics},
     keywords = {Gamma-ray bursts, Shocks, Relativistic jets, Computational methods, Open source software, High Energy Astrophysical Phenomena},
         year = 2026,
        month = feb,
       volume = {50},
          eid = {100490},
        pages = {100490},
          doi = {10.1016/j.jheap.2025.100490},
archivePrefix = {arXiv},
       eprint = {2507.10829},
 primaryClass = {astro-ph.HE},
       adsurl = {https://ui.adsabs.harvard.edu/abs/2026JHEAp..5000490W},
      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}

Zhang, B. (2018). The Physics of Gamma-Ray Bursts. Cambridge University Press. ADS

@book{Zhang2018,
  author    = {Zhang, Bing},
  title     = {{The Physics of Gamma-Ray Bursts}},
  publisher = {Cambridge University Press},
  year      = {2018},
  doi       = {10.1017/9781139226530}
}

---