VegasAfterglow 2.0.0rc1

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

VegasAfterglow

Latest Release Notes | Validation Report (PDF) | 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.

---

Table of Contents

• VegasAfterglow
  • Table of Contents
  • Features
  • Performance Highlights
  • Installation
    • Python Installation
    • C++ Installation
  • Usage
    • Quick Start
    • Light Curve & Spectrum Calculation
    • Internal Quantities Evolution
    • MCMC Parameter Fitting
    • Parameter Fitting with Redback
  • 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.

---

Performance Highlights

VegasAfterglow delivers exceptional computational performance through deep optimization of its core algorithms. A 100-point single-frequency light curve (forward shock and synchrotron emission only) for a structured jet viewed off-axis can be generated in approximately 1 millisecond on a single core of an Apple M2 chip. This performance enables comprehensive Bayesian inference on standard laptop hardware in seconds to minutes:

• Tophat jet: On-axis forward synchrotron only, 10,000 MCMC steps (320,000 likelihood evaluations) with 8 parameters and 15 data points completes in ~15 seconds on an Apple M2 laptop (8 cores).
• Structured jet: The same MCMC run completes in ~1 minute.

This acceleration allows rapid iteration over different physical models and makes VegasAfterglow well suited for both detailed analyses of individual GRB events and large-scale population studies.

---

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]

VegasAfterglow requires Python 3.8 or higher.

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

For cases where pip installation is not viable or when the development version is required:

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

Get your first light curve in under 10 lines:

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

model = Model(TophatJet(0.1, 1e52, 300), ISM(1), Observer(1e26, 0.1, 0), Radiation(0.1, 1e-3, 2.3))
times, bands = np.logspace(2, 8, 100), [1e9, 1e14, 1e17]
results = model.flux_density_grid(times, bands)
[plt.loglog(times, results.total[i,:]) for i in range(len(bands))]
plt.xlabel('Time (s)'); plt.ylabel('Flux (erg/cm²/s/Hz)'); plt.show()

Read the detailed sections below for more information.

---

We provide basic example scripts (script/quick.ipynb, script/details.ipynb and script/vegas-mcmc.ipynb) that demonstrate how to set up and run afterglow simulations. This section shows how to calculate light curves and spectra for a simple GRB afterglow model without the need for observational data and perform MCMC parameter fitting with observational data. The notebook can be run using either Jupyter Notebook or VSCode with the Jupyter extension.

To avoid conflicts when updating the repository in the future, make a copy of the example notebook in the same directory and work with the copy instead of the original.

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)

First, let's set up the physical components of our afterglow model, including the environment, jet, observer, and radiation parameters:

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

# 1. Define the circumburst environment (constant density ISM)
medium = ISM(n_ism=1) #in cgs unit

# 2. Configure the jet structure (top-hat with opening angle, energy, and Lorentz factor)
jet = TophatJet(theta_c=0.1, E_iso=1e52, Gamma0=300) #in cgs unit

# 3. Set observer parameters (distance, redshift, viewing angle)
obs = Observer(lumi_dist=1e26, z=0.1, theta_obs=0) #in cgs unit

# 4. Define radiation microphysics parameters
rad = Radiation(eps_e=1e-1, eps_B=1e-3, p=2.3)

# 5. Combine all components into a complete afterglow model
model = Model(jet=jet, medium=medium, observer=obs, fwd_rad=rad)

Light Curve Calculation (click to expand/collapse)

Now, let's compute and plot multi-wavelength light curves to see how the afterglow evolves over time:

# 1. Create logarithmic time array from 10² to 10⁸ seconds (100s to ~3yrs)
times = np.logspace(2, 8, 100)

# 2. Define observing frequencies (radio, optical, X-ray bands in Hz)
bands = np.array([1e9, 1e14, 1e17])

# 3. Calculate the afterglow emission at each time and frequency
# NOTE: times array must be in ascending order, frequencies can be in random order
results = model.flux_density_grid(times, bands)

