# %% [markdown]
# # FRCFD Deep Dive: Does the 200.2 Hz Anchor "Breathe" with the Orbital Motion?
#
# - Installs `gwpy` if not already present
# - Extracts amplitude envelope of the 200.2 Hz resonance from LIGO data
# - Compares envelope modulation with expected orbital frequency (chirp mass)
# - Cross‑correlates H1 and L1 envelopes to identify astrophysical signal
# - Checks Phase‑Amplitude Coupling (PAC) as evidence of "Time Dilation Friction"
# ---------------------------
# 0. Install gwpy (if needed)
# ---------------------------
import subprocess
import sys
try:
import gwpy
print("gwpy already installed.")
except ImportError:
print("Installing gwpy...")
subprocess.check_call([sys.executable, "-m", "pip", "install", "gwpy"])
import gwpy
print("gwpy installed successfully.")
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt, hilbert, correlate
from gwpy.timeseries import TimeSeries
# ---------------------------
# 1. Parameters (Refined)
# ---------------------------
# GW150914 chirp mass final estimate (Abbott et al. 2016)
M_CHIRP = 28.1 # solar masses
G = 6.67430e-11
c = 299792458
M_sun = 1.989e30
M_chirp_kg = M_CHIRP * M_sun
F_ANCHOR = 200.2 # Hz (FRCFD invariant)
BANDWIDTH = 10 # Hz (ultra‑narrow: 195–205 Hz)
DURATION = 1.0 # seconds of data before merger (ends 50 ms before peak)
T0 = 1126259462.423 # merger peak GPS time
START_TIME = T0 - DURATION
END_TIME = T0 - 0.05 # avoid strongest merger non‑stationarity
# Light travel time between LIGO Hanford (H1) and Livingston (L1)
DELAY_H1_L1 = 0.0069 # seconds (positive = H1 signal arrives later)
# ---------------------------
# 2. Helper: Orbital frequency from chirp mass
# ---------------------------
def orbital_frequency(t_relative_to_merger):
"""
t_relative_to_merger: negative seconds (e.g., -0.5 = 0.5 s before merger)
Returns f_orb in Hz.
"""
tau = -t_relative_to_merger
tau = np.maximum(tau, 1e-6)
GMc3 = (G * M_chirp_kg) / c**3
factor = (5 / 256) * (1 / tau)
f_orb = (1 / np.pi) * (factor)**(3/8) * (GMc3)**(-5/8)
return f_orb
# ---------------------------
# 3. Function to extract amplitude envelope for a given detector
# ---------------------------
def extract_envelope(detector, start_gps, end_gps, f_anchor, bandwidth, fs_target=None):
"""
Fetches LIGO data, bandpass filters around f_anchor, and returns:
time_axis (seconds relative to T0), amplitude envelope, strain data
"""
data = TimeSeries.fetch_open_data(detector, start_gps, end_gps, verbose=True)
fs = data.sample_rate.value
if fs_target is not None and fs != fs_target:
data = data.resample(fs_target)
fs = fs_target
# Bandpass around anchor
nyquist = fs / 2
low = max(f_anchor - bandwidth/2, 0.1)
high = f_anchor + bandwidth/2
b, a = butter(4, [low/nyquist, high/nyquist], btype='band')
filtered = filtfilt(b, a, data.value)
# Heterodyne demodulation (complex downconversion)
t = np.arange(len(filtered)) / fs
analytic = filtered * np.exp(-1j * 2 * np.pi * f_anchor * t)
# Low‑pass to keep modulation (cutoff = 200 Hz, covers orbital sweep)
b_lp, a_lp = butter(4, 200/nyquist, btype='low')
envelope_complex = filtfilt(b_lp, a_lp, analytic)
amplitude_envelope = np.abs(envelope_complex)
# Time axis relative to merger (T0)
time_rel = t + start_gps - T0 # negative before merger
return time_rel, amplitude_envelope, filtered, fs
# ---------------------------
# 4. Extract envelopes for H1 and L1 (time‑aligned to same GPS)
# ---------------------------
print("Fetching H1 data...")
t_h1, env_h1, strain_h1, fs = extract_envelope('H1', START_TIME, END_TIME, F_ANCHOR, BANDWIDTH)
print("Fetching L1 data...")
t_l1, env_l1, strain_l1, _ = extract_envelope('L1', START_TIME, END_TIME, F_ANCHOR, BANDWIDTH)
# Check that time axes match (they should, because same GPS window)
if not np.allclose(t_h1, t_l1):
print("Warning: time axes differ – resampling may be needed.")
# Simple resample to common time grid (using linear interpolation)
common_t = np.linspace(max(t_h1[0], t_l1[0]), min(t_h1[-1], t_l1[-1]), min(len(t_h1), len(t_l1)))
env_h1 = np.interp(common_t, t_h1, env_h1)
env_l1 = np.interp(common_t, t_l1, env_l1)
t = common_t
else:
t = t_h1
# Remove DC offset to see fluctuations
env_h1_ac = env_h1 - np.mean(env_h1)
env_l1_ac = env_l1 - np.mean(env_l1)
# ---------------------------
# 5. Cross‑correlation of H1 and L1 envelopes
# ---------------------------
corr = correlate(env_h1_ac, env_l1_ac, mode='same')
lags = np.linspace(-len(env_h1_ac)//2, len(env_h1_ac)//2, len(corr)) / fs
plt.figure(figsize=(12, 5))
plt.plot(lags * 1000, corr) # lags in ms
plt.axvline(DELAY_H1_L1 * 1000, color='r', linestyle='--', label=f'Expected delay H1→L1 = {DELAY_H1_L1*1000:.1f} ms')
plt.axvline(-DELAY_H1_L1 * 1000, color='r', linestyle='--')
plt.xlabel('Time lag (ms)')
plt.ylabel('Cross‑correlation')
plt.title('H1 vs L1 – 200.2 Hz Envelope Cross‑correlation')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()
# Find actual peak lag
peak_idx = np.argmax(np.abs(corr))
peak_lag_ms = lags[peak_idx] * 1000
print(f"Peak cross‑correlation at lag = {peak_lag_ms:.2f} ms")
print(f"Expected |lag| = {DELAY_H1_L1*1000:.1f} ms")
if np.abs(peak_lag_ms - DELAY_H1_L1*1000) < 5:
print("✅ Envelope modulation is consistent with astrophysical light travel time.")
else:
print("⚠️ Peak lag does not match expected delay – modulation may be local noise.")
# ---------------------------
# 6. Plot envelopes and orbital frequency
# ---------------------------
f_orb_expected = orbital_frequency(t)
fig, ax1 = plt.subplots(figsize=(14, 6))
ax1.plot(t, env_h1, label='H1 envelope', color='blue', alpha=0.7)
ax1.plot(t, env_l1, label='L1 envelope', color='orange', alpha=0.7)
ax1.set_xlabel('Time relative to merger (s)')
ax1.set_ylabel('200.2 Hz Amplitude (strain)')
ax1.legend(loc='upper left')
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-DURATION, 0)
ax2 = ax1.twinx()
ax2.plot(t, f_orb_expected, 'r--', label='Expected f_orb', linewidth=1.5)
ax2.set_ylabel('Orbital frequency (Hz)', color='red')
ax2.tick_params(axis='y', labelcolor='red')
ax2.legend(loc='upper right')
plt.title(f'200.2 Hz Amplitude Envelope vs. Orbital Frequency (GW150914)')
plt.tight_layout()
plt.show()
# ---------------------------
# 7. Phase‑Amplitude Coupling (PAC) check
# ---------------------------
# Use the main GW signal (broadband 30–400 Hz) to get the instantaneous phase
# of the orbital motion. Then bin the 200.2 Hz envelope by that phase.
data_gw = TimeSeries.fetch_open_data('H1', START_TIME, END_TIME, verbose=False)
b_broad, a_broad = butter(4, [30/fs, 400/fs], btype='band')
gw_broad = filtfilt(b_broad, a_broad, data_gw.value)
phase_gw = np.angle(hilbert(gw_broad)) # instantaneous phase of the GW
# Resample envelope to same time grid (if needed)
if len(t) != len(phase_gw):
phase_gw_resampled = np.interp(t, np.arange(len(phase_gw))/fs + START_TIME - T0, phase_gw)
else:
phase_gw_resampled = phase_gw
# Bin envelope amplitude by phase (0 to 2π)
n_bins = 36
bins = np.linspace(-np.pi, np.pi, n_bins+1)
bin_centers = (bins[:-1] + bins[1:]) / 2
mean_amp = np.zeros(n_bins)
for i in range(n_bins):
mask = (phase_gw_resampled >= bins[i]) & (phase_gw_resampled < bins[i+1])
if np.any(mask):
mean_amp[i] = np.mean(env_h1[mask])
plt.figure(figsize=(8, 5))
plt.polar(bin_centers, mean_amp / np.max(mean_amp))
plt.title('Phase‑Amplitude Coupling: 200.2 Hz Envelope vs. GW Phase')
plt.grid(True)
plt.show()
# Simple modulation depth
mod_depth = (np.max(mean_amp) - np.min(mean_amp)) / np.mean(mean_amp)
print(f"Modulation depth vs. GW phase: {mod_depth:.2f}")
if mod_depth > 0.2:
print("✅ Significant phase‑amplitude coupling detected – supports FRCFD regulator.")
else:
print("⚠️ Weak phase‑amplitude coupling – regulator may be subtle or absent.")
# ---------------------------
# 8. Optional: FFT of envelope to look for orbital frequency peak
# ---------------------------
# Use the H1 envelope (common features are astrophysical)
envelope_segment = env_h1[t < -0.05] # avoid near‑merger edge
if len(envelope_segment) > 100:
fft_env = np.fft.rfft(envelope_segment - np.mean(envelope_segment))
freqs_env = np.fft.rfftfreq(len(envelope_segment), d=1/fs)
plt.figure(figsize=(10, 4))
plt.semilogy(freqs_env, np.abs(fft_env))
plt.xlim(0, 200)
plt.xlabel('Modulation frequency (Hz)')
plt.ylabel('Power')
plt.title('Power spectrum of 200.2 Hz amplitude envelope (H1)')
plt.grid(True, alpha=0.3)
# Overlay expected orbital frequency range at the time
f_min = orbital_frequency(-0.9)
f_max = orbital_frequency(-0.05)
plt.axvspan(f_min, f_max, alpha=0.3, color='red', label=f'Expected f_orb range: {f_min:.0f}–{f_max:.0f} Hz')
plt.legend()
plt.show()
peak_mod_freq = freqs_env[np.argmax(np.abs(fft_env)[1:])+1]
print(f"Dominant envelope modulation frequency: {peak_mod_freq:.1f} Hz")
if f_min <= peak_mod_freq <= f_max:
print("✅ Modulation frequency falls within expected orbital sweep.")
else:
print("⚠️ Modulation frequency outside expected orbital range.")
else:
print("Not enough envelope data for FFT.")
