Examples

Basic Usage

Setting up a simple afterglow model

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

# Define the circumburst environment (constant density ISM)
medium = ISM(n_ism=1)

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

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

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

# Combine all components into a complete afterglow model
model = Model(jet=jet, medium=medium, observer=obs, fwd_rad=rad, resolutions=(0.3,1,10))

# Define time range for light curve calculation
times = np.logspace(2, 8, 200)

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

# Calculate the afterglow emission at each time and frequency
# NOTE that the times array needs to be in ascending order
results = model.flux_density_grid(times, bands)

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

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

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

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

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

# Calculate spectra at each epoch
results = model.flux_density_grid(epochs, frequencies)

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

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

Calculate flux on time-frequency pairs

Suppose you want to calculate the flux at specific time-frequency pairs (t_i, nu_i) instead of a grid (t_i, nu_j), you can use the following method:

# Define time range for light curve calculation
times = np.logspace(2, 8, 200)

# Define observing frequencies (must be the same length as times)
bands = np.logspace(9,17, 200)

results = model.flux_density(times, bands) #times array could be random order

# the returned results is a FluxDict object with arrays of the same shape as the input times and bands.

Calculate flux with exposure time averaging

For observations with finite exposure times, you can calculate time-averaged flux by sampling multiple points within each exposure:

# Define observation times (start of exposure)
times = np.logspace(2, 8, 50)

# Define observing frequencies (must be the same length as times)
bands = np.logspace(9, 17, 50)

# Define exposure times for each observation (in seconds)
expo_time = np.ones_like(times) * 100  # 100-second exposures

# Calculate time-averaged flux with 20 sample points per exposure
results = model.flux_density_exposures(times, bands, expo_time, num_points=20)

# The returned results is a FluxDict object with arrays of the same shape as input times and bands
# Each flux value represents the average over the corresponding exposure time

Note

The function samples num_points evenly spaced within each exposure time and averages the computed flux. Higher num_points gives more accurate time averaging but increases computation time. The minimum value is 2.

Ambient Media Models

Wind Medium

from VegasAfterglow import Wind

# Create a stellar wind medium
wind = Wind(A_star=0.1)  # A* parameter

#..other settings
model = Model(medium=wind, ...)

Stratified Medium

from VegasAfterglow import Wind

# Create a stratified stellar wind medium;
# smooth transited stratified medium. Inner region, n(r) = n0, middle region n(r) \propto 1/r^2, outer region n(r)=n_ism
# A = 0 (default): fallback to n = n_ism
# n0 = inf (default): wind bubble, from wind profile to ism profile
# A = 0 & n0 = inf: pure wind;
wind = Wind(A_star=0.1, n_ism = 1, n0 = 1e-3)

#..other settings
model = Model(medium=wind, ...)

User-Defined Medium

from VegasAfterglow import Medium

mp = 1.67e-24 # proton mass in gram

# Define a custom density profile function
def density(phi, theta, r):# r in cm, phi and theta in radians
    return mp # n_ism =  1 cm^-3
    #return whatever density profile (cm^-3) you want as a function of phi, theta, and r

# Create a user-defined medium
medium = Medium(rho=density)

#..other settings
model = Model(medium=medium, ...)

Jet Models

Gaussian Jet

from VegasAfterglow import GaussianJet

# Create a structured jet with Gaussian energy profile
jet = GaussianJet(
    theta_c=0.05,         # Core angular size (radians)
    E_iso=1e53,           # Isotropic-equivalent energy (ergs)
    Gamma0=300            # Initial Lorentz factor
)

#..other settings
model = Model(jet=jet, ...)

Power-Law Jet

from VegasAfterglow import PowerLawJet

# Create a power-law structured jet
jet = PowerLawJet(
    theta_c=0.05,         # Core angular size (radians)
    E_iso=1e53,           # Isotropic-equivalent energy (ergs)
    Gamma0=300,           # Initial Lorentz factor
    k_e=2.0,              # Power-law index for energy angular dependence
    k_g=2.0               # Power-law index for Lorentz factor angular dependence
)

