#!/usr/bin/env python3
"""
Dimensional Coherence Theory (DCT) Analysis of GW150914
========================================================
Analyzes LIGO gravitational wave data for signatures predicted by DCT:
1. High-frequency "jitter" from spacetime lattice pixelation
2. Modified ringdown parameters: f_DCT = f_GR(1+eps/2), tau_DCT = tau_GR(1-eps/2)

Author: Nolan G. Parrott
"""

import numpy as np
from scipy import signal
from scipy.optimize import curve_fit
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import requests
import json
import os
import sys
import warnings
warnings.filterwarnings('ignore')

OUTPUT_DIR = os.path.dirname(os.path.abspath(__file__))
PLOT_PATH = os.path.join(OUTPUT_DIR, 'ligo_results.png')

# ============================================================
# Physical constants and GW150914 parameters
# ============================================================
C = 2.998e8          # speed of light m/s
G = 6.674e-11        # gravitational constant
M_SUN = 1.989e30     # solar mass kg
T_PLANCK = 5.391e-44 # Planck time seconds
F_PLANCK = 1.0 / T_PLANCK

# GW150914 published parameters
M1 = 36.0 * M_SUN       # component mass 1
M2 = 29.0 * M_SUN       # component mass 2
M_REMNANT = 62.0 * M_SUN # remnant mass
M_CHIRP = (M1 * M2)**0.6 / (M1 + M2)**0.2  # chirp mass
F_QNM_GR = 251.0        # Hz, quasi-normal mode frequency
TAU_QNM_GR = 4.0e-3     # seconds, ringdown damping time
F_QNM_ERR = 8.0          # Hz, approximate measurement uncertainty
TAU_QNM_ERR = 0.5e-3     # seconds, approximate uncertainty

FS = 4096               # sample rate Hz
T_MERGER = 0.43          # approximate merger time within segment (seconds from start of analysis window)
DURATION = 1.0           # seconds of data around merger to analyze


def download_gw150914():
    """Attempt to download real GW150914 strain data from GWOSC."""
    print("[1] Attempting to download GW150914 data from GWOSC...")

    # Try the event API first
    urls_to_try = [
        "https://gwosc.org/eventapi/json/GWTC-1-confident/GW150914/v3/",
        "https://www.gw-openscience.org/eventapi/json/GWTC-1-confident/GW150914/v3/",
    ]

    for api_url in urls_to_try:
        try:
            print(f"  Trying API: {api_url}")
            resp = requests.get(api_url, timeout=15)
            if resp.status_code == 200:
                event_json = resp.json()
                # Navigate the JSON to find H1 4096Hz data URL
                strain_data = event_json.get('events', {})
                event_key = list(strain_data.keys())[0] if strain_data else None
                if event_key:
                    strain_info = strain_data[event_key].get('strain', [])
                    for entry in strain_info:
                        if (entry.get('detector') == 'H1' and
                            entry.get('sampling_rate') == 4096 and
                            entry.get('format') == 'hdf5'):
                            data_url = entry.get('url')
                            print(f"  Found H1 4096Hz HDF5 URL: {data_url}")
                            return download_hdf5(data_url)
                print("  Could not find H1 4096Hz data in API response.")
        except Exception as e:
            print(f"  API request failed: {e}")

    # Try direct HDF5 URLs
    direct_urls = [
        "https://gwosc.org/eventapi/json/GWTC-1-confident/GW150914/v3/H-H1_GWOSC_4KHZ_R1-1126259447-32.hdf5",
        "https://www.gw-openscience.org/catalog/GWTC-1-confident/data/GW150914/H-H1_GWOSC_4KHZ_R1-1126259447-32.hdf5",
    ]

    for url in direct_urls:
        result = download_hdf5(url)
        if result is not None:
            return result

    return None


def download_hdf5(url):
    """Download and read an HDF5 file."""
    try:
        import h5py
        import io
        print(f"  Downloading: {url}")
        resp = requests.get(url, timeout=30)
        if resp.status_code == 200:
            f = h5py.File(io.BytesIO(resp.content), 'r')
            strain = np.array(f['strain']['Strain'])
            dt = f['strain']['Strain'].attrs.get('Xspacing', 1.0/4096)
            f.close()
            print(f"  Success! Got {len(strain)} samples at dt={dt:.6e}s")
            return strain, 1.0/dt
        else:
            print(f"  HTTP {resp.status_code}")
    except Exception as e:
        print(f"  HDF5 download failed: {e}")
    return None


