from importlib.resources import files
import numpy as np
import lal
from scipy.interpolate import interp1d
from pycbc.conversions import get_final_from_initial, get_lm_f0tau
from pycbc.waveform import get_td_waveform_modes
from pycbc.types import (TimeSeries, complex64, zeros)
[docs]
def gen_nrsurqnm(**kwds):
hlm = get_td_waveform_modes(approximant='NRSur7dq4', mass1=kwds['mass1'], mass2=kwds['mass2'],
spin1x=kwds['spin1x'], spin1y=kwds['spin1y'], spin1z=kwds['spin1z'],
spin2x=kwds['spin2x'], spin2y=kwds['spin2y'], spin2z=kwds['spin2z'],
distance = kwds['distance'],
delta_t=kwds['delta_t'], f_lower=kwds['f_lower'],
mode_array=['22','21','20','33','32','31','30','44','43','42','41','40'])
h = 0
for l in range(2,5):
for m in range(-1*l, l+1):
if l==4 and abs(m)==4:
continue # skip 44 mode for ringdown treatment later
h_modes = hlm[(l,m)][0] + 1j * hlm[(l,m)][1]
Y_lm = lal.SpinWeightedSphericalHarmonic(kwds['inclination'], np.pi/2 - kwds['coa_phase'], -2, l, m)
h += h_modes * Y_lm
if 'ringdown_mode' not in kwds or kwds['ringdown_mode'] is None:
return h.real(), -h.imag()
h44 = hlm[(4,4)][0] + 1j * hlm[(4,4)][1]
qnm_par = {}
qnm_par['final_mass'], qnm_par['final_spin'] = get_final_from_initial(
kwds['mass1'], kwds['mass2'],
kwds['spin1x'], kwds['spin1y'], kwds['spin1z'],
kwds['spin2x'], kwds['spin2y'], kwds['spin2z'],
approximant='NRSur7dq4',
f_ref=kwds['f_ref'])
qnm_par['freq'] = {}
qnm_par['tau'] = {}
for mode in kwds['ringdown_mode']:
if len(mode) == 3:
l = int(mode[0])
m = int(mode[1])
n = int(mode[2])
qnm_par['freq'][mode], qnm_par['tau'][mode] = get_lm_f0tau(qnm_par['final_mass'],
qnm_par['final_spin'],
l,m,n)
elif len(mode) == 6:
l1 = int(mode[0])
m1 = int(mode[1])
n1 = int(mode[2])
l2 = int(mode[3])
m2 = int(mode[4])
n2 = int(mode[5])
f1, tau1 = get_lm_f0tau(qnm_par['final_mass'], qnm_par['final_spin'], l1,m1,n1)
f2, tau2 = get_lm_f0tau(qnm_par['final_mass'], qnm_par['final_spin'], l2,m2,n2)
qnm_par['freq'][mode] = f1 + f2
qnm_par['tau'][mode] = 1 / (1/tau1 + 1/tau2)
else:
raise ValueError("Invalid mode format in rindown_mode")
#print("QNM parameters:", qnm_par)
t_final = max(qnm_par['tau'][m] for m in kwds['ringdown_mode']) * np.log(1000)
ringdown_start_time = kwds['toffset']
start_idx = int(np.floor(float(ringdown_start_time - h44.start_time) * h44.sample_rate))
end_idx = int(float(ringdown_start_time + t_final - h44.start_time) * h44.sample_rate)
start_time = h44.sample_times[start_idx]
end_time = h44.sample_times[end_idx]
h44_slice = h44.time_slice(start_time, end_time)
N = len(h44_slice)
qnm = {}
for m in kwds['ringdown_mode']:
qnm[m] = TimeSeries(zeros(N, dtype=complex64), delta_t=kwds['delta_t'], epoch = start_time)
sample_times = qnm[m].sample_times.numpy()
omega = 2 * np.pi * qnm_par['freq'][m] - 1j / qnm_par['tau'][m]
qnm[m].data = 1e-22 * np.exp(-1j * omega * (sample_times - start_time) )
M = len(kwds['ringdown_mode'])
G = np.zeros((N, M), dtype=complex)
#print(N,M)
for k, m in enumerate(kwds['ringdown_mode']):
G[:, k] = qnm[m].data
G_H = G.conj().T
A = np.linalg.solve(G_H @ G, G_H @ h44_slice.data)
allqnm = 0
for i, m in enumerate(kwds['ringdown_mode']):
delta_param = 'delta_' + m
if delta_param in kwds and kwds[delta_param] is not None:
#print("Applying deviation", delta_param, kwds[delta_param])
A[i] *= (1 + kwds[delta_param])
qnm[m].data *= A[i]
allqnm += qnm[m]
#print(allqnm.data)
Y_44 = lal.SpinWeightedSphericalHarmonic(kwds['inclination'], np.pi/2 - kwds['coa_phase'], -2, 4, 4)
Y_4m4 = lal.SpinWeightedSphericalHarmonic(kwds['inclination'], np.pi/2 - kwds['coa_phase'], -2, 4, -4)
qnm44 = allqnm * Y_44 + np.conj(allqnm) * Y_4m4
h.data[start_idx:start_idx+len(qnm44)] += qnm44
#print(qnm44.data)
return h.real(), -h.imag()
[docs]
def gen_nrsur7dq4_tdtaper(**kwds):
from pycbc.waveform import get_td_waveform
import pycbc.waveform.utils
if 'approximant' in kwds:
kwds.pop("approximant")
hp, hc = get_td_waveform(approximant='NRSur7dq4', **kwds)
if 'window' not in kwds or kwds['window'] is None:
window = 0.05
else:
window = kwds['window']
hp_taper = pycbc.waveform.utils.td_taper(hp, hp.start_time, hp.start_time + window)
hc_taper = pycbc.waveform.utils.td_taper(hc, hc.start_time, hc.start_time + window)
return hp_taper, hc_taper
[docs]
def qnm_decomposition(qnm_modes, qnm_par, ringdown_start_time, h_target):
'''
Decompose the target waveform h_target into the given QNM modes starting from ringdown_start_time
'''
t_duration = max(qnm_par['tau'][m] for m in qnm_modes) * np.log(1000)
start_idx = int(np.floor(float(ringdown_start_time - h_target.start_time) * h_target.sample_rate))
end_idx = int(float(ringdown_start_time + t_duration - h_target.start_time) * h_target.sample_rate)
start_time = h_target.sample_times[start_idx]
end_time = h_target.sample_times[end_idx]
h_target_slice = h_target.time_slice(start_time, end_time)
N = len(h_target_slice)
qnm = {}
dynamical_range = 1e-22
for m in qnm_modes:
qnm[m] = TimeSeries(zeros(N, dtype=complex64), delta_t=h_target.delta_t, epoch = start_time)
sample_times = qnm[m].sample_times.numpy()
omega = 2 * np.pi * qnm_par['freq'][m] - 1j / qnm_par['tau'][m]
qnm[m].data = dynamical_range * np.exp(-1j * omega * (sample_times - start_time) )
M = len(qnm_modes)
G = np.zeros((N, M), dtype=complex)
for k, m in enumerate(qnm_modes):
G[:, k] = qnm[m].data
G_H = G.conj().T
A = np.linalg.solve(G_H @ G, G_H @ h_target_slice.data)
A_modes = {}
for i, m in enumerate(qnm_modes):
A_modes[m] = A[i]
allqnm = 0
for m in qnm_modes:
#delta_param = 'delta_' + m
#if delta_param in kwds and kwds[delta_param] is not None:
# #print("Applying deviation", delta_param, kwds[delta_param])
# A_modes[m] *= (1 + kwds[delta_param])
qnm[m].data *= A_modes[m]
allqnm += qnm[m]
for m in qnm_modes:
A_modes[m] = A_modes[m] * dynamical_range
omega = 2 * np.pi * qnm_par['freq'][m] - 1j / qnm_par['tau'][m]
A_modes[m] *= np.exp(-1j * omega * (ringdown_start_time - start_time) )
return A_modes, allqnm, h_target_slice
[docs]
def get_qnmpar(qnm_modes, **kwds):
'''Get the QNM parameters (frequency and damping time) for the given modes based on
the remnant black hole properties predicted by NRSur7dq4. The remnant properties are
calculated from the initial binary parameters using NRSur7dq4 fits.
The input kwds should contain the initial binary parameters and f_ref for NRSur7dq4.
The qnm_modes should be a list of mode labels, e.g. ["220","221"] for linear modes,
or ["220220", "220221"] for quadratic modes.
'''
qnm_par = {}
qnm_par['final_mass'], qnm_par['final_spin'] = \
get_final_from_initial(kwds['mass1'], kwds['mass2'],
kwds['spin1x'], kwds['spin1y'], kwds['spin1z'],
kwds['spin2x'], kwds['spin2y'], kwds['spin2z'],
approximant='NRSur7dq4',
f_ref=kwds['f_ref'])
qnm_par['freq'] = {}
qnm_par['tau'] = {}
for mode in qnm_modes:
if len(mode) == 3:
l = int(mode[0])
m = int(mode[1])
n = int(mode[2])
qnm_par['freq'][mode], qnm_par['tau'][mode] = \
get_lm_f0tau(qnm_par['final_mass'], qnm_par['final_spin'],
l,m,n)
elif len(mode) == 6:
l1 = int(mode[0])
m1 = int(mode[1])
n1 = int(mode[2])
l2 = int(mode[3])
m2 = int(mode[4])
n2 = int(mode[5])
f1, tau1 = get_lm_f0tau(qnm_par['final_mass'],qnm_par['final_spin'],l1,m1,n1)
f2, tau2 = get_lm_f0tau(qnm_par['final_mass'],qnm_par['final_spin'],l2,m2,n2)
qnm_par['freq'][mode] = f1 + f2
qnm_par['tau'][mode] = 1 / (1/tau1 + 1/tau2)
else:
raise ValueError("Invalid mode format in rindown_mode")
return qnm_par
QUADRATIC_MODES = ['220220', '220221','221221','220222','221222','222222']
_interpolation_cache = {}
for mode in QUADRATIC_MODES:
filename = files(__package__).joinpath('data', f'spin_{mode}_interpolation.npy')
data = np.load(filename, allow_pickle=True).item()
spin_loaded = data['spin']
complex_ratio_loaded = data[mode+'_complex_ratio']
interp_loaded = interp1d(spin_loaded, complex_ratio_loaded, kind='cubic', fill_value="extrapolate")
_interpolation_cache[mode] = interp_loaded
[docs]
def load_interpolation_function(label):
"""
load interpolation function for quadratic modes from cache
"""
return _interpolation_cache[label]
[docs]
def gen_nrsur_remove_qqnm(**kwds):
'''
Generate a NRSur7dq4 waveform with quadratic modes being removed from (4,4) mode.
The quadratic modes are constructed from (2,2) mode with amplitude ratio from theory predictions.
The ringdown part of (4,4) mode is then replaced by the residual after subtracting the quadratic modes.
kwds should contain:
- mass1, mass2, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z, distance, delta_t, f_lower, f_ref, inclination, coa_phase,
- mode22: the overtone mode used for quadratic mode),
- mode_quadratic: the list of quadratic modes to be removed, e.g. "220220 220221"
- toffset: the start time for ringdown treatment, e.g. 0.002
- quadratic_tgr: the amplitude deviation parameter for quadratic modes, e.g. 0. Optional
- quadratic_tgr_phase: the phase deviation (in radians) for quadratic modes, e.g. 0. Optional
- qqnm_deltaf: the fractional deviation in frequency for the nonGR quadratic modes, e.g. 0. Optional
- qqnm_deltatau: the fractional deviation in damping time for the nonGR quadratic modes, e.g. 0. Optional
'''
hlm = get_td_waveform_modes(approximant='NRSur7dq4',
mass1=kwds['mass1'], mass2=kwds['mass2'],
spin1x=kwds['spin1x'], spin1y=kwds['spin1y'], spin1z=kwds['spin1z'],
spin2x=kwds['spin2x'], spin2y=kwds['spin2y'], spin2z=kwds['spin2z'],
distance = kwds['distance'],
delta_t=kwds['delta_t'],
f_lower=kwds['f_lower'],
f_ref=kwds['f_ref'],
mode_array=['22','21','20','33','32','31','30','44','43','42','41','40'])
h = 0
for l in range(2,5):
for m in range(-1*l, l+1):
if l==4 and abs(m)==4:
continue # skip 44 mode for ringdown treatment later
h_modes = hlm[(l,m)][0] + 1j * hlm[(l,m)][1]
Y_lm = lal.SpinWeightedSphericalHarmonic(kwds['inclination'], np.pi/2 - kwds['coa_phase'], -2, l, m)
h += h_modes * Y_lm
h22 = hlm[(2,2)][0] + 1j * hlm[(2,2)][1]
mode22 = kwds['mode22'].split()
kwds_quadratic_modes = kwds['mode_quadratic'].split()
qnm_par = get_qnmpar(mode22 + QUADRATIC_MODES, **kwds)
# construct QNM from (2,2) mode
qnm_start_time = kwds['toffset']
if qnm_start_time >= 0.002:
A_modes_22, _, _ = qnm_decomposition(mode22, qnm_par, qnm_start_time, h22)
else:
fit_A_modes_22 = {}
fit_A_modes_22['220'] = []
fit_A_modes_22['221'] = []
fit_A_modes_22['222'] = []
for t_fit in np.arange(0.002, 0.00367, 0.0003333):
this_A, _, _ = qnm_decomposition(mode22, qnm_par, t_fit, h22)
# shift to qnm_start_time
for m in ['220','221','222']:
omega = 2 * np.pi * qnm_par['freq'][m] - 1j / qnm_par['tau'][m]
this_A[m] *= np.exp(-1j * omega * (qnm_start_time - t_fit))
fit_A_modes_22[m].append(this_A[m])
A_modes_22 = {}
for m in ['220','221','222']:
real_parts = [z.real for z in fit_A_modes_22[m]]
imag_parts = [z.imag for z in fit_A_modes_22[m]]
real_mean = np.mean(real_parts)
real_std = np.std(real_parts)
imag_mean = np.mean(imag_parts)
imag_std = np.std(imag_parts)
real_samples = np.random.normal(real_mean, real_std, 1)
imag_samples = np.random.normal(imag_mean, imag_std, 1)
A_modes_22[m] = real_samples + 1j * imag_samples
# getting quadratic modes amplitude
A_modes_quadratic = {}
for mode in kwds_quadratic_modes:
load_interp = load_interpolation_function(mode)
chi_f = qnm_par['final_spin']
ratio = load_interp(chi_f)
mode1 = mode[:3]
mode2 = mode[3:]
conversion_factor = kwds['distance'] * 1e6 * lal.PC_SI / qnm_par['final_mass'] / lal.MRSUN_SI
A_modes_quadratic[mode] = A_modes_22[mode1] * A_modes_22[mode2] * ratio * conversion_factor
# subtract quadratic modes from (4,4) mode
h44 = hlm[(4,4)][0] + 1j * hlm[(4,4)][1]
start_idx = int(np.floor(float(qnm_start_time - h44.start_time) * h44.sample_rate))
start_time = h44.sample_times[start_idx]
end_time = h44.sample_times[-1]
h44_slice = h44.time_slice(start_time, end_time)
N = len(h44_slice)
qnm = {}
qnm_nongr_freqtau = {}
# testingGR amplitude
if 'quadratic_tgr' in kwds and kwds['quadratic_tgr'] is not None:
tgr_amplitude = kwds['quadratic_tgr']
else:
tgr_amplitude = 1.0
# testingGR phase (in radians); GR value is 0
if 'quadratic_tgr_phase' in kwds and kwds['quadratic_tgr_phase'] is not None:
tgr_phase = kwds['quadratic_tgr_phase']
else:
tgr_phase = 0.0
# complex deviation factor: amplitude rescaling and phase shift
tgr_factor = tgr_amplitude * np.exp(1j * tgr_phase)
# testingGR frequency and damping time
has_freqtau_tgr = (
('qqnm_deltaf' in kwds and kwds['qqnm_deltaf'] is not None)
or ('qqnm_deltatau' in kwds and kwds['qqnm_deltatau'] is not None)
)
qqnm_deltaf = kwds.get('qqnm_deltaf')
qqnm_deltatau = kwds.get('qqnm_deltatau')
if qqnm_deltaf is None:
qqnm_deltaf = 0.0
if qqnm_deltatau is None:
qqnm_deltatau = 0.0
for m in kwds_quadratic_modes:
# QQNM waveform
qnm[m] = TimeSeries(zeros(N, dtype=complex64), delta_t=h44.delta_t, epoch = start_time)
sample_times = qnm[m].sample_times.numpy()
omega = 2 * np.pi * qnm_par['freq'][m] - 1j / qnm_par['tau'][m]
qnm[m].data = -A_modes_quadratic[m] * np.exp(-1j * omega * (sample_times - qnm_start_time) ) # -1 to account for 'm/2 pi + pi' conversion to NRSur
# Subtract the quadratic mode contribution from (4,4) mode
h44_slice -= tgr_factor * qnm[m]
# Add the quadratic mode with nonGR frequency and damping time
if has_freqtau_tgr:
qnm_nongr_freqtau[m] = TimeSeries(zeros(N, dtype=complex64), delta_t=h44.delta_t, epoch = start_time)
qqnm_nongr_f = qnm_par['freq'][m] * (1 + qqnm_deltaf)
qqnm_nongr_tau = qnm_par['tau'][m] * (1 + qqnm_deltatau)
omega_nongr = 2 * np.pi * qqnm_nongr_f - 1j / qqnm_nongr_tau
qnm_nongr_freqtau[m].data = -A_modes_quadratic[m] * np.exp(-1j * omega_nongr * (sample_times - qnm_start_time) )
h44_slice += qnm_nongr_freqtau[m]
# add the residual (4,4) mode back to the waveform
Y_44 = lal.SpinWeightedSphericalHarmonic(kwds['inclination'], np.pi/2 - kwds['coa_phase'], -2, 4, 4)
Y_4m4 = lal.SpinWeightedSphericalHarmonic(kwds['inclination'], np.pi/2 - kwds['coa_phase'], -2, 4, -4)
pol44 = h44_slice * Y_44 + np.conj(h44_slice) * Y_4m4
h.data[start_idx:start_idx+len(pol44)] += pol44
return h.real(), -h.imag()