"""
dct_nanograv_spectrum.py

DCT NANOGrav 15-yr stochastic-GW spectrum test.
Measured spectral index gamma_PTA = 3.2 +/- 0.6.
The supermassive black-hole binary (SMBHB) inspiral prediction is gamma = 13/3.
DCT inherits the SMBHB inspiral exactly, plus a scalar-mode contribution at
f ~ 7.6e-20 Hz which lies BELOW NANOGrav's 1/yr (~ 3.17e-8 Hz) sensitivity.
Both DCT and SMBHB sit ~ 1.9 sigma from the central value: status = NEUTRAL.

Reproduces: gamma_DCT = gamma_SMBHB = 13/3 ~ 4.33; sigma_PTA ~ 1.89.
Public archive: NANOGrav 15-yr posterior on gamma_PTA (Agazie et al. 2023).
Run with: python3 dct_nanograv_spectrum.py
"""

import numpy as np


# -------------------------------------------------------------------
# Canonical DCT constants (verbatim).
# -------------------------------------------------------------------
P_0       = 0.851
H0_PLANCK = 67.36         # km/s/Mpc
GAMMA_PTA_OBS   = 3.2
GAMMA_PTA_SIGMA = 0.6
GAMMA_SMBHB     = 13.0 / 3.0   # standard inspiral spectral index

# DCT scalar-mode characteristic frequency (audit-locked).
F_SCALAR_HZ     = 7.6e-20
F_NANOGRAV_LO   = 1.0 / (365.25 * 24.0 * 3600.0)   # 1/yr ~ 3.17e-8 Hz


def main() -> None:
    sigma_smbhb = abs(GAMMA_SMBHB - GAMMA_PTA_OBS) / GAMMA_PTA_SIGMA
    # DCT inherits SMBHB exactly (the scalar mode is below sensitivity).
    gamma_dct = GAMMA_SMBHB
    sigma_dct = abs(gamma_dct - GAMMA_PTA_OBS) / GAMMA_PTA_SIGMA

    print("=" * 60)
    print("DCT NANOGrav 15-yr spectrum test")
    print("=" * 60)
    print(f"Measured     gamma_PTA = {GAMMA_PTA_OBS} +/- {GAMMA_PTA_SIGMA}")
    print(f"SMBHB pred.  gamma     = 13/3 = {GAMMA_SMBHB:.4f}")
    print(f"DCT  pred.   gamma     = {gamma_dct:.4f}  (inherits SMBHB)")
    print()
    print(f"Tension vs SMBHB : {sigma_smbhb:.2f} sigma")
    print(f"Tension vs DCT   : {sigma_dct:.2f} sigma")
    print()
    print("DCT scalar-mode characteristic frequency:")
    print(f"  f_scalar  = {F_SCALAR_HZ:.2e} Hz")
    print(f"  f_NG_low  = {F_NANOGRAV_LO:.2e} Hz (1/yr cadence floor)")
    print(f"  ratio     = {F_SCALAR_HZ / F_NANOGRAV_LO:.2e}  (well below sensitivity)")
    print()
    print("Status: NEUTRAL.  DCT and SMBHB both at ~1.9 sigma; data cannot")
    print("        distinguish them, and the scalar mode is undetectable to PTAs.")
    print("=" * 60)


if __name__ == "__main__":
    main()