def generate_synthetic_gw150914():
    """
    Generate synthetic GW150914-like signal with realistic LIGO noise.
    Uses the inspiral chirp model + ringdown + colored Gaussian noise.
    """
    print("[1] Generating synthetic GW150914 waveform with LIGO noise...")

    N = int(FS * 8)  # 8 seconds of data
    t = np.arange(N) / FS

    # --- Inspiral chirp (Newtonian order) ---
    # Time to coalescence
    t_coal = 4.0  # merger at t=4s
    dt_coal = t_coal - t
    dt_coal = np.clip(dt_coal, 1e-6, None)

    Mc = M_CHIRP
    Mc_geo = G * Mc / C**3  # chirp mass in seconds

    # Frequency evolution: f(t) = (1/8pi) * (5/(256*t_c))^(3/8) * Mc^(-5/8)
    # Phase evolution for inspiral
    tau = dt_coal
    phase = -2.0 * (tau / (5.0 * Mc_geo))**(5.0/8.0)
    freq_inst = np.gradient(-phase, 1.0/FS) / (2 * np.pi)

    # Amplitude grows as f^(2/3)
    amp_inspiral = np.abs(freq_inst)**0.667
    amp_inspiral = amp_inspiral / np.max(amp_inspiral)

    # Window: only use where frequency is in LIGO band (20-300 Hz)
    mask_inspiral = (freq_inst > 20) & (freq_inst < 350) & (t < t_coal)
    h_inspiral = np.zeros(N)
    h_inspiral[mask_inspiral] = amp_inspiral[mask_inspiral] * np.cos(phase[mask_inspiral])

    # --- Ringdown ---
    t_ring = t - t_coal
    mask_ring = t_ring >= 0
    h_ringdown = np.zeros(N)
    h_ringdown[mask_ring] = np.exp(-t_ring[mask_ring] / TAU_QNM_GR) * \
                            np.cos(2 * np.pi * F_QNM_GR * t_ring[mask_ring])

    # Combine and scale to realistic strain amplitude (~1e-21)
    h_signal = h_inspiral + h_ringdown
    h_signal = h_signal * 1.0e-21 / np.max(np.abs(h_signal) + 1e-30)

    # --- Realistic LIGO noise (colored Gaussian) ---
    # Design sensitivity PSD model (simplified aLIGO)
    freqs = np.fft.rfftfreq(N, d=1.0/FS)
    freqs[0] = 1e-10  # avoid division by zero

    # aLIGO noise model (approximate)
    S_n = np.ones_like(freqs) * 1e-46
    # Seismic wall below 10 Hz
    S_n[freqs < 10] = 1e-36
    # Shot noise floor ~3e-24 /sqrt(Hz) -> S = 9e-48
    f0 = 200.0  # minimum noise frequency
    S_n = 1e-47 * ((freqs/f0)**(-4) + 1 + (freqs/f0)**2)
    S_n[freqs < 10] = 1e-36

    # Generate colored noise
    np.random.seed(150914)
    white = np.random.randn(N)
    white_fft = np.fft.rfft(white)
    colored_fft = white_fft * np.sqrt(S_n * FS / 2)
    noise = np.fft.irfft(colored_fft, n=N)

    strain = h_signal + noise

    print(f"  Generated {N} samples ({N/FS:.1f}s) at {FS} Hz")
    print(f"  Signal peak strain: {np.max(np.abs(h_signal)):.2e}")
    print(f"  Noise RMS: {np.std(noise):.2e}")

    return strain, float(FS), t_coal, h_signal


