added functions to compute lensed waveform

utils/helper_functions.py

+import numpy as np
+from pycbc.waveform import get_fd_waveform
+import pycbc.filter
+from scipy.interpolate import interp1d
+
+def L2(original, prediction):
+    diff = np.sqrt(np.sum(np.abs(original - prediction)**2))
+    norm = np.sqrt(np.sum(np.abs(original)**2))
+    return diff/norm
+
+def ampFac(t, Ft):
+    """
+    Returns the angular frequency (dimensionless)
+    and amplification factor
+    Input params:
+    ---------------
+    t: time array
+    Ft: time series data
+    """
+    fft_ = np.fft.rfft(Ft) * (t[1] - t[0])
+    freq_ = np.fft.rfftfreq(len(t), t[1] - t[0]) * 2.0 * np.pi
+    ampfac = fft_ * freq_ * 1j 
+    ampfac = np.conj(ampfac) * np.exp(1j * freq_ * t[0])
+    ampfac += Ft[-1]
+    return freq_, ampfac
+
+def unlensed_fd_waveform(**params):
+    """
+    Function returning unlensed waveforms
+    """
+    hp, hc = get_fd_waveform(
+        approximant=params["approx"],
+        mass1=params["m1"],
+        mass2=params["m2"],
+        distance=params["d"],
+        spin1z=params["chi1"],
+        spin2z=params["chi2"],
+        inclination=params["iota"],
+        coa_phase=params["phi0"],
+        delta_f=params["delta_f"],
+        f_lower=params["f_low"],
+        f_final=params["f_final"],
+    )
+
+    hp = hp[: int(params["f_final"] / params["delta_f"]) + 1]
+    hc = hc[: int(params["f_final"] / params["delta_f"]) + 1]
+    f = hp.get_sample_frequencies().data
+    return f, hp, hc
+
+def lensed_fd_waveform(**params):
+    """
+    Function which returns lensed waveforms
+    """
+    f, hp, hc = unlensed_fd_waveform(
+        approx=params["approx"],
+        m1=params["m1"],
+        m2=params["m2"],
+        d=params["d"],
+        chi1=params["chi1"],
+        chi2=params["chi2"],
+        iota=params["iota"],
+        phi0=params["phi0"],
+        delta_f=params["delta_f"],
+        f_low=params["f_low"],
+        f_final=params["f_final"],
+    )
+
+    freqs_, ampfac_ = params["f"], params["Ff"]
+    ampfac_abs = np.abs(ampfac_)
+    ampfac_phase = np.angle(ampfac_)
+    Ff_abs_interp = interp1d(
+        freqs_, ampfac_abs, fill_value="extrapolate", kind="linear"
+    )
+    Ff_phase_interp = interp1d(
+        freqs_, np.angle(ampfac_), fill_value="extrapolate", kind="linear"
+    )
+
+    lensed_hp = Ff_abs_interp(f) * np.exp(1j * Ff_phase_interp(f)) * hp
+    lensed_hp = lensed_hp.cyclic_time_shift(3)
+
+    lensed_hc = Ff_abs_interp(f) * np.exp(1j * Ff_phase_interp(f)) * hc
+    lensed_hc = lensed_hc.cyclic_time_shift(3)
+    return f, lensed_hp, lensed_hc
+
+def CalcLogMisMatch(xf1, xf2, Sh, f_min, f_max):
+    match, shift_idx = pycbc.filter.matchedfilter.match(
+        xf1,
+        xf2,
+        psd=Sh,
+        low_frequency_cutoff=f_min,
+        high_frequency_cutoff=f_max,
+        v1_norm=None,
+        v2_norm=None,
+    )
+    if match >= 1:
+        return -100
+    else:
+        return np.log10(1 - match)
\ No newline at end of file