"import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from scipy.integrate import solve_ivp
from scipy import signal

# ============================================================
#  CONSTANTS  (SI)
# ============================================================
C      = 2.99792458e8       # m/s
G      = 6.67430e-11        # m³ kg⁻¹ s⁻²
HBAR   = 1.054571817e-34    # J s
KB     = 1.380649e-23       # J/K
L_PL   = np.sqrt(HBAR * G / C**3)  # Planck length ~1.616e-35 m

# ============================================================
#  UNIVERSAL SNR²  —  matched filter for a known waveform
# ============================================================
def snr_squared(h_t, dt, psd_func):
    """
    Matched-filter SNR² for a deterministic waveform h(t).

    Parameters
    ----------
    h_t : array_like
        Time-domain signal samples.
    dt : float
        Sample spacing (s).
    psd_func : callable
        Function returning one-sided PSD S_n(f) [strain²/Hz].

    Returns
    -------
    snr2 : float
        4 ∫ |h̃(f)|² / S_n(f) df
    """
    h = np.asarray(h_t, dtype=float)
    h = h - np.mean(h)
    N = len(h)
    H = np.fft.rfft(h) * dt          # Fourier transform (integral, not sum)
    f = np.fft.rfftfreq(N, dt)
    S_n = psd_func(f)
    S_n = np.maximum(S_n, 1e-300)    # avoid division by zero
    # One-sided PSD → integral over positive frequencies only
    integrand = np.abs(H)**2 / S_n
    return 4.0 * np.trapz(integrand, f)"