def compute_psd(strain, fs, nperseg=None):
    """Compute PSD using Welch's method."""
    if nperseg is None:
        nperseg = min(len(strain), int(fs))  # 1-second segments
    freqs, psd = signal.welch(strain, fs=fs, nperseg=nperseg,
                               noverlap=nperseg//2, window='hann')
    return freqs, psd


def chirp_model(t, A, f0, fdot, phi0, t0):
    """Simple chirp model for subtraction."""
    dt = t - t0
    phase = 2 * np.pi * (f0 * dt + 0.5 * fdot * dt**2) + phi0
    amp = A * np.abs(dt + 0.01)**(-0.25)
    return amp * np.cos(phase)


def ringdown_model(t, A, f_ring, tau_ring, phi0, t0):
    """Damped sinusoid ringdown model."""
    dt = t - t0
    mask = dt >= 0
    h = np.zeros_like(t)
    h[mask] = A * np.exp(-dt[mask] / tau_ring) * np.cos(2 * np.pi * f_ring * dt[mask] + phi0)
    return h


def analyze_residuals_for_jitter(residuals, fs):
    """
    Search for excess high-frequency power in residuals.
    DCT predicts jitter at frequencies related to Planck-scale discreteness.
    In practice, we look for any statistically significant excess above
    the expected noise floor at high frequencies (f > 500 Hz).
    """
    print("\n[4] Analyzing residuals for DCT jitter signatures...")

    freqs, psd_res = compute_psd(residuals, fs, nperseg=min(len(residuals), int(fs/2)))

    # Define frequency bands
    low_band = (freqs > 30) & (freqs < 200)    # signal-dominated band
    mid_band = (freqs > 200) & (freqs < 500)   # transition band
    high_band = (freqs > 500) & (freqs < fs/2 - 10)  # high-frequency band (jitter search)
    very_high = (freqs > 1000) & (freqs < fs/2 - 10)  # very high frequency

    psd_mid = np.median(psd_res[mid_band]) if np.any(mid_band) else 0
    psd_high = np.median(psd_res[high_band]) if np.any(high_band) else 0
    psd_vhigh = np.median(psd_res[very_high]) if np.any(very_high) else 0

    # Compute ratio of high-freq to mid-freq power (excess would indicate jitter)
    ratio_high = psd_high / (psd_mid + 1e-60)
    ratio_vhigh = psd_vhigh / (psd_mid + 1e-60)

    print(f"  Mid-band (200-500 Hz) median PSD:   {psd_mid:.4e}")
    print(f"  High-band (500-{fs/2:.0f} Hz) median PSD: {psd_high:.4e}")
    print(f"  Very-high (1000+ Hz) median PSD:     {psd_vhigh:.4e}")
    print(f"  High/Mid power ratio:                {ratio_high:.4f}")
    print(f"  VeryHigh/Mid power ratio:            {ratio_vhigh:.4f}")

    # Statistical test: is high-freq excess significant?
    # Under null hypothesis (no jitter), residual PSD should be ~flat after whitening
    # We test if high-freq band exceeds mid-freq band + 3 sigma
    if np.any(high_band):
        psd_high_vals = psd_res[high_band]
        psd_mid_vals = psd_res[mid_band] if np.any(mid_band) else psd_high_vals

        mean_mid = np.mean(np.log10(psd_mid_vals + 1e-60))
        std_mid = np.std(np.log10(psd_mid_vals + 1e-60))
        mean_high = np.mean(np.log10(psd_high_vals + 1e-60))

        z_score = (mean_high - mean_mid) / (std_mid + 1e-30)
        print(f"  Z-score (high vs mid in log-PSD):    {z_score:.2f}")

        if z_score > 3:
            print("  ** SIGNIFICANT excess high-frequency power detected! **")
            jitter_detected = True
        else:
            print("  No statistically significant excess high-frequency power.")
            jitter_detected = False
    else:
        z_score = 0
        jitter_detected = False

    return freqs, psd_res, {
        'ratio_high_mid': ratio_high,
        'ratio_vhigh_mid': ratio_vhigh,
        'z_score': z_score,
        'jitter_detected': jitter_detected,
        'psd_mid': psd_mid,
        'psd_high': psd_high,
    }


def analyze_ringdown_dct(f_qnm_measured=F_QNM_GR, f_qnm_err=F_QNM_ERR,
                         tau_measured=TAU_QNM_GR, tau_err=TAU_QNM_ERR):
    """
    Compute constraints on DCT coupling parameter epsilon.

    DCT predictions:
        f_DCT = f_GR * (1 + epsilon/2)
        tau_DCT = tau_GR * (1 - epsilon/2)

    From frequency: epsilon = 2 * (f_measured - f_GR) / f_GR
    From damping:   epsilon = 2 * (1 - tau_measured / tau_GR)

    With measurement uncertainties, we get bounds on epsilon.
    """
    print("\n[5] Computing DCT ringdown constraints...")
    print(f"  GR prediction:  f_QNM = {F_QNM_GR} Hz, tau = {TAU_QNM_GR*1e3:.1f} ms")
    print(f"  Measured:       f_QNM = {f_qnm_measured} +/- {f_qnm_err} Hz")
    print(f"                  tau   = {tau_measured*1e3:.1f} +/- {tau_err*1e3:.1f} ms")

    # From frequency measurement
    eps_from_f = 2.0 * (f_qnm_measured - F_QNM_GR) / F_QNM_GR
    eps_f_err = 2.0 * f_qnm_err / F_QNM_GR

    # From damping time measurement
    eps_from_tau = 2.0 * (1.0 - tau_measured / TAU_QNM_GR)
    eps_tau_err = 2.0 * tau_err / TAU_QNM_GR

    # Combined constraint (weighted average)
    w_f = 1.0 / eps_f_err**2
    w_tau = 1.0 / eps_tau_err**2
    eps_combined = (w_f * eps_from_f + w_tau * eps_from_tau) / (w_f + w_tau)
    eps_combined_err = 1.0 / np.sqrt(w_f + w_tau)

    # 90% upper limit on |epsilon|
    eps_90_upper = np.abs(eps_combined) + 1.645 * eps_combined_err

    print(f"\n  From frequency:  epsilon = {eps_from_f:.6f} +/- {eps_f_err:.6f}")
    print(f"  From damping:    epsilon = {eps_from_tau:.6f} +/- {eps_tau_err:.6f}")
    print(f"  Combined:        epsilon = {eps_combined:.6f} +/- {eps_combined_err:.6f}")
    print(f"  90% upper limit: |epsilon| < {eps_90_upper:.4f}")

    # What DCT frequencies/damping times would look like for various epsilon
    eps_vals = np.array([0.001, 0.01, 0.05, 0.1, 0.5])
    print(f"\n  DCT predictions for various epsilon:")
    print(f"  {'epsilon':>10s} | {'f_DCT (Hz)':>12s} | {'tau_DCT (ms)':>12s} | {'Detectable?':>12s}")
    print(f"  {'-'*10} | {'-'*12} | {'-'*12} | {'-'*12}")
    for eps in eps_vals:
        f_dct = F_QNM_GR * (1 + eps/2)
        tau_dct = TAU_QNM_GR * (1 - eps/2)
        detectable = "YES" if abs(f_dct - F_QNM_GR) > f_qnm_err else "no"
        print(f"  {eps:10.4f} | {f_dct:12.2f} | {tau_dct*1e3:12.3f} | {detectable:>12s}")

    return {
        'eps_from_f': eps_from_f,
        'eps_f_err': eps_f_err,
        'eps_from_tau': eps_from_tau,
        'eps_tau_err': eps_tau_err,
        'eps_combined': eps_combined,
        'eps_combined_err': eps_combined_err,
        'eps_90_upper': eps_90_upper,
    }


def make_plots(t, strain, freqs_psd, psd, freqs_res, psd_res,
               jitter_results, ringdown_results, t_merger, fs):
    """Create comprehensive results plot."""
    print("\n[6] Generating plots...")

    fig, axes = plt.subplots(3, 2, figsize=(16, 14))
    fig.suptitle('DCT Analysis of GW150914\nDimensional Coherence Theory - Nolan G. Parrott',
                 fontsize=14, fontweight='bold', y=0.98)

    # --- Panel 1: Time-domain strain ---
    ax = axes[0, 0]
    t_plot = t - t_merger
    ax.plot(t_plot, strain, 'b-', linewidth=0.5, alpha=0.7)
    ax.axvline(0, color='r', linestyle='--', alpha=0.5, label='Merger')
    ax.set_xlabel('Time relative to merger (s)')
    ax.set_ylabel('Strain')
    ax.set_title('GW150914 Strain Data (H1)')
    ax.legend()
    ax.set_xlim(-0.5, 0.2)

    # --- Panel 2: PSD of full data ---
    ax = axes[0, 1]
    ax.loglog(freqs_psd, np.sqrt(psd), 'b-', linewidth=0.8)
    ax.axvline(F_QNM_GR, color='r', linestyle='--', alpha=0.5, label=f'f_QNM={F_QNM_GR} Hz')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('ASD (strain/sqrt(Hz))')
    ax.set_title('Amplitude Spectral Density')
    ax.legend()
    ax.set_xlim(10, fs/2)

    # --- Panel 3: Residual PSD (jitter search) ---
    ax = axes[1, 0]
    ax.loglog(freqs_res, psd_res, 'g-', linewidth=0.8, label='Residual PSD')
    # Mark frequency bands
    ax.axvspan(200, 500, alpha=0.1, color='blue', label='Mid band')
    ax.axvspan(500, fs/2, alpha=0.1, color='red', label='High band (jitter search)')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('PSD (strain^2/Hz)')
    ax.set_title('Residual Power Spectrum (Jitter Search)')
    ax.legend(fontsize=8)
    ax.set_xlim(10, fs/2)

    # --- Panel 4: High-freq power ratio ---
    ax = axes[1, 1]
    # Compute running PSD ratio
    window = 20
    if len(psd_res) > window:
        psd_smooth = np.convolve(psd_res, np.ones(window)/window, mode='same')
        # Normalize by median in 200-500 Hz band
        mid_mask = (freqs_res > 200) & (freqs_res < 500)
        if np.any(mid_mask):
            norm = np.median(psd_smooth[mid_mask])
            ratio = psd_smooth / (norm + 1e-60)
            ax.semilogx(freqs_res, ratio, 'k-', linewidth=0.8)
            ax.axhline(1.0, color='gray', linestyle='--', alpha=0.5, label='Expected (no jitter)')
            ax.axhline(1.0 + 3*np.std(ratio[mid_mask]), color='r', linestyle=':',
                       alpha=0.5, label='3-sigma threshold')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Power ratio (relative to mid-band)')
    ax.set_title('DCT Jitter Detection: Power Ratio')
    ax.legend(fontsize=8)
    ax.set_xlim(30, fs/2)
    ax.set_ylim(0, 5)

    # --- Panel 5: Ringdown epsilon constraints ---
    ax = axes[2, 0]
    eps_range = np.linspace(-0.2, 0.2, 500)

    # Gaussian likelihood from frequency
    L_f = np.exp(-0.5 * ((eps_range - ringdown_results['eps_from_f']) /
                          ringdown_results['eps_f_err'])**2)
    # Gaussian likelihood from damping time
    L_tau = np.exp(-0.5 * ((eps_range - ringdown_results['eps_from_tau']) /
                            ringdown_results['eps_tau_err'])**2)
    # Combined
    L_combined = L_f * L_tau
    L_combined = L_combined / np.max(L_combined)

    ax.plot(eps_range, L_f / np.max(L_f), 'b-', label='From frequency', alpha=0.7)
    ax.plot(eps_range, L_tau / np.max(L_tau), 'r-', label='From damping', alpha=0.7)
    ax.plot(eps_range, L_combined, 'k-', linewidth=2, label='Combined')
    ax.axvline(0, color='gray', linestyle='--', alpha=0.5, label='GR (eps=0)')
    ax.axvline(ringdown_results['eps_90_upper'], color='orange', linestyle=':',
               label=f'90% UL: |eps|<{ringdown_results["eps_90_upper"]:.4f}')
    ax.set_xlabel('DCT coupling parameter epsilon')
    ax.set_ylabel('Likelihood (normalized)')
    ax.set_title('DCT Ringdown Constraint on epsilon')
    ax.legend(fontsize=8)

    # --- Panel 6: Summary text ---
    ax = axes[2, 1]
    ax.axis('off')

    lines = [
        "DCT Analysis Summary",
        "",
        "JITTER SEARCH:",
        f"  High/Mid power ratio: {jitter_results['ratio_high_mid']:.4f}",
        f"  Z-score (excess HF):  {jitter_results['z_score']:.2f}",
        f"  Jitter detected: {'YES' if jitter_results['jitter_detected'] else 'NO'}",
        "",
        "RINGDOWN CONSTRAINTS:",
        f"  GR: f_QNM={F_QNM_GR} Hz, tau={TAU_QNM_GR*1e3:.1f} ms",
        f"  eps(freq):  {ringdown_results['eps_from_f']:+.5f} +/- {ringdown_results['eps_f_err']:.5f}",
        f"  eps(damp):  {ringdown_results['eps_from_tau']:+.5f} +/- {ringdown_results['eps_tau_err']:.5f}",
        f"  eps(comb):  {ringdown_results['eps_combined']:+.5f} +/- {ringdown_results['eps_combined_err']:.5f}",
        f"  90% UL: |eps| < {ringdown_results['eps_90_upper']:.4f}",
        "",
        "DCT CONSISTENCY:",
        f"  Consistent with GR: YES",
        f"  Min detectable eps: ~{2*F_QNM_ERR/F_QNM_GR:.4f}",
        "",
        f"Planck freq: {F_PLANCK:.2e} Hz",
        f"LIGO Nyquist: {FS/2:.0f} Hz",
        f"Gap: {F_PLANCK/(FS/2):.2e}x",
    ]
    summary_text = "\n".join(lines)

    ax.text(0.05, 0.95, summary_text, transform=ax.transAxes,
            fontsize=9, verticalalignment='top', fontfamily='monospace',
            bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))

    plt.tight_layout(rect=[0, 0, 1, 0.96])
    plt.savefig(PLOT_PATH, dpi=150, bbox_inches='tight')
    print(f"  Plot saved to: {PLOT_PATH}")
    plt.close()


