1. Intro to Waveforms

This tutorial demonstrates how to use PyART’s waveform class

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
import seaborn as sns

# logging configuration, default INFO level
from PyART.logging_config import setup_logging
setup_logging()

# Set up Seaborn aesthetics
sns.set_context('talk') 
sns.set_theme(font_scale=1.2)
sns.set_style('ticks')

# Update matplotlib rcParams
rcParams.update(
    {
        'text.usetex': False,
        'font.family': 'stixgeneral',
        'mathtext.fontset': 'stix',
        'axes.grid': True,
        'grid.linestyle': ':',
        'grid.color': '#bbbbbb',
        'axes.linewidth': 1,
    }
)

We now create an empty Waveform object and inspect its attributes. This class is inherited by all waveform catalogs and models, meaning that they will all have these attributes and methods.
For more information, see the Waveform class documentation.

from PyART.waveform import Waveform

# Create an empty Waveform object
waveform = Waveform()

attributes = [attr for attr in dir(waveform) if not attr.startswith('_') and not callable(getattr(waveform, attr))]
methods = [method for method in dir(waveform) if not method.startswith('_') and callable(getattr(waveform, method))]

print("Attributes of the Waveform object:")
for attr in attributes:
    print(f"- {attr}")

print("\nMethods of the Waveform object:")
for method in methods:
    print(f"- {method}")

Attributes of the Waveform object:
- dothlm
- dyn
- f
- hc
- hlm
- hp
- kind
- psi4lm
- t
- t_psi4
- u

Methods of the Waveform object:
- compute_dothlm
- compute_hphc
- compute_psi4lm
- cut
- dynamics_from_hlm
- ej_from_hlm
- find_max
- integrate_data
- interpolate_hlm
- multiply_by
- to_frequency

We now load a waveform from the SXS catalog and plot its modes \(h_{\ell m}\).
They are stored as a dictionary of the form: hlm[(l,m)]. Each dictionary entry has the amplitude: A and phase p of the mode, as well as the real and imag part of \(h_{\ell,m}\).

Note that, by default, the modes are not interpolated to a uniform time grid. Therefore, before plotting, we:

  • Interpolate to a uniform time grid with waveform.interpolate_hlm

  • Shift the time array so that the (2,2) mode peaks at \(t=0M\), using the find_max method.

from PyART.catalogs import sxs

# Utility function to check if an array is uniformly spaced
def is_uniform(arr, rtol=1e-10, atol=1e-12):
    delta = np.diff(arr)
    return np.allclose(delta, delta[0], rtol=rtol, atol=atol)

sxs_waveform = sxs.Waveform_SXS(ID='0180',
                                download=True,
                                ignore_deprecation=True,
                                downloads=["hlm", "metadata", "horizons"],
                                load=["hlm", "metadata", "horizons"]
                                )


# verify that the time array is uniform
print("Is the time array uniform?", is_uniform(sxs_waveform.u))

# interpolate to uniform time gri, dT = 1M
t_interp, hlm_interp = sxs_waveform.interpolate_hlm(1)

# overwrite time array and modes in the waveform object
sxs_waveform._t = t_interp
sxs_waveform._u = t_interp
sxs_waveform._hlm = hlm_interp

# verify that the time array is uniform
print("Is the time array uniform?", is_uniform(sxs_waveform.u))

t_max, _, _, _, mrg_idx = sxs_waveform.find_max(return_idx=True, mode=(2,2))

fig, ax = plt.subplots()
modes = [(2, 2), (2, 1), (3, 3), (4, 4), (2, -2), (2, -1), (3, -3), (4, -4)]
for mode in modes:
    ax.semilogy(sxs_waveform.u-t_max, sxs_waveform.hlm[mode]['A'], label=f'{mode}')

# legend, out of box
ax.legend(ncol=1, bbox_to_anchor=(1.05, 1), loc='upper left')
ax.axvline(0, color='k', linestyle='--', alpha=0.5)
ax.set_xlabel('t/M')
ax.set_ylabel(r'$h_{\ell m}/\nu$')
plt.show()

2026-03-17 19:34:27 The path ../dat/SXS/SXS_BBH_0180 does not exist or contains no 'Lev*' directory.

