VegasAfterglow 2.0.0rc2
High-performance C++ framework with Python interface for modeling Gamma-Ray Burst (GRB) afterglows.
License: BSD-3-Clause
Requires Python: >=3.8

# VegasAfterglow

<img align="left" src="https://github.com/YihanWangAstro/VegasAfterglow/raw/main/assets/logo.svg" alt="VegasAfterglow Logo" width="350"/>

[![C++ Version](https://img.shields.io/badge/C%2B%2B-20-blue.svg)](https://en.cppreference.com/w/cpp/20)
[![PyPI version](https://img.shields.io/pypi/v/VegasAfterglow.svg)](https://pypi.org/project/VegasAfterglow/)
[![Build Status](https://github.com/YihanWangAstro/VegasAfterglow/actions/workflows/PyPI-build.yml/badge.svg)](https://github.com/YihanWangAstro/VegasAfterglow/actions/workflows/PyPI-build.yml)
[![License](https://img.shields.io/badge/License-BSD--3--Clause-green.svg)](LICENSE)
[![Platform](https://img.shields.io/badge/Platform-Linux%20|%20macOS%20|%20Windows-lightgrey.svg)]()
[![Python Version](https://img.shields.io/badge/Python-3.8%2B-blue.svg)](https://www.python.org/)
[![Documentation](https://img.shields.io/badge/Documentation-Online-brightgreen.svg)](https://yihanwangastro.github.io/VegasAfterglow/docs/index.html)
[![Validation Report](https://img.shields.io/badge/Validation-Report-blue.svg)](https://yihanwangastro.github.io/VegasAfterglow/reports/latest/comprehensive_report.pdf)

<div align="left">

**[Latest Release Notes](CHANGELOG.md#v200-beta---2026-02-02)** | **[Validation Report (PDF)](https://yihanwangastro.github.io/VegasAfterglow/reports/latest/comprehensive_report.pdf)** | **[Install Now](#installation)**
</div>

<p align="justify">
<b>VegasAfterglow</b> 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.
</p>

<h3 align="center">Ecosystem & Integration</h3>

* **Advanced Inference:** While VegasAfterglow includes a standalone fitter, its models are also integrated into [**Redback**](https://github.com/nikhil-sarin/redback) for more advanced Bayesian inference, model comparison, and multi-messenger analysis.
* **Community Tools:** We also recommend checking out [**PyFRS**](https://github.com/leiwh/PyFRS), [**Jetsimpy**](https://github.com/haowang-astro/jetsimpy), [**ASGARD**](https://github.com/mikuru1096/ASGARD_GRBAfterglow) as tool for afterglow modeling.
<br clear="left"/>

---

## Table of Contents

- [VegasAfterglow](#vegasafterglow)
  - [Table of Contents](#table-of-contents)
  - [Features](#features)
  - [Performance Highlights](#performance-highlights)
  - [Installation](#installation)
    - [Python Installation](#python-installation)
    - [C++ Installation](#c-installation)
  - [Usage](#usage)
    - [Quick Start](#quick-start)
    - [Light Curve \& Spectrum Calculation](#light-curve--spectrum-calculation)
    - [Internal Quantities Evolution](#internal-quantities-evolution)
    - [MCMC Parameter Fitting](#mcmc-parameter-fitting)
    - [Parameter Fitting with **Redback**](#parameter-fitting-with-redback)
  - [Validation \& Testing](#validation--testing)
  - [Documentation](#documentation)
  - [Contributing](#contributing)
  - [License](#license)
  - [Acknowledgments \& Citation](#acknowledgments--citation)

---

## Features

<h3 align="center">Shock Dynamics</h3>

<img align="right" src="https://github.com/YihanWangAstro/VegasAfterglow/raw/main/assets/shock_dynamics.svg" width="450"/>

- **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](https://yihanwangastro.github.io/VegasAfterglow/docs/examples.html#user-defined-medium)**.
- **Energy and Mass Injection:** Supports user-defined profiles for continuous energy and/or mass injection into the blast wave.

<br clear="right"/>

<h3 align="center">Jet Structure & Geometry</h3>

<img align="right" src="https://github.com/YihanWangAstro/VegasAfterglow/raw/main/assets/jet_geometry.svg" width="450"/>

- **Structured Jet Profiles:** Allows **[user-defined angular profiles](https://yihanwangastro.github.io/VegasAfterglow/docs/examples.html#user-defined-jet)** 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.

<br clear="right"/>

<h3 align="center">Radiation Mechanisms</h3>

<img align="right" src="https://github.com/YihanWangAstro/VegasAfterglow/raw/main/assets/radiation_mechanisms.svg" width="450"/>

- **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.

<br clear="right"/>

---

## Performance Highlights

<img align="left" src="https://github.com/YihanWangAstro/VegasAfterglow/raw/main/assets/convergence_plot.png" width="450"/>

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.

<br clear="left"/>

---

## Installation

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

### Python Installation

To install VegasAfterglow using pip:

```bash
pip install VegasAfterglow
```

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

```bash
pip install VegasAfterglow[mcmc]
```

VegasAfterglow requires Python 3.8 or higher.

<details>
<summary><b>Alternative: Install from Source</b> <i>(click to expand/collapse)</i></summary>
<br>

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

1. Clone this repository:

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

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

```bash
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).
</details>

### C++ Installation

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

<details>
<summary><b>Instructions for C++ Installation</b> <i>(click to expand/collapse)</i></summary>
<br>

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

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

2. Compile and run tests:

```bash
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++](#creating-custom-problem-generators-with-c) section or examine the example problem generators in `tests/demo/`.

<details>
<summary><b>Build Prerequisites</b> <i>(click to expand for dependency information)</i></summary>
<br>

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)

</details>
</details>

---

## Usage

### Quick Start

**Get your first light curve in under 10 lines:**

```python
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`.

<details>
<summary><b>Model Setup</b> <i>(click to expand/collapse)</i></summary>
<br>

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

```python
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)
```

</details>

<details>
<summary><b>Light Curve Calculation</b> <i>(click to expand/collapse)</i></summary>
<br>

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

```python
# 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')
```

<div align="center">
<img src="assets/quick-lc.png" alt="Afterglow Light Curves" width="600"/>

Running the light curve script will produce this figure showing the afterglow evolution across different frequencies.
</div>
</details>

<details>
<summary><b>Spectrum Analysis</b> <i>(click to expand/collapse)</i></summary>
<br>

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

```python
# 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')
```

<div align="center">
<img src="assets/quick-spec.png" alt="Broadband Spectra" width="600"/>

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.
</div>
</details>

<details>
<summary><b>Time-Frequency Pairs Calculation</b> <i>(click to expand/collapse)</i></summary>
<br>

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:

```python
# 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)

</details>



### 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`.

<details>
<summary><b>Model Setup</b> <i>(click to expand/collapse)</i></summary>
<br>

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:

```python
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.15,0.5,10))
```

</details>

<details>
<summary><b>Get the simulation quantities</b> <i>(click to expand/collapse)</i></summary>
<br>

Now, let's get the internal simulation quantities:

```python

# 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

</details>

<details>
<summary><b>Checking the evolution of various parameters</b> <i>(click to expand/collapse)</i></summary>
<br>

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:

```python
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')
```

<div align="center">
<img src="assets/shock_quantities.png" alt="Shock evolution" width="1000"/>
</div>

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

```python
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')
```

<div align="center">
<img src="assets/electron_quantities.png" alt="Shock evolution" width="1000"/>
</div>

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

```python
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')
```
<div align="center">
<img src="assets/photon_quantities.png" alt="Shock evolution" width="1000"/>
</div>

**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:

```python
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')
```
<div align="center">
<img src="assets/doppler.png" alt="Shock evolution" width="600"/>
</div>

**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:

```python
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')
```
<div align="center">
<img src="assets/EAT.png" alt="Shock evolution" width="600"/>

</div>
</details>


### 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.

<details>
<summary><b>1. Configuring the Model</b> <i>(click to expand/collapse)</i></summary>
<br>

```python
# 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](https://yihanwangastro.github.io/VegasAfterglow/docs/index.html) for all available options.
</details>

<details>
<summary><b>2. Defining Parameters and Running MCMC</b> <i>(click to expand/collapse)</i></summary>
<br>

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

```python
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 configuration and add data directly:

```python
# 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.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 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.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 threads
    top_k=10,                      # Number of best-fit parameters to return
)
```

**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)

</details>

<details>
<summary><b>3. Analyzing Results and Generating Predictions</b> <i>(click to expand/collapse)</i></summary>
<br>

Check the top-k parameters and their uncertainties:

```python
# 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

```python
# 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, bands)

# 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:

```python
# 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:

```python
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)
```

</details>

### Parameter Fitting with [**Redback**](https://github.com/nikhil-sarin/redback)

For parameter estimation and Bayesian inference, VegasAfterglow models are integrated into [**redback**](https://github.com/nikhil-sarin/redback), which provides a unified interface for fitting all transient types.

<details>
<summary><b> Loading Data & Fitting:</b> <i>(click to expand/collapse)</i></summary>

```python
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


</details>

**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**](https://redback.readthedocs.io/en/latest/?badge=latest).
<br>


---

## Validation & Testing

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

<details>
<summary><b>Benchmark Tests</b> <i>(click to expand/collapse)</i></summary>
<br>

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 [Validation Report (PDF)](https://yihanwangastro.github.io/VegasAfterglow/reports/latest/comprehensive_report.pdf).

</details>

<details>
<summary><b>Regression Tests</b> <i>(click to expand/collapse)</i></summary>
<br>

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

</details>

<details>
<summary><b>Running Validation</b> <i>(click to expand/collapse)</i></summary>
<br>

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

```bash
pip install -e ".[test]" --config-settings=cmake.define.AFTERGLOW_PROFILE=O