检查真实 GW150914 探测器应变

from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from scipy.signal import welch, butter, sosfiltfilt
from IPython.display import display as _display
def display(value):
    if isinstance(value,pd.DataFrame):
        table=value.to_html(max_rows=12,escape=True)
        style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
        return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
    return _display(value)
RECIPE_ID='ph05-gw150914'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
                      "axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")

Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6
Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads
Core data mode: embedded snapshot; no network requests

def load_snapshot(recipe_id, expected_sha256, snapshot_key):
    """Read embedded data from this notebook, including a renamed copy."""
    for path in sorted(Path.cwd().glob("*.ipynb")):
        try:
            notebook = json.loads(path.read_text(encoding="utf-8"))
            metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
            blob = metadata.get("snapshots", {}).get(snapshot_key)
            if metadata.get("recipe_id") != recipe_id or not blob:
                continue
            if blob.get("sha256") != expected_sha256:
                continue
            raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
            if hashlib.sha256(raw).hexdigest() != expected_sha256:
                raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
            return json.loads(raw)
        except (OSError, json.JSONDecodeError):
            continue
    raise FileNotFoundError(
        "Save/open the complete .ipynb in the active working folder. "
        "The embedded snapshot could not be located for " + recipe_id)

SNAPSHOT_HASHES={'gw150914': '4068109f074163cc9a5b99edc94606ede7862f9b5907e7903c307fe7fd713edf'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

Embedded snapshots: ['gw150914']

gw=snapshot("gw150914");strain=np.asarray(gw["strain"],float);fs=gw["sample_rate"]
assert len(strain)==fs*gw["duration"] and np.isfinite(strain).all()
relative_time=gw["gps_start"]+np.arange(len(strain))/fs-gw["event_gps"]
centered=strain-strain.mean()
print("Release:",gw["version"],"; detector:",gw["detector"],"; samples:",len(strain),"; rate:",fs,"Hz")
fig,ax=plt.subplots(figsize=(9,4));ax.plot(relative_time,centered,lw=.5)
ax.axvline(0,color="red",linestyle="--",label="Catalogued event time")
ax.set(xlabel="Time relative to catalogued event (s)",ylabel="H1 strain (dimensionless)",title="Frozen detector measurements");ax.legend();plt.tight_layout();plt.show()

Release: GW150914 v3 / 4KHZ_R1 ; detector: H1 ; samples: 131072 ; rate: 4096 Hz

frequency,psd=welch(centered,fs=fs,nperseg=fs*4,noverlap=fs*2,scaling="density")
positive=frequency>0
fig,ax=plt.subplots(figsize=(8,5));ax.loglog(frequency[positive],np.sqrt(psd[positive]))
ax.axvspan(35,350,color="orange",alpha=.12,label="Declared display band")
ax.set(xlim=(10,1500),xlabel="Frequency (Hz)",ylabel="ASD (strain / sqrt(Hz))",title="Welch spectrum of the frozen record");ax.legend();plt.tight_layout();plt.show()

sos=butter(4,[35,350],btype="bandpass",fs=fs,output="sos")
filtered=sosfiltfilt(sos,centered)
window=(relative_time>=-.25)&(relative_time<=.15)
fig,ax=plt.subplots(figsize=(9,4));ax.plot(relative_time[window],1e21*filtered[window],lw=1)
ax.axvline(0,color="gray",linestyle="--")
ax.set(xlabel="Time relative to catalogued event (s)",ylabel="Filtered strain (10⁻²¹)",title="Declared band-pass; interpretation requires the published analysis");plt.tight_layout();plt.show()

transform=np.fft.rfft(centered);squared=np.abs(transform)**2
fft_variance=(squared[0]+squared[-1]+2*squared[1:-1].sum())/len(centered)**2
assert np.isclose(fft_variance,np.mean(centered**2),rtol=1e-12,atol=0)
assert np.isfinite(filtered).all() and np.all(psd>=0)
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"relative_time_s":relative_time,"centered_strain":centered,"filtered_strain":filtered}).to_csv(OUTPUT_DIR/"strain-display.csv",index=False)
pd.DataFrame({"frequency_Hz":frequency,"PSD_strain_squared_per_Hz":psd}).to_csv(OUTPUT_DIR/"welch-spectrum.csv",index=False)
(OUTPUT_DIR/"processing.json").write_text(json.dumps({"detector":gw["detector"],"release":gw["version"],"event_GPS":gw["event_gps"],"sample_rate_Hz":fs,"band_Hz":[35,350],"filter":"4th-order Butterworth SOS, forward/backward","inference":"display processing; no independent detection significance"},indent=2))
print("Parseval check passed; source credit: GWOSC, LIGO, Virgo and KAGRA. Exports:",OUTPUT_DIR)

Parseval check passed; source credit: GWOSC, LIGO, Virgo and KAGRA. Exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ph05-gw150914

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://gwosc.org/eventapi/json/GWTC-1-confident/GW150914/v3/',timeout=20)
        response.raise_for_status()
        print("Source:",response.url)
        print(response.text[:700])
    except requests.RequestException as error:
        print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
    print("Online extension skipped; embedded core data unchanged")

Online extension skipped; embedded core data unchanged