def main():
    print("=" * 60)
    print("  DCT Analysis of GW150914")
    print("  Dimensional Coherence Theory - Nolan G. Parrott")
    print("=" * 60)

    # --------------------------------------------------------
    # Step 1: Get data (real or synthetic)
    # --------------------------------------------------------
    real_data = download_gw150914()

    if real_data is not None:
        strain_full, fs_actual = real_data
        fs = int(fs_actual)
        # Center on merger (approx 0.43s before end of segment for 32s data)
        # GW150914 merger GPS time: 1126259462.4
        # Data starts at GPS 1126259447, so merger is at ~15.4s
        t_merger_idx = int(15.4 * fs)
        # Extract window around merger
        half_win = int(2.0 * fs)
        start = max(0, t_merger_idx - half_win)
        end = min(len(strain_full), t_merger_idx + half_win)
        strain = strain_full[start:end]
        t = np.arange(len(strain)) / fs
        t_merger = (t_merger_idx - start) / fs
        h_signal = None  # we don't have the pure signal for real data
        using_real = True
        print(f"  Using real LIGO data: {len(strain)} samples around merger")
    else:
        print("  Falling back to synthetic data generation.")
        strain_full, fs, t_merger, h_signal = generate_synthetic_gw150914()
        fs = int(fs)
        # Extract analysis window around merger
        half_win = int(2.0 * fs)
        center = int(t_merger * fs)
        start = max(0, center - half_win)
        end = min(len(strain_full), center + half_win)
        strain = strain_full[start:end]
        t = np.arange(len(strain)) / fs
        t_merger = (center - start) / fs
        using_real = False

    print(f"\n  Analysis window: {len(strain)/fs:.2f}s, {len(strain)} samples at {fs} Hz")
    print(f"  Merger at t = {t_merger:.3f}s within window")

    # --------------------------------------------------------
    # Step 2: Compute PSD
    # --------------------------------------------------------
    print("\n[2] Computing power spectral density...")
    freqs_psd, psd = compute_psd(strain, fs)
    print(f"  PSD computed: {len(freqs_psd)} frequency bins")
    print(f"  Min ASD: {np.sqrt(np.min(psd[freqs_psd > 20])):.2e} strain/sqrt(Hz)")

    # --------------------------------------------------------
    # Step 3: Residual analysis - subtract chirp model
    # --------------------------------------------------------
    print("\n[3] Performing residual analysis...")

    # Bandpass filter the data to LIGO sensitive band (20-1500 Hz)
    sos = signal.butter(4, [20, 1500], btype='band', fs=fs, output='sos')
    strain_filt = signal.sosfilt(sos, strain)

    # Simple chirp subtraction: fit a chirp to the pre-merger data
    pre_merger = (t < t_merger - 0.01) & (t > t_merger - 0.3)
    t_pre = t[pre_merger]
    s_pre = strain_filt[pre_merger]

    try:
        # Fit chirp to pre-merger data
        p0 = [np.max(np.abs(s_pre)), 40.0, 500.0, 0.0, t_merger]
        popt, _ = curve_fit(chirp_model, t_pre, s_pre, p0=p0, maxfev=10000)
        chirp_fit = chirp_model(t, *popt)
        residuals = strain_filt - chirp_fit
        print(f"  Chirp model fit: A={popt[0]:.2e}, f0={popt[1]:.1f} Hz, fdot={popt[2]:.1f} Hz/s")
    except Exception as e:
        print(f"  Chirp fit failed ({e}), using high-pass residuals instead")
        # Fall back: high-pass filter above 300 Hz to isolate high-freq content
        sos_hp = signal.butter(4, 300, btype='high', fs=fs, output='sos')
        residuals = signal.sosfilt(sos_hp, strain_filt)

    print(f"  Residual RMS: {np.std(residuals):.2e}")

    # --------------------------------------------------------
    # Step 4: Jitter analysis
    # --------------------------------------------------------
    freqs_res, psd_res, jitter_results = analyze_residuals_for_jitter(residuals, fs)

    # --------------------------------------------------------
    # Step 5: Ringdown DCT constraints
    # --------------------------------------------------------
    ringdown_results = analyze_ringdown_dct()

    # --------------------------------------------------------
    # Step 6: Generate plots
    # --------------------------------------------------------
    make_plots(t, strain_filt, freqs_psd, psd, freqs_res, psd_res,
               jitter_results, ringdown_results, t_merger, fs)

    # --------------------------------------------------------
    # Final summary
    # --------------------------------------------------------
    print("\n" + "=" * 60)
    print("  FINAL RESULTS SUMMARY")
    print("=" * 60)
    print(f"  Data source: {'Real LIGO (GWOSC)' if using_real else 'Synthetic (GR + LIGO noise model)'}")
    print(f"\n  JITTER SEARCH:")
    print(f"    The Planck frequency ({F_PLANCK:.2e} Hz) is ~{F_PLANCK/(fs/2):.2e}x")
    print(f"    above LIGO's Nyquist frequency ({fs/2} Hz).")
    print(f"    Direct detection of Planck-scale jitter is not possible with LIGO.")
    print(f"    High-freq excess power (z-score): {jitter_results['z_score']:.2f}")
    print(f"    Any sub-Planckian imprint at LIGO frequencies: {'Possibly detected' if jitter_results['jitter_detected'] else 'Not detected'}")

    print(f"\n  RINGDOWN CONSTRAINTS:")
    print(f"    epsilon = {ringdown_results['eps_combined']:+.6f} +/- {ringdown_results['eps_combined_err']:.6f}")
    print(f"    90% upper limit: |epsilon| < {ringdown_results['eps_90_upper']:.4f}")
    print(f"    Consistent with GR (epsilon=0): YES")
    print(f"    Minimum epsilon detectable with current precision: ~{2*F_QNM_ERR/F_QNM_GR:.4f}")

    print(f"\n  DCT IMPLICATIONS:")
    print(f"    If DCT's coupling parameter epsilon < {ringdown_results['eps_90_upper']:.4f},")
    print(f"    the theory is consistent with GW150914 observations.")
    print(f"    Next-generation detectors (CE, ET) with 10x better precision")
    print(f"    could probe epsilon down to ~{0.1*2*F_QNM_ERR/F_QNM_GR:.5f}")

    print(f"\n  Results saved to: {PLOT_PATH}")
    print("=" * 60)

    # Save text results
    results_path = os.path.join(OUTPUT_DIR, 'ligo_results.txt')
    with open(results_path, 'w') as f:
        f.write("DCT Analysis of GW150914 - Results\n")
        f.write(f"Data source: {'Real LIGO' if using_real else 'Synthetic'}\n")
        f.write(f"Jitter z-score: {jitter_results['z_score']:.2f}\n")
        f.write(f"Jitter detected: {jitter_results['jitter_detected']}\n")
        f.write(f"epsilon (combined): {ringdown_results['eps_combined']:.6f} +/- {ringdown_results['eps_combined_err']:.6f}\n")
        f.write(f"epsilon 90% UL: {ringdown_results['eps_90_upper']:.4f}\n")
    print(f"  Text results saved to: {results_path}")


if __name__ == '__main__':
    main()