# 4. Visualize the multi-wavelength light curves
plt.figure(figsize=(4.8, 3.6),dpi=200)

# 5. Plot each frequency band
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')

def add_note(plt):
    plt.annotate('jet break',xy=(3e4, 1e-26), xytext=(3e3, 5e-28), arrowprops=dict(arrowstyle='->'))
    plt.annotate(r'$\nu_m=\nu_a$',xy=(8e5, 2e-25), xytext=(7.5e4, 5e-24), arrowprops=dict(arrowstyle='->'))
    plt.annotate(r'$\nu=\nu_a$',xy=(4e6, 4e-25), xytext=(7.5e5, 5e-24), arrowprops=dict(arrowstyle='->'))

add_note(plt)
plt.xlabel('Time (s)')
plt.ylabel('Flux Density (erg/cm²/s/Hz)')
plt.legend()
plt.title('Light Curves')
plt.savefig('quick-lc.png',dpi=300,bbox_inches='tight')

Running the light curve script will produce this figure showing the afterglow evolution across different frequencies.

Spectrum Analysis (click to expand/collapse)

We can also examine how the broadband spectrum evolves at different times after the burst:

# 1. Define broad frequency range (10⁵ to 10²² Hz)
frequencies = np.logspace(5, 22, 100)

# 2. Select specific time epochs for spectral snapshots
epochs = np.array([1e2, 1e3, 1e4, 1e5 ,1e6, 1e7, 1e8])

# 3. Calculate spectra at each epoch
# NOTE: epochs array must be in ascending order, frequencies can be in random order
results = model.flux_density_grid(epochs, frequencies)

# 4. Plot broadband spectra at each epoch
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')

# 5. Add vertical lines marking the bands from the light curve plot
for i, band in enumerate(bands):
    exp = int(np.floor(np.log10(band)))
    base = band / 10**exp
    plt.axvline(band,ls='--',color='C'+str(i))

plt.xlabel('frequency (Hz)')
plt.ylabel('flux density (erg/cm²/s/Hz)')
plt.legend(ncol=2)
plt.title('Synchrotron Spectra')
plt.savefig('quick-spec.png',dpi=300,bbox_inches='tight')

The spectral analysis code will generate this visualization showing spectra at different times, with vertical lines indicating the frequencies calculated in the light curve example.

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

If you want to calculate flux at specific time-frequency pairs (t_i, nu_i) instead of a grid (t_i, nu_j), you can use the alternative series interfaces:

# Define time and frequency arrays (must be the same length)
times = np.logspace(2, 8, 200)
frequencies = np.logspace(9, 17, 200)

# For time-frequency pairs (times array must be in ascending order)
results = model.flux_density(times, frequencies)

# The returned results is a FluxDict object with different flux components
print("Result attributes:", dir(results))  # Shows .fwd, .rvs, .total attributes
print("Total flux shape:", results.total.shape)  # Same shape as input arrays
print("Forward shock shape:", results.fwd.sync.shape)  # Forward shock synchrotron component

Key differences:

• flux_density_grid(): Calculates flux on a time-frequency grid (NxM output from N times and M frequencies)
• flux_density(): Calculates flux at paired time-frequency points (N output from N time-frequency pairs), requires ascending order time arrays
• flux_density_exposures(): Same as above but with exposure time averaging for realistic observational scenarios

Return value structure: All flux calculation methods return a FluxDict object with:

• .total: Combined flux from all components
• .fwd: Forward shock flux (has .sync and .ssc attributes)
• .rvs: Reverse shock flux (has .sync and .ssc attributes)

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)

Same as for light curve generation, let's set up the physical components of our afterglow model, including the environment, jet, observer, and radiation parameters:

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

medium = ISM(n_ism=1)

jet = TophatJet(theta_c=0.3, E_iso=1e52, Gamma0=100)

z = 0.1
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.1,5,10))

Get the simulation quantities (click to expand/collapse)

Now, let's get the internal simulation quantities:

# Get the simulation details over a time range
details = model.details(t_min=1e0, t_max=1e8)