2026-03-17 19:34:27 Downloading the simulation from the SXS catalog.

2026-03-17 19:34:28 Setting the download (cache) directory to ../dat/SXS/

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[3], line 8
      5     delta = np.diff(arr)
      6     return np.allclose(delta, delta[0], rtol=rtol, atol=atol)
----> 8 sxs_waveform = sxs.Waveform_SXS(ID='0180',
      9                                 download=True,
     10                                 ignore_deprecation=True,
     11                                 downloads=["hlm", "metadata", "horizons"],
     12                                 load=["hlm", "metadata", "horizons"]
     13                                 )
     16 # verify that the time array is uniform
     17 print("Is the time array uniform?", is_uniform(sxs_waveform.u))

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/catalogs/sxs.py:135, in Waveform_SXS.__init__(self, path, ID, order, level, cut_N, cut_U, ellmax, load, download, downloads, load_m0, nu_rescale, src, ignore_deprecation, basename)
    131     logging.info(
    132         f"The path {self.sxs_data_path} does not exist or contains no 'Lev*' directory."
    133     )
    134     logging.info("Downloading the simulation from the SXS catalog.")
--> 135     self.download_simulation(
    136         ID=self.ID,
    137         path=path,
    138         downloads=downloads,
    139         level=self.level,
    140         ignore_deprecation=ignore_deprecation,
    141         extrapolation_order=order,
    142     )
    143 else:
    144     logging.warning(
    145         "Use download=True to download the simulation from the SXS catalog."
    146     )

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/PyART/catalogs/sxs.py:290, in Waveform_SXS.download_simulation(self, ID, path, downloads, level, ignore_deprecation, extrapolation_order)
    288 sys.stdout = open(os.devnull, "w")
    289 sys.stderr = open(os.devnull, "w")
--> 290 sxs_sim = sxsmod.load(
    291     name_level,
    292     extrapolation_order=extrapolation_order,
    293     ignore_deprecation=ignore_deprecation,
    294     progress=True,
    295 )
    296 logging.info(f"Loaded SXS simulation {name_level}.")
    298 # Set Level if not already set

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/sxs/handlers.py:288, in load(location, download, cache, progress, truepath, **kwargs)
    285     return sxscatalog.load(location, download=download, **kwargs)
    287 elif sxs_id_version_lev_exact_re.match(location):
--> 288     return Simulation(location, download=download, cache=cache, progress=progress, **kwargs)
    290 else:
    291     # Try to find an appropriate SXS file in the simulations
    292     simulations = Simulations.load(
    293         download=download,
    294         local=kwargs.get("local", False),
    295         annex_dir=kwargs.get("annex_dir", None)
    296     )

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/sxs/simulations/simulation.py:135, in Simulation(location, *args, **kwargs)
    132 from ..metadata.metric import MetadataMetric
    134 # Load the simulation catalog
--> 135 simulations = load("simulations")
    136 v = Version(simulations.tag)
    137 latest_version = f"v{v.major}.{v.minor}"

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/sxs/handlers.py:285, in load(location, download, cache, progress, truepath, **kwargs)
    282     return Catalog.load(download=download)
    284 elif location in ["simulations", "dataframe"]:
--> 285     return sxscatalog.load(location, download=download, **kwargs)
    287 elif sxs_id_version_lev_exact_re.match(location):
    288     return Simulation(location, download=download, cache=cache, progress=progress, **kwargs)

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/sxscatalog/__init__.py:16, in load(location, download, **kwargs)
     13 from .simulations import Simulations
     15 if location == "simulations":
---> 16     return Simulations.load(
     17         download=download,
     18         tag=kwargs.get("tag", ""),
     19         local=kwargs.get("local", False),
     20         annex_dir=kwargs.get("annex_dir", None),
     21         output_file=kwargs.get("output_file", None),
     22         compute_md5=kwargs.get("compute_md5", False),
     23         show_progress=kwargs.get("show_progress", False),
     24         ignore_cached=kwargs.get("ignore_cached", False),
     25     )
     27 elif location == "dataframe":
     28     return load("simulations", download=download, **kwargs).dataframe

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/sxscatalog/simulations/simulations.py:415, in Simulations.load(cls, download, tag, local, annex_dir, output_file, compute_md5, show_progress, ignore_cached)
    413     if temp_json.exists():
    414         with zipfile.ZipFile(temp_zip, "w", compression=zipfile.ZIP_BZIP2) as simulations_zip:
--> 415             simulations_zip.write(temp_json, arcname="simulations.json")
    416         temp_zip.replace(cache_path)
    417 finally:

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/zipfile.py:1828, in ZipFile.write(self, filename, arcname, compress_type, compresslevel)
   1825     zinfo._compresslevel = self.compresslevel
   1827 with open(filename, "rb") as src, self.open(zinfo, 'w') as dest:
-> 1828     shutil.copyfileobj(src, dest, 1024*8)

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/shutil.py:200, in copyfileobj(fsrc, fdst, length)
    198 if not buf:
    199     break
--> 200 fdst_write(buf)

File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/zipfile.py:1194, in _ZipWriteFile.write(self, data)
   1192 self._crc = crc32(data, self._crc)
   1193 if self._compressor:
-> 1194     data = self._compressor.compress(data)
   1195     self._compress_size += len(data)
   1196 self._fileobj.write(data)

KeyboardInterrupt:

We now use the modes to compute the energy and angular momentum radiated as a function of time. To do so we use the ej_from_hlm method of the Waveform class.

Note that this method requires knowledge of the initial ADM mass and angular momentum of the system, which are stored as attributes of the Waveform object (and automatically loaded, for most catalogs).

M_adm_0 = sxs_waveform.metadata['E0byM']
J_adm_0 = sxs_waveform.metadata['Jz0']
m1 = sxs_waveform.metadata['m1']
m2 = sxs_waveform.metadata['m2']

eb, e, jorb = sxs_waveform.ej_from_hlm(M_adm_0, J_adm_0, m1, m2, modes=modes)

# identify merger
e_mrg = eb[mrg_idx]
j_mrg = jorb[mrg_idx]

fig, ax = plt.subplots()
ax.plot(jorb, eb)
ax.set_xlabel(r'$J/M^2$')
ax.set_ylabel(r'$E_b/M$')
ax.scatter(j_mrg, e_mrg, color='k', label='Merger')
ax.set_title('Binding energy vs Angular Momentum')
plt.show()
../_images/74c71ba965d60e521d04ac2ce80c034212f53eb5bfadbddbe17cf21441ad38b3.png

We will now look at the dynamics. This is not available for all catalogs/models. In the case of SXS waveforms, we do load Horizon data when we download the waveform and can access it.

# plot the mass / spin evolution
fig, axs = plt.subplots(2, 1, figsize=(12, 8), sharex=True)
axs[0].semilogy(sxs_waveform.dyn['t'], abs(sxs_waveform.dyn['m1']-m1), label=r'$m_1$')
axs[0].semilogy(sxs_waveform.dyn['t'], abs(sxs_waveform.dyn['m2']-m2), label=r'$m_2$', linestyle='--')
axs[1].semilogy(sxs_waveform.dyn['t'], sxs_waveform.dyn['chi1'], label=r'S_1')
axs[1].semilogy(sxs_waveform.dyn['t'], sxs_waveform.dyn['chi2'], label=r'S_2', linestyle='--')
axs[0].set_ylabel(r'$\Delta m_i$')
axs[1].set_ylabel(r'$|S_i|$')
axs[1].set_xlabel('t/M')
axs[0].legend()
plt.show()


# plot the trajectories
xA = sxs_waveform.dyn['x1']
xB = sxs_waveform.dyn['x2']
fig, ax = plt.subplots(figsize=(8, 8))
ax.plot(xA[:, 0], xA[:, 1], label='BH 1')
ax.plot(xB[:, 0], xB[:, 1], label='BH 2')
ax.set_xlabel('x/M')
ax.set_ylabel('y/M')
ax.set_title('Black Hole Trajectories')
ax.legend()
plt.show()
../_images/dc1bef1a2c90f60b7583b23a64b14b0085135be5c53dd13543ec2ff5c1df4f84.png ../_images/bd2db77485439e71e83a1593f72014469cbccddceb75c73c90fd4b0a32b67418.png