"""Closure comparison of the model section: hydrogen stretch share of S in four published spectra and the change in Tc when S is redistributed isotropically.

Run with the project environment: .venv/bin/python numerics/derived/closure.py
The printed output is stored beside this file as closure.out.txt.
"""
import pathlib
import sys, math
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
from eliashberg import tc, moments
from spectra import SPECTRA
K=11.6045
for name in ["Mg2IrH6 (Sanna et al.)","Mg2IrH6 (Dolui et al.)","Mg2RhH6 (Sanna et al.)","Mg2PtH6 (Sanna et al.)"]:
    bins = SPECTRA[name]["bins"]
    # hydrogen bins: above 40 meV; stretches: >=150 meV
    H = [(l,w) for l,w in bins if w>=40]
    M = [(l,w) for l,w in bins if w<40]
    SH = sum(l*w*w for l,w in H)
    Sst = sum(l*w*w for l,w in H if w>=150)
    share = Sst/SH
    t0 = tc(bins,0.13,10.0)*K
    lam0 = sum(l for l,_ in bins)
    # isotropic: stretches get 1/3 of hydrogen S, rest pro rata
    new = list(M)
    for l,w in H:
        if w>=150: new.append((l*(1/3)/share, w))
        else: new.append((l*(2/3)/(1-share), w))
    lam1 = sum(l for l,_ in new)
    S0 = sum(l*w*w for l,w in bins); S1=sum(l*w*w for l,w in new)
    t1 = tc(new,0.13,10.0)*K
    print(f"{name:28s} stretch share of H S {share*100:.1f}%  lam {lam0:.2f}->{lam1:.2f}  Tc {t0:.1f}->{t1:.1f} K ({(t1/t0-1)*100:+.1f}%)  S ratio {S1/S0:.4f}")