Working with ICCUB Catalog and Waveform Integration¶
This tutorial demonstrates how to work with waveforms from the ICCUB (Institute of Cosmos Sciences) numerical relativity catalog and perform waveform integration from Ψ₄ (Weyl scalar) to strain h.
ICCUB simulations often provide Ψ₄ data, which needs to be integrated twice to obtain the gravitational wave strain. PyART provides tools for this integration.
Setup¶
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
import numpy as np
from PyART.catalogs.icc import Waveform_ICC
from PyART.analysis.integrate_multipole import Multipole
from PyART.analysis.integrate_wave import IntegrateMultipole
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 7
5 import numpy as np
6 from PyART.catalogs.icc import Waveform_ICC
----> 7 from PyART.analysis.integrate_multipole import Multipole
8 from PyART.analysis.integrate_wave import IntegrateMultipole
ModuleNotFoundError: No module named 'PyART.analysis.integrate_multipole'
Load ICCUB Waveform¶
First, we load a waveform from the ICCUB catalog. We’ll load it in two ways:
Without integration (raw Ψ₄ data)
With integration (computed strain)
# Path to ICCUB catalog data
# Note: Adjust this path to your local data directory
path = './local_data/icc/catalog'
ID = 1 # Simulation ID
# Integration options
integr_opts = {
'f0': 0.002, # Low-frequency cutoff
'extrap_psi4': True, # Extrapolate Psi4 to infinity
'method': 'FFI', # Fixed Frequency Integration method
'window': [20, -20] # Time window for integration
}
# Load without integration (raw data)
print("Loading raw Psi4 data...")
icc = Waveform_ICC(path=path, ID=ID)
# Load with integration
print("Loading and integrating to strain...")
iicc = Waveform_ICC(
path=path,
ID=ID,
integrate=True,
integr_opts=integr_opts,
load_puncts=True # Also load puncture trajectory data
)
print("\nWaveform loaded successfully!")
Display Metadata¶
Let’s examine the simulation parameters:
print("Simulation Metadata:")
print('=' * 60)
for k, v in icc.metadata.items():
print(f'{k:20s} : {v}')
print('=' * 60)
Compare Stored vs Integrated Waveforms¶
If the simulation provides both Ψ₄ and pre-computed strain, we can compare with our integration:
# Find merger times
try:
iicc_tmrg, _, _, _ = iicc.find_max(height=0.3)
icc_tmrg, _, _, _ = icc.find_max(height=0.3)
except:
iicc_tmrg = 0
icc_tmrg = 0
# Define waveform types to plot
waves = ['hlm', 'dothlm', 'psi4lm']
ylabs = [r'$h_{22}$', r'$\dot{h}_{22}$', r'$\psi_4^{22}$']
npanels = len(waves)
plt.figure(figsize=(12, 3*npanels))
plt.suptitle(f'ICCUB Simulation ID: {ID:04d}', fontsize=18, y=0.995)
for i, wave_name in enumerate(waves):
plt.subplot(npanels, 1, i+1)
# Get integrated waveform time
if wave_name == 'psi4lm':
iicc_t = iicc.t_psi4
else:
iicc_t = iicc.u - iicc_tmrg
# Plot stored waveform if available
icc_wave = getattr(icc, wave_name)
if icc_wave:
if wave_name == 'psi4lm':
icc_t = icc.t_psi4
else:
icc_t = icc.u - icc_tmrg
plt.plot(icc_t, icc_wave[(2,2)]['A'],
lw=1.5, c='red', alpha=0.7, label='Stored')
plt.plot(icc_t, icc_wave[(2,2)]['real'],
lw=1.0, c='red', alpha=0.5)
plt.plot(icc_t, icc_wave[(2,2)]['imag'],
lw=1.0, c='red', ls=':', alpha=0.5)
# Plot integrated waveform
iicc_wave = getattr(iicc, wave_name)
plt.plot(iicc_t, iicc_wave[(2,2)]['A'],
lw=1.5, c='blue', alpha=0.7, label='Integrated')
plt.plot(iicc_t, iicc_wave[(2,2)]['real'],
lw=1.0, c='blue', alpha=0.5)
plt.plot(iicc_t, iicc_wave[(2,2)]['imag'],
lw=1.0, c='blue', ls=':', alpha=0.5)
plt.ylabel(ylabs[i], fontsize=14)
plt.xlabel(r'$t-t_{\rm mrg}$ (M)', fontsize=14) if i == len(waves)-1 else None
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Visualize Puncture Trajectories¶
For scattering simulations, we can plot the trajectories of the two black holes:
if iicc.metadata['scat_angle'] is not None:
plt.figure(figsize=(10, 8))
# Plot both puncture trajectories
plt.plot(iicc.puncts['x0'], iicc.puncts['y0'],
linewidth=2, label='BH 1', marker='o', markersize=3)
plt.plot(iicc.puncts['x1'], iicc.puncts['y1'],
linewidth=2, label='BH 2', marker='s', markersize=3)
plt.xlabel(r'$x$ (M)', fontsize=16)
plt.ylabel(r'$y$ (M)', fontsize=16)
plt.title(f'Black Hole Trajectories (Scattering Angle: {iicc.metadata["scat_angle"]:.2f}°)',
fontsize=16)
plt.legend(fontsize=14)
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.tight_layout()
plt.show()
else:
print("This is not a scattering simulation - no trajectory plot available.")
Testing the Integration Method¶
We can explicitly test the integration by comparing old and new integration methods:
# Get Psi4 data
t = icc.t_psi4
psi4 = icc.psi4lm[(2,2)]['z']
l, m = 2, 2
M = 1
# New integration method
mode_new = IntegrateMultipole(
l, m, t, psi4,
**integr_opts,
mass=M,
radius=icc.r_extr,
integrand='psi4'
)
# Old integration method
mode_old = Multipole(l, m, t, psi4, mass=M, radius=icc.r_extr, integrand='psi4')
mode_old.hlm, mode_old.dothlm = mode_old.integrate_wave(integr_opts=integr_opts)
# Plot comparison
plt.figure(figsize=(14, 10))
# Strain h
plt.subplot(3, 2, 1)
plt.plot(t, mode_old.hlm.real, c='r', ls='-', label='Old method', linewidth=2)
plt.plot(t, mode_new.h.real, c='b', ls='--', label='New method', linewidth=2, alpha=0.7)
plt.ylabel(r'Re[$h_{22}$]', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.title('Strain Comparison', fontsize=14)
plt.subplot(3, 2, 2)
plt.plot(t, np.abs(mode_old.hlm - mode_new.h), c='green', label='Difference')
plt.ylabel(r'$|h_{\rm old} - h_{\rm new}|$', fontsize=14)
plt.yscale('log')
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.title('Absolute Difference', fontsize=14)
# Strain derivative
plt.subplot(3, 2, 3)
plt.plot(t, mode_old.dothlm.real, c='r', ls='-', linewidth=2)
plt.plot(t, mode_new.dh.real, c='b', ls='--', linewidth=2, alpha=0.7)
plt.ylabel(r'Re[$\dot{h}_{22}$]', fontsize=14)
plt.grid(True, alpha=0.3)
plt.title('Strain Derivative', fontsize=14)
plt.subplot(3, 2, 4)
plt.plot(t, np.abs(mode_old.dothlm - mode_new.dh), c='green')
plt.ylabel(r'$|\dot{h}_{\rm old} - \dot{h}_{\rm new}|$', fontsize=14)
plt.grid(True, alpha=0.3)
plt.title('Derivative Difference', fontsize=14)
# Psi4
plt.subplot(3, 2, 5)
plt.plot(t, mode_old.psi.real, c='r', ls='-', linewidth=2)
plt.plot(t, mode_new.psi4.real, c='b', ls='--', linewidth=2, alpha=0.7)
plt.xlabel('Time (M)', fontsize=14)
plt.ylabel(r'Re[$\psi_4^{22}$]', fontsize=14)
plt.grid(True, alpha=0.3)
plt.title('Weyl Scalar', fontsize=14)
plt.subplot(3, 2, 6)
plt.plot(t, np.abs(mode_old.psi - mode_new.psi4), c='green')
plt.xlabel('Time (M)', fontsize=14)
plt.ylabel(r'$|\psi_4^{\rm old} - \psi_4^{\rm new}|$', fontsize=14)
plt.grid(True, alpha=0.3)
plt.title('Psi4 Difference', fontsize=14)
plt.tight_layout()
plt.show()
print(f"\nMaximum relative difference in h: {np.max(np.abs(mode_old.hlm - mode_new.h)/np.max(np.abs(mode_old.hlm))):.2e}")
Summary¶
This tutorial demonstrated:
Loading waveforms from the ICCUB catalog
Integrating Ψ₄ (Weyl scalar) to obtain gravitational wave strain
Comparing stored and integrated waveforms
Visualizing puncture trajectories for scattering simulations
Testing and validating integration methods
Key Points¶
ICCUB provides high-quality numerical relativity data including scattering scenarios
The Fixed Frequency Integration (FFI) method provides accurate strain from Ψ₄
Integration requires careful choices of f₀ (low-frequency cutoff) and time windows
PyART’s integration tools handle extrapolation and windowing automatically
Next Steps¶
Explore different integration methods
Compute scattering angles from trajectory data
Compare ICCUB waveforms with other catalogs
Perform mismatch calculations with EOB models