#..other settings
model = Model(jet=jet, ...)

Two-Component Jet

from VegasAfterglow import TwoComponentJet

# Create a two-component jet
jet = TwoComponentJet(
    theta_c=0.05,        # Narrow component angular size (radians)
    E_iso=1e53,          # Isotropic-equivalent energy of the narrow component (ergs)
    Gamma0=300,          # Initial Lorentz factor of the narrow component
    theta_w=0.1,         # Wide component angular size (radians)
    E_iso_w=1e52,        # Isotropic-equivalent energy of the wide component (ergs)
    Gamma0_w=100         # Initial Lorentz factor of the wide component
)

#..other settings
model = Model(jet=jet, ...)

Step Power-Law Jet

from VegasAfterglow import StepPowerLawJet

# Create a step power-law structured jet (uniform core with sharp transition)
jet = StepPowerLawJet(
    theta_c=0.05,        # Core angular size (radians)
    E_iso=1e53,          # Isotropic-equivalent energy of the core component (ergs)
    Gamma0=300,          # Initial Lorentz factor of the core component
    E_iso_w=1e52,        # Isotropic-equivalent energy of the wide component (ergs)
    Gamma0_w=100,        # Initial Lorentz factor of the wide component
    k_e=2.0,             # Power-law index for energy angular dependence
    k_g=2.0              # Power-law index for Lorentz factor angular dependence
)

#..other settings
model = Model(jet=jet, ...)

Jet with Spreading

from VegasAfterglow import TophatJet

jet = TophatJet(
    theta_c=0.05,
    E_iso=1e53,
    Gamma0=300,
    spreading=True       # Enable spreading
)

#..other settings
model = Model(jet=jet, ...)

Note

The jet spreading (Lateral Expansion) is experimental and only works for the top-hat jet, Gaussian jet, and power-law jet with a jet core. The spreading prescription may not work for arbitrary user-defined jet structures.

Magnetar Spin-down

from VegasAfterglow import Magnetar

# Create a tophat jet with magnetar spin-down energy injection; Luminosity 1e46 erg/s, t_0 = 100 seconds, and q = 2
jet = TophatJet(theta_c=0.05, E_iso=1e53, Gamma0=300, magnetar=Magnetar(L0=1e46, t0=100, q=2))

Note

The magnetar spin-down injection is implemented in the default form L0*(1+t/t0)^(-q) for theta < theta_c. You can pass the magnetar argument to the power-law and Gaussian jet as well.

User-Defined Jet

You may also define your own jet structure by providing the energy and lorentz factor profile. Those two profiles are required to complete a jet structure. You may also provide the magnetization profile, enregy injection profile, and mass injection profile. Those profiles are optional and will be set to zero function if not provided.

from VegasAfterglow import Ejecta

# Define a custom energy profile function, required to complete the jet structure
def E_iso_profile(phi, theta):
    return 1e53  # E_iso = 1e53 erg isotropic fireball
    #return whatever energy profile you want as a function of phi and theta in unit of erg [not erg per solid angle]

# Define a custom lorentz factor profile function, required to complete the jet structure
def Gamma0_profile(phi, theta):
    return 300 # Gamma0 = 300
    #return whatever lorentz factor profile you want as a function of phi and theta

# Define a custom magnetization profile function, optional
def sigma0_profile(phi, theta):
    return 0.1 # sigma = 0.1
    #return whatever magnetization profile you want as a function of phi and theta

# Define a custom energy injection profile function, optional
def E_dot_profile(phi, theta, t):
    return 1e46 * (1 + t / 100)**(-2) # L = 1e46 erg/s, t0 = 100 seconds
    #return whatever energy injection  profile you want as a function of phi, theta, and time in unit of erg/s [not erg/s per solid angle]

# Define a custom mass injection profile function, optional
def M_dot_profile(phi, theta, t):
    #return whatever mass injection profile you want as a function of phi, theta, and time in unit of g/s [not g/s per solid angle]

