"""Eliashberg Tc for published hydride spectral functions, coarse-grained into Einstein modes.

The bins (coupling, frequency) were read from the cumulative lambda curves printed in the
cited papers (see research/gap-ambient-hopfield-eta.md, facts H7, H19 to H21). For each
spectrum this reports lambda, omega_log, the second moment, the Hopfield sum
S = lambda <omega^2> in hydrogen units (eV/A^2), the asymptote 0.1827 sqrt(S), the
Eliashberg Tc, and the efficiency Phi = Tc / asymptote.
"""

import json
import math
from pathlib import Path

from eliashberg import moments, tc

HERE = Path(__file__).parent
K_PER_MEV = 11.6045
M_H_OMEGA2_1000K = 1.7904
CM = 0.1239842  # meV per cm^-1

SPECTRA = {
    "Mg2IrH6 (Sanna et al.)": {"bins": [(0.75, 19), (0.32, 55), (0.76, 85), (0.17, 161), (0.07, 175), (0.06, 192)], "published": "77 K (first-principles Coulomb), lambda 2.1"},
    "Mg2IrH6 (Dolui et al.)": {"bins": [(0.76, 22), (0.45, 50), (0.90, 88), (0.20, 170), (0.21, 193)], "published": "175, 170, 160 K at mu* = 0.10, 0.125, 0.16"},
    "Mg2RhH6 (Sanna et al.)": {"bins": [(0.30, 24), (0.27, 52), (0.48, 88), (0.09, 168), (0.04, 175), (0.07, 183), (0.02, 190)], "published": "48.5 K (first-principles Coulomb), lambda 1.3"},
    "Mg2PtH6 (Sanna et al.)": {"bins": [(0.21, 22), (0.97, 72), (0.11, 119), (0.08, 138), (0.01, 158)], "published": "80.4 K (first-principles Coulomb), lambda 1.4"},
    "PdH, anharmonic (Errea et al.)": {"bins": [(0.10, 130 * CM), (0.22, 555 * CM), (0.05, 707 * CM), (0.02, 775 * CM), (0.01, 870 * CM)], "published": "5.0 K at mu* = 0.085, lambda 0.40"},
}


def phi_einstein(lam, mu):
    return tc([(lam, 1.0)], mu, 10.0) / (0.1827 * math.sqrt(lam))


if __name__ == "__main__":
    out = {"spectra": [], "phi_einstein": []}
    for name, s in SPECTRA.items():
        modes = [(l, w) for l, w in s["bins"]]
        lam, w_log, w2 = moments(modes)
        S = lam * (w2 * K_PER_MEV / 1000.0) ** 2 * M_H_OMEGA2_1000K
        stiff = sum(l * w * w for l, w in modes if w >= 150) / sum(l * w * w for l, w in modes)
        lam_stiff = sum(l for l, w in modes if w >= 150)
        row = {"name": name, "lambda": lam, "omega_log_meV": w_log, "omega_2_meV": w2, "shape": w_log / w2,
               "S_eV_A2": S, "asymptote_K": 182.7 * math.sqrt(S / M_H_OMEGA2_1000K),
               "lambda_above_150meV": lam_stiff, "share_of_S_above_150meV": stiff, "published": s["published"], "tc": {}}
        for mu in (0.10, 0.13, 0.16):
            t = tc(modes, mu, 10.0) * K_PER_MEV
            row["tc"][f"{mu:.2f}"] = t
        row["phi_0.13"] = row["tc"]["0.13"] / row["asymptote_K"]
        row["phi_einstein_same_lambda_0.13"] = phi_einstein(lam, 0.13)
        out["spectra"].append(row)
    for mu in (0.10, 0.13, 0.16):
        out["phi_einstein"].append({"mu_star": mu, **{str(l): phi_einstein(l, mu) for l in (1, 1.5, 2, 3, 4, 5, 10)}})
    # sensitivity of Tc to lambda at fixed S (Einstein): d ln Tc / d ln lambda
    out["dlnTc_dlnlambda_fixed_S"] = {}
    for l in (1.5, 2.0, 3.0, 4.0):
        f = lambda x: tc([(x, 1.0 / math.sqrt(x))], 0.13, 10.0)
        out["dlnTc_dlnlambda_fixed_S"][str(l)] = (math.log(f(l * 1.02)) - math.log(f(l / 1.02))) / (2 * math.log(1.02))
    (HERE / "spectra.json").write_text(json.dumps(out, indent=1))
    for r in out["spectra"]:
        print(f"{r['name']:32s} lam={r['lambda']:.2f} wlog={r['omega_log_meV']:.0f} w2={r['omega_2_meV']:.0f} meV shape={r['shape']:.2f} S={r['S_eV_A2']:.2f} asym={r['asymptote_K']:.0f} K  Tc={r['tc']['0.10']:.0f}/{r['tc']['0.13']:.0f}/{r['tc']['0.16']:.0f} K  Phi={r['phi_0.13']:.2f} (Einstein {r['phi_einstein_same_lambda_0.13']:.2f})  lam>150meV={r['lambda_above_150meV']:.2f} share of S {r['share_of_S_above_150meV']:.2f} | {r['published']}")
    for p in out["phi_einstein"]:
        print("Phi Einstein mu*=", p["mu_star"], {k: round(v, 2) for k, v in p.items() if k != "mu_star"})
    print("d ln Tc / d ln lambda at fixed S:", {k: round(v, 2) for k, v in out["dlnTc_dlnlambda_fixed_S"].items()})