# Thanks to ChatGPT and the Fortran-90 code from arxiv:1303.5961,
#     https://code.google.com/archive/p/cosmonom/downloads
# here's how to get cosmological look-back time from redshift in Python:

import matplotlib.pyplot as plt
from numpy import arcsinh

Om = 0.3153  # from Planck Collaboration 2018: https://arxiv.org/abs/1807.06209
furthest_z = 14.44  # 2025 record: https://arxiv.org/abs/2505.11263

def make_curve(H0, color_main, color_dark, label):
    # Radiation density parameter from Planck Collaboration 2018 using CMB
    # temperature from Fixsen 2009, N_eff = 3.046, and Omega_r * h^2 ≈ 4.15e-5:
    OL = 1.0 - Om - 0.415 / (H0**2)

    def age_at_z(z):
        return (2/3) * arcsinh(((OL / (Om * (1 + z)**3))**0.5)) / (H0 * OL**0.5) * 977.8
    age0 = age_at_z(0)
    def zt(z): return age0 - age_at_z(z)

    rs = [z * 20 / 299 for z in range(300)]
    lb = [zt(z) for z in rs]
    # Split into pre- and post-MoM-z14 ranges
    x1, y1 = zip(*[(x, y) for x, y in zip(rs, lb) if x <= furthest_z])
    x2, y2 = zip(*[(x, y) for x, y in zip(rs, lb) if x > furthest_z])
    plt.plot(x1, y1, color=color_main, label=f"{label}: $H_0$={H0} km/s/Mpc")
    plt.plot(x2, y2, color=color_dark)

plt.figure(figsize=(9,6))
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

make_curve(67.36, 'cyan', 'teal', 'Planck Collaboration')
make_curve(69.32, 'red', 'darkred', 'Consensus compromise')
make_curve(70, 'orange', 'darkorange', 'Single digit precision')
make_curve(73.2, 'blue', 'midnightblue', 'SH0ES Team')

plt.title('Look-back Time by Redshift for Different Hubble Constants')
plt.xlabel('Redshift z: (observed λ - expected λ) / expected λ')
plt.ylabel('Look-back Time in billion years (Gyr)')
legend = plt.legend(title="Assuming $Ω_m$=0.3153, as per arXiv:1807.06209")
plt.grid(True, color='lightgray')

plt.ylim(0, 13.9)
plt.xlim(0, 20)
plt.xticks(range(21))
plt.yticks(range(14))

plt.text(0.5, 13.78, "— Big Bang: 13.78 Gyr (as per Planck Collaboration) —", va='center')

eqs = (r"$\mathrm{ageAtRedshift}(z)=\int_{z}^{\infty}\frac{dz'}{(1+z')\sqrt{\Omega_{\Lambda}+\Omega_{m}(1+z')^{3}}}\,\frac{977.8}{H_{0}}$" "\n"
    r"$=\frac{2\,\sinh^{-1}\!\left(\dfrac{\sqrt{\Omega_{\Lambda}/\Omega_{m}}}{(1+z)^{3/2}}\right)\,977.8}{3H_{0}\sqrt{\Omega_{\Lambda}}}$ Gyr, as per arXiv:gr-qc/0508073." "\n\n"
    r"$\mathrm{lookBackTime}(z)=\mathrm{ageAtRedshift}(0)-\mathrm{ageAtRedshift}(z)$." "\n\n"
    r"$\Omega_\Lambda = 1.0 - \Omega_m - \frac{0.415}{H_0^2}$, per Planck 2018 and Fixsen 2009")
fig = plt.gcf()
fig.text(0.5, 0.5, eqs, ha='center', va='center', fontsize=10,
    bbox=dict(boxstyle='round,pad=0.4', facecolor='white', alpha=0.9, linewidth=0))

plt.text(furthest_z, 10.5, 'Furthest observation as of 2025:\n'
    'the galaxy MoM-z14, at z=14.44,\nor about 13.5 billion years ago',
    ha='center')

outpath_multi = 'lookback-time-by-redshift-H0-comparison.png'
plt.savefig(outpath_multi, bbox_inches='tight')
plt.show()  # https://i.ibb.co/GfLjrWN3/lookback-time-by-redshift-H0-comparison-2.png