# Create a user-defined jet
jet = Ejecta(E_iso=E_iso_profile, Gamma0=Gamma0_profile, sigma0=sigma0_profile, E_dot=E_dot_profile, M_dot=M_dot_profile)

#..other settings

#if your jet is not axisymmetric, set axisymmetric to False
model = Model(jet=jet, ..., axisymmetric=False, resolutions=(0.3, 1, 10))

# the user-defined jet structure could be spiky, the default resolution may not resolve the jet structure. if that is the case, you can try a finer resolution (phi_ppd, theta_ppd, t_ppd)
# where phi_ppd is the number of points per degree in the phi direction, theta_ppd is the number of points per degree in the theta direction, and t_ppd is the number of points per decade in the time direction    .

Note

Setting user-defined structured jet in the Python level is OK for light curve and spectrum calculation. However, it is not recommended for MCMC parameter fitting if you do care about the performance. The reason is that setting user-defined profiles in the Python level leads to a large overhead due to the Python-C++ inter-process communication. Users are recommended to set up the user-defined jet structure in the C++ level for MCMC parameter fitting for better performance, if you want the best performance.

Radiation Processes

Reverse Shock Emission

from VegasAfterglow import Radiation

#set the jet duration to be 100 seconds, the default is 1 second. The jet duration affects the reverse shock thickness (thin shell or thick shell).
jet = TophatJet(theta_c=0.1, E_iso=1e52, Gamma0=300, duration = 100)

# Create a radiation model with both forward and reverse shock synchrotron radiation
fwd_rad = Radiation(eps_e=1e-1, eps_B=1e-3, p=2.3)
rvs_rad = Radiation(eps_e=1e-2, eps_B=1e-4, p=2.4)

#..other settings
model = Model(fwd_rad=fwd_rad, rvs_rad=rvs_rad, resolutions=(0.5, 1, 10),...)

times = np.logspace(2, 8, 200)

bands = np.array([1e9, 1e14, 1e17])

results = model.flux_density_grid(times, bands)

plt.figure(figsize=(4.8, 3.6),dpi=200)

# 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.fwd.sync[i,:], label=fr'${base:.1f} \times 10^{{{exp}}}$ Hz (fwd)')
    plt.loglog(times, results.rvs.sync[i,:], label=fr'${base:.1f} \times 10^{{{exp}}}$ Hz (rvs)')#reverse shock synchrotron

Note

You may increase the resolution of the grid to improve the accuracy of the reverse shock synchrotron radiation if you see spiky features.

Inverse Compton Cooling

from VegasAfterglow import Radiation

# Create a radiation model with inverse Compton cooling (without Klein-Nishina correction) on synchrotron radiation
rad = Radiation(eps_e=1e-1, eps_B=1e-3, p=2.3, ssc_cooling=True, kn=False)

#..other settings
model = Model(fwd_rad=rad, ...)

Self-Synchrotron Compton Radiation

from VegasAfterglow import Radiation

# Create a radiation model with self-Compton radiation
rad = Radiation(eps_e=1e-1, eps_B=1e-3, p=2.3, ssc=True, kn=True, ssc_cooling=True)

#..other settings
model = Model(fwd_rad=rad, ...)

times = np.logspace(2, 8, 200)

bands = np.array([1e9, 1e14, 1e17])

results = model.flux_density_grid(times, bands)

plt.figure(figsize=(4.8, 3.6),dpi=200)

# 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.fwd.sync[i,:], label=fr'${base:.1f} \times 10^{{{exp}}}$ Hz (sync)')#synchrotron
    plt.loglog(times, results.fwd.ssc[i,:], label=fr'${base:.1f} \times 10^{{{exp}}}$ Hz (SSC)')#SSC

Note

(ssc_cooling = False, kn = False, ssc = True): The IC radiation is calculated based on synchrotron spectrum without IC cooling.

(ssc_cooling = True, kn = False, ssc = True): The IC radiation is calculated based on synchrotron spectrum with IC cooling, but without Klein-Nishina correction.

(ssc_cooling = True, kn = True, ssc = True): The IC radiation is calculated based on synchrotron spectrum with both IC cooling and Klein-Nishina correction.