# Print the available attributes
print("Simulation details attributes:", dir(details))
print("Forward shock attributes:", dir(details.fwd))

You will get a SimulationDetails object with the following structure:

Main grid coordinates:

• details.phi: 1D numpy array of azimuthal angles in radians
• details.theta: 1D numpy array of polar angles in radians
• details.t_src: 3D numpy array of source frame times on coordinate (phi_i, theta_j, t_k) grid in seconds

Forward shock details (accessed via details.fwd):

• details.fwd.t_comv: 3D numpy array of comoving times for the forward shock in seconds
• details.fwd.t_obs: 3D numpy array of observer times for the forward shock in seconds
• details.fwd.Gamma: 3D numpy array of downstream Lorentz factors for the forward shock
• details.fwd.Gamma_th: 3D numpy array of thermal Lorentz factors for the forward shock
• details.fwd.r: 3D numpy array of lab frame radii in cm
• details.fwd.B_comv: 3D numpy array of downstream comoving magnetic field strengths for the forward shock in Gauss
• details.fwd.theta: 3D numpy array of polar angles for the forward shock in radians
• details.fwd.N_p: 3D numpy array of downstream shocked proton number per solid angle for the forward shock
• details.fwd.N_e: 3D numpy array of downstream synchrotron electron number per solid angle for the forward shock
• details.fwd.gamma_a: 3D numpy array of comoving frame self-absorption Lorentz factors for the forward shock
• details.fwd.gamma_m: 3D numpy array of comoving frame injection Lorentz factors for the forward shock
• details.fwd.gamma_c: 3D numpy array of comoving frame cooling Lorentz factors for the forward shock
• details.fwd.gamma_M: 3D numpy array of comoving frame maximum Lorentz factors for the forward shock
• details.fwd.nu_a: 3D numpy array of comoving frame self-absorption frequencies for the forward shock in Hz
• details.fwd.nu_m: 3D numpy array of comoving frame injection frequencies for the forward shock in Hz
• details.fwd.nu_c: 3D numpy array of comoving frame cooling frequencies for the forward shock in Hz
• details.fwd.nu_M: 3D numpy array of comoving frame maximum frequencies for the forward shock in Hz
• details.fwd.I_nu_max: 3D numpy array of comoving frame synchrotron maximum specific intensities for the forward shock in erg/cm²/s/Hz
• details.fwd.Doppler: 3D numpy array of Doppler factors for the forward shock

Reverse shock details (accessed via details.rvs, if reverse shock is enabled):

• Similar attributes as forward shock but for the reverse shock component

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

To analyze the temporal evolution of physical parameters across different reference frames, we can visualize how key quantities evolve in the source, comoving, and observer frames. The following analysis demonstrates the comprehensive tracking of shock dynamics and microphysical parameters throughout the afterglow evolution:

Multi-parameter evolution visualization: This code creates a comprehensive multi-panel figure displaying the temporal evolution of fundamental shock parameters (Lorentz factor, magnetic field, particle numbers, radius, and peak synchrotron power) across all three reference 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)))

#plot the evolution of various parameters for phi = 0 and theta = 0 (so the first two indexes are 0)
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')

Electron energy distribution analysis: This visualization focuses specifically on the characteristic electron energies (self-absorption, injection, and cooling) across all three reference frames:

frames = ['t_src', 't_comv', 't_obs']
xlabels = [r'$t_{\rm src}$ [s]', r'$t^\prime$ [s]', r'$t_{\rm obs}$ [s]']
plt.figure(figsize= (4.2*len(frames), 3.6))

