PUH-BrianMartell puh_kelly_dcbh_eccentric_gw_sim.py- Updated New Py Code

puh_kelly_dcbh_eccentric_gw_sim.py

import numpy as np
import matplotlib.pyplot as plt

# PUH v11: Kelly DCBH Eccentric GW Sim — PSR Merger Inspiral
t = np.linspace(0, 10, 1000)  # Time to merger (arb.)
e0 = 0.8  # Initial eccentricity (PSR asymmetry)
a0 = 100  # Semi-major (knot cluster)

# Eccentric decay de/dt ~ - (304/15) e^5 / (1-e^2)^{5/2} (toy Peters)
de_dt = - (304/15) * e0**5 / (1 - e0**2)**(5/2)
e = e0 + de_dt * t  # Simplified decay

# h ~ (G M / c^2 r) ~ 1/a (arb.)
h = 1 / a0 * (1 - e**2)**(-3/2)  # Quadrupole approx

plt.figure(figsize=(10,6))
plt.subplot(1,2,1)
plt.plot(t, e, label='Eccentricity Decay', color='cyan', lw=2)
plt.axhline(0, color='red', ls=':', label='Merger e=0')
plt.xlabel('Time t (arb.)'); plt.ylabel('Eccentricity e')
plt.title('PUH v11: DCBH Eccentric Decay')

plt.subplot(1,2,2)
plt.plot(t, h, label='Strain h(t)', color='gold', lw=2)
plt.xlabel('Time t (arb.)'); plt.ylabel('Strain h (arb.)')
plt.title('GW Burst Signature')
plt.tight_layout()
plt.savefig('puh_kelly_dcbh_eccentric_gw_sim.png', dpi=300)
plt.show()

print(f"Peak h ~{np.max(h):.2f} at e~{e[np.argmax(h)]:.2f} — LISA detectable.")