for i, frame in enumerate(frames):
    plt.subplot(1, len(frames), i + 1)
    if frame == 't_src':
        t = getattr(details, frame)
    else:
        t = getattr(details.fwd, frame)
    plt.loglog(t[0, 0, :], details.fwd.gamma_a[0, 0, :],label=r'$\gamma_a^\prime$',c='firebrick')
    plt.loglog(t[0, 0, :], details.fwd.gamma_m[0, 0, :],label=r'$\gamma_m^\prime$',c='yellowgreen')
    plt.loglog(t[0, 0, :], details.fwd.gamma_c[0, 0, :],label=r'$\gamma_c^\prime$',c='royalblue')
    plt.loglog(t[0, 0, :], details.fwd.gamma_a[0, 0, :]*details.fwd.Doppler[0,0,:]/(1+z),label=r'$\gamma_a$',ls='--',c='firebrick')
    plt.loglog(t[0, 0, :], details.fwd.gamma_m[0, 0, :]*details.fwd.Doppler[0,0,:]/(1+z),label=r'$\gamma_m$',ls='--',c='yellowgreen')
    plt.loglog(t[0, 0, :], details.fwd.gamma_c[0, 0, :]*details.fwd.Doppler[0,0,:]/(1+z),label=r'$\gamma_c$',ls='--',c='royalblue')
    plt.xlabel(xlabels[i])
    plt.ylabel(r'$\gamma_e^\prime$')
    plt.legend(ncol=2)
plt.tight_layout()
plt.savefig('electron_quantities.png', dpi=300,bbox_inches='tight')

Synchrotron frequency evolution: This analysis tracks the evolution of characteristic synchrotron frequencies:

frames = ['t_src', 't_comv', 't_obs']
xlabels = [r'$t_{\rm src}$ [s]', r'$t^\prime$ [s]', r'$t_{\rm obs}$ [s]']
plt.figure(figsize= (4.2*len(frames), 3.6))

for i, frame in enumerate(frames):
    plt.subplot(1, len(frames), i + 1)
    if frame == 't_src':
        t = getattr(details, frame)
    else:
        t = getattr(details.fwd, frame)
    plt.loglog(t[0, 0, :], details.fwd.nu_a[0, 0, :],label=r'$\nu_a^\prime$',c='firebrick')
    plt.loglog(t[0, 0, :], details.fwd.nu_m[0, 0, :],label=r'$\nu_m^\prime$',c='yellowgreen')
    plt.loglog(t[0, 0, :], details.fwd.nu_c[0, 0, :],label=r'$\nu_c^\prime$',c='royalblue')
    plt.loglog(t[0, 0, :], details.fwd.nu_a[0, 0, :]*details.fwd.Doppler[0,0,:]/(1+z),label=r'$\nu_a$',ls='--',c='firebrick')
    plt.loglog(t[0, 0, :], details.fwd.nu_m[0, 0, :]*details.fwd.Doppler[0,0,:]/(1+z),label=r'$\nu_m$',ls='--',c='yellowgreen')
    plt.loglog(t[0, 0, :], details.fwd.nu_c[0, 0, :]*details.fwd.Doppler[0,0,:]/(1+z),label=r'$\nu_c$',ls='--',c='royalblue')
    plt.xlabel(xlabels[i])
    plt.ylabel(r'$\nu$ [Hz]')
    plt.legend(ncol=2)
plt.tight_layout()
plt.savefig('photon_quantities.png', dpi=300,bbox_inches='tight')

Doppler factor spatial distribution: This polar plot visualizes the spatial distribution of the Doppler factor across the jet structure, showing how relativistic beaming varies with angular position and radial distance:

plt.figure(figsize=(6,6))
ax = plt.subplot(111, polar=True)

theta = details.fwd.theta[0,:,:]
r     = details.fwd.r[0,:,:]
D     = details.fwd.Doppler[0,:,:]

# Polar contour plot
scale = 3.0
c = ax.contourf(theta*scale, r, np.log10(D), levels=30, cmap='viridis')

ax.set_rscale('log')
true_ticks = np.linspace(0, 0.3, 6)
ax.set_xticks(true_ticks * scale)
ax.set_xticklabels([f"{t:.2f}" for t in true_ticks])
ax.set_xlim(0,0.3*scale)
ax.set_ylabel(r'$\theta$ [rad]')
ax.set_xlabel(r'$r$ [cm]')

plt.colorbar(c, ax=ax, label=r'$\log_{10} D$')

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

Equal arrival time (EAT) surface visualization: This final visualization maps the equal arrival time surfaces in polar coordinates, illustrating how light from different parts of the jet reaches the observer at the same time:

plt.figure(figsize=(6,6))
ax = plt.subplot(111, polar=True)

theta = details.fwd.theta[0,:,:]
r     = details.fwd.r[0,:,:]
t_obs = details.fwd.t_obs[0,:,:]

scale = 3.0
c = ax.contourf(theta*scale, r, np.log10(t_obs), levels=30, cmap='viridis')

ax.set_rscale('log')
true_ticks = np.linspace(0, 0.3, 6)
ax.set_xticks(true_ticks * scale)
ax.set_xticklabels([f"{t:.2f}" for t in true_ticks])
ax.set_xlim(0,0.3*scale)
ax.set_ylabel(r'$\theta$ [rad]')
ax.set_xlabel(r'$r$ [cm]')

plt.colorbar(c, ax=ax, label=r'$\log_{10} (t_{\rm obs}/s)$')

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

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. Preparing Data and 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 ObsData, Setups, Fitter, ParamDef, Scale

VegasAfterglow provides flexible options for loading observational data through the ObsData class. You can add light curves (specific flux vs. time) and spectra (specific flux vs. frequency) in multiple ways.

# Create an instance to store observational data
data = ObsData()

# Method 1: Add data directly from lists or numpy arrays

# For light curves
t_data = [1e3, 2e3, 5e3, 1e4, 2e4]  # Time in seconds
flux_data = [1e-26, 8e-27, 5e-27, 3e-27, 2e-27]  # Specific flux in erg/cm²/s/Hz
flux_err = [1e-28, 8e-28, 5e-28, 3e-28, 2e-28]  # Specific flux error in erg/cm²/s/Hz
data.add_flux_density(nu=4.84e14, t=t_data, f_nu=flux_data, err=flux_err)  # All quantities in CGS units
# You can also assign weights to each data point to account for systematic uncertainties or correlations. You don't need to worry about the weights' normalization, the code will normalize them automatically.
#data.add_flux_density(nu=4.84e14, t=t_data, f_nu=flux_data, err=flux_err, weights=np.ones_like(t_data))

# For spectra
nu_data = [...]  # Frequencies in Hz
spectrum_data = [...] # Specific flux values in erg/cm²/s/Hz
spectrum_err = [...]   # Specific flux errors in erg/cm²/s/Hz
data.add_spectrum(t=3000, nu=nu_data, f_nu=spectrum_data, err=spectrum_err, weights=np.ones_like(nu_data))  # All quantities in CGS units

# Method 2: Load from CSV files

data = ObsData()
# Define your bands and files
bands = [2.4e17, 4.84e14, 1.4e14]  # Example: X-ray, optical R-band
lc_files = ["data/ep.csv", "data/r.csv", "data/vt-r.csv"]

# Load light curves from files
for nu, fname in zip(bands, lc_files):
    df = pd.read_csv(fname)
    data.add_flux_density(nu=nu, t=df["t"], f_nu=df["Fv_obs"], err=df["Fv_err"])  # All quantities in CGS units

times = [3000] # Example: time in seconds
spec_files = ["data/ep-spec.csv"]

# Load spectra from files
for t, fname in zip(times, spec_files):
    df = pd.read_csv(fname)
    data.add_spectrum(t=t, nu=df["nu"], f_nu=df["Fv_obs"], err=df["Fv_err"])  # All quantities in CGS units

Note: The ObsData interface is designed to be flexible. You can mix and match different data sources, and add multiple light curves at different frequencies as well as multiple spectra at different times.

The Setups class defines the global properties and environment for your model. These settings remain fixed during the MCMC process. Check the documentation for all available options.

cfg = Setups()

# Source properties
cfg.lumi_dist = 3.364e28    # Luminosity distance [cm]
cfg.z = 1.58               # Redshift

# Physical model configuration
cfg.medium = "wind"        # Ambient medium: "wind", "ism", etc. (see documentation)
cfg.jet = "powerlaw"       # Jet structure: "powerlaw", "gaussian", "tophat", etc. (see documentation)

These settings affect how the model is calculated but are not varied during the MCMC process.

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

The ParamDef class is used to define the parameters for MCMC exploration. Each parameter requires a name, prior range, and sampling scale:

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
]

Scale Types:

• Scale.LOG: Sample in logarithmic space (log10) - ideal for parameters spanning multiple orders of magnitude
• Scale.LINEAR: Sample in linear space - appropriate for parameters with narrower ranges
• Scale.FIXED: Keep parameter fixed at the initial value - use for parameters you don't want to vary

Parameter Choices: The parameters you include depend on your model configuration (See documentation for all options):

• For "wind" medium: use A_star parameter
• For "ISM" medium: use n_ism parameter instead
• Different jet structures may require different parameters

Initialize the Fitter class with your data and configuration, then run the MCMC process:

# Create the fitter object
fitter = Fitter(data, cfg)

# Option 1: Nested sampling with dynesty (computes evidence, robust for multimodal posteriors)
result = fitter.fit(
    mc_params,
    resolution=(0.15, 0.5, 10),    # 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 processes
    top_k=10,                      # Number of best-fit parameters to return
    outdir="bilby_output",         # Output directory (default)
    label="afterglow_fit",         # Run label (default: "afterglow")
)

# Option 2: MCMC with emcee (faster, good for unimodal posteriors, not optimal for multimodal posteriors)
result = fitter.fit(
    mc_params,
    resolution=(0.15, 0.5, 10),    # Grid resolution (phi, theta, t)
    sampler="emcee",               # MCMC sampler
    nsteps=50000,                  # Number of steps per walker
    nburn=10000,                   # Burn-in steps to discard
    thin=1,                        # Save every nth sample
    npool=8,                       # Number of parallel processes
    top_k=10,                      # Number of best-fit parameters to return
    outdir="bilby_output",         # Output directory (default)
    label="afterglow_fit",         # Run label (default: "afterglow")
)

Important Notes:

• Parameters with Scale.LOG are sampled as log10_<name> (e.g., log10_E_iso)
• The sampler works in log10 space for LOG-scale parameters, then transforms back
• Use npool to parallelize likelihood evaluations across multiple CPU cores

The result object contains:

• samples: The posterior samples (shape: [n_samples, 1, n_params])
• labels: Parameter names (with log10_ prefix for LOG-scale params)
• latex_labels: LaTeX-formatted labels for plotting (e.g., $\log_{10}(E_{\rm iso})$)
• top_k_params: Top-k maximum likelihood parameter values
• top_k_log_probs: Log probabilities for top-k samples
• bilby_result: Full bilby Result object (for advanced diagnostics)

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

Check the top-k parameters and their uncertainties:

# Print top-k parameters (maximum likelihood)
print("Top-k parameters:")
header = f"{'Rank':>4s}  {'chi^2':>10s}  " + "  ".join(f"{name:>10s}" for name in result.labels)
print(header)
print("-" * len(header))
for i in range(result.top_k_params.shape[0]):
    chi2 = -2 * result.top_k_log_probs[i]
    vals = "  ".join(f"{val:10.4f}" for val in result.top_k_params[i])
    print(f"{i+1:4d}  {chi2:10.2f}  {vals}")

Use the best-fit parameters to generate model predictions

# Define time and frequency ranges for predictions
t_out = np.logspace(2, 9, 150)

nu_out = np.logspace(16,20,150)

best_params = result.top_k_params[0]

# Generate model light curves at the specified bands using the best-fit parameters
lc = fitter.flux_density_grid(best_params, t_out, band)

# Generate model spectra at the specified times using the best-fit parameters
spec = fitter.flux_density_grid(best_params, times, nu_out)

Now you can plot the best-fit model:

# Function to plot model light curves along with observed data
def draw_bestfit(t,lc_fit, nu, spec_fit):
    fig =plt.figure(figsize=(4.5, 7.5))

    ax1 = fig.add_subplot(211)
    ax2 = fig.add_subplot(212)

    shift = [1,1,200]
    colors = ['blue', 'orange', 'green']
    for i, file, sft, c in zip(range(len(lc_files)), lc_files, shift, colors ):
        df = pd.read_csv(file)
        ax1.errorbar(df["t"], df["Fv_obs"]*sft, df["Fv_err"]*sft, fmt='o',markersize=4,label=file, color=c,markeredgecolor='k', markeredgewidth=.4)
        ax1.plot(t, np.array(lc_fit[i,:])*sft, color=c,lw=1)

    ax1.set_xscale('log')
    ax1.set_yscale('log')
    ax1.set_xlabel('t [s]')
    ax1.set_ylabel(r'$F_\nu$ [erg/cm$^2$/s/Hz]')
    ax1.legend()

    for i, file, sft, c in zip(range(len(spec_files)), spec_files, shift, colors ):
        df = pd.read_csv(file)
        ax2.errorbar(df["nu"], df["Fv_obs"]*sft, df["Fv_err"]*sft, fmt='o',markersize=4,label=file, color=c,markeredgecolor='k', markeredgewidth=.4)
        ax2.plot(nu, np.array(spec_fit[:,i])*sft, color=c,lw=1)

    ax2.set_xscale('log')
    ax2.set_yscale('log')
    ax2.set_xlabel(r'$\nu$ [Hz]')
    ax2.set_ylabel(r'$F_\nu$ [erg/cm$^2$/s/Hz]')
    ax2.legend()
    plt.tight_layout()

draw_bestfit(t_out, lc, nu_out, spec)

Corner plots are essential for visualizing parameter correlations and posterior distributions:

def plot_corner(flat_chain, labels, filename="corner_plot.png"):
    fig = corner.corner(
        flat_chain,
        labels=labels,
        quantiles=[0.16, 0.5, 0.84],  # For median and ±1σ
        show_titles=True,
        title_kwargs={"fontsize": 14},
        label_kwargs={"fontsize": 14},
        truths=np.median(flat_chain, axis=0),  # Show median values
        truth_color='red',
        bins=30,
        smooth=1,
        fill_contours=True,
        levels=[0.16, 0.5, 0.68],  # 1σ and 2σ contours
        color='k'
    )
    fig.savefig(filename, dpi=300, bbox_inches='tight')

# Create the corner plot
flat_chain = result.samples.reshape(-1, result.samples.shape[-1])
plot_corner(flat_chain, result.latex_labels)

Parameter Fitting with Redback

For parameter estimation and Bayesian inference, VegasAfterglow models are integrated into redback, which provides a unified interface for fitting all transient types.

Loading Data & Fitting: (click to expand/collapse)

import redback

# Multiple ways to load data:
GRB = 'GRB070809'

# Method 1: From Swift BAT+XRT (recommended for Swift GRBs)
redback.get_data.get_bat_xrt_afterglow_data_from_swift(grb=GRB, data_mode="flux")
afterglow = redback.afterglow.SGRB.from_swift_grb(name=GRB, data_mode='flux')

# Method 2: From open access catalogs
afterglow = redback.afterglow.Afterglow.from_open_access_catalogue(
    GRB, data_mode='flux_density'
)

# Method 3: From your own data files
import pandas as pd
data = pd.read_csv('my_grb_data.csv')
afterglow = redback.transient.Afterglow(
    name=GRB,
    data_mode='flux_density',
    time=data['time'].values,
    flux_density=data['flux'].values,
    flux_density_err=data['flux_err'].values,
    frequency=data['frequency'].values
)

# Fit with VegasAfterglow tophat model
result = redback.fit_model(
    transient=afterglow,
    model='vegas_tophat',  # VegasAfterglow tophat jet model
    sampler='dynesty',
    nlive=1000
)

result.plot_corner()
result.plot_lightcurve()
result.plot_residuals()

Available VegasAfterglow models in redback:

• vegas_tophat - Tophat jet + ISM
• vegas_gaussian - Gaussian jet + ISM
• vegas_powerlaw - Power-law jet + ISM
• vegas_two_component - Two-component jet + ISM
• vegas_tophat_wind - Tophat jet + Wind medium
• ... and more

Why use the Redback interface?

• Data Management: Seamlessly load data from Swift, Fermi, BATSE, or custom files. Supports flux, flux_density, magnitude, and luminosity.
• Advanced Statistics: Access multiple samplers (dynesty, emcee, ultranest), parallel sampling, and Bayesian model comparison.
• Visualization: Generate publication-ready corner plots and multi-band light curves automatically.

For complete documentation on the API, visit the redback documentation.

---

Validation & Testing

VegasAfterglow includes a comprehensive validation framework to ensure numerical accuracy and physical correctness. The validation suite consists of two main components:

Benchmark Tests (click to expand/collapse)

Benchmark tests measure computational performance and verify numerical convergence across resolution parameters:

• Performance Timing: Measures execution time for various jet/medium/radiation configurations
• Resolution Convergence: Tests convergence in phi (azimuthal), theta (polar), and time dimensions
• Pass Criteria: Mean error < 5% and max error < 10% at fiducial resolution

The default fiducial resolution (0.15, 0.5, 10) has been tested to converge for most configurations. For details, see the convergence report PDF generated by the validation suite (validation/output/validation_report.pdf).

Regression Tests (click to expand/collapse)

Regression tests verify that simulation outputs match theoretical predictions from GRB afterglow theory:

• Shock Dynamics: Validates power-law scaling of Lorentz factor, radius, magnetic field, and particle number
• Characteristic Frequencies: Verifies evolution of injection (nu_m) and cooling (nu_c) frequencies
• Spectral Shapes: Checks power-law indices across different frequency regimes (I-V)
• Evolutionary Phases: Tests coasting, Blandford-McKee, and Sedov-Taylor phases for ISM and wind media

Running Validation (click to expand/collapse)

Prerequisite: Install with validation support, which includes per-stage CPU profiling and PDF report dependencies:

pip install -e ".[test]" --config-settings=cmake.define.AFTERGLOW_PROFILE=ON

# Run full validation suite (benchmark + regression + PDF report)
python validation/run_validation.py --all

# Run only benchmark tests
python validation/run_validation.py --benchmark

# Run only regression tests
python validation/run_validation.py --regression

# Check existing results without re-running tests
python validation/run_validation.py --check-only

The validation runner generates a comprehensive PDF report at validation/comprehensive_report.pdf with convergence plots, summary grids, and detailed diagnostics. The overview page includes a stacked bar chart breaking down CPU time by internal C++ computation stage.

To return to the normal (zero-overhead) build:

pip install -e .

Powered by Claude Code

---

Documentation

Comprehensive documentation is available at Documentation including:

• Installation Guide: Detailed instructions for setting up VegasAfterglow
• Quick Start: Get up and running with basic examples
• Examples: Practical examples showing common use cases
• MCMC Fitting: Complete guide to Bayesian parameter estimation with bilby integration
• Parameter Reference: Comprehensive reference for all physical and numerical parameters
• Validation & Testing: Benchmark and regression test framework for verifying numerical accuracy
• Python API Reference: Complete documentation of the Python interface
• C++ API Reference: Detailed documentation of C++ classes and functions
• Troubleshooting: Common issues and solutions
• Contributing Guide: Information for developers who wish to contribute

For a complete history of changes and new features, see our 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.

The BSD 3-Clause License is a permissive open source license that allows you to:

• Freely use, modify, and distribute the software in source and binary forms
• Use the software for commercial purposes
• Integrate the software into proprietary applications

Requirements:

• You must include the original copyright notice and the license text
• You cannot use the names of the authors or contributors to endorse derived products
• The license provides no warranty or liability protection

For the full license text, see the LICENSE file in the repository.

---

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 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}
}

---

Parts of the documentation and code comments were generated with the assistance of Claude Code.