"""What a phonon-mediated superconductor needs for a given Tc.

Uses the Einstein-mode Eliashberg solution Tc = omega * f(lambda, mu*) from eliashberg.py
to give, for each lambda, the mode frequency and the hydrogen Hopfield parameter
eta = lambda * M_H * omega^2 that a target Tc requires. Also tabulates isotherms in the
(lambda, omega) plane for the figures.
"""

import json
import math
from pathlib import Path

import numpy as np

from eliashberg import tc

HERE = Path(__file__).parent
KB = 1.380649e-23
HBAR = 1.054571817e-34
EV = 1.602176634e-19
AMU = 1.66053906660e-27
M_H = 1.00784 * AMU
K_PER_MEV = 11.6045


def m_omega2(omega_kelvin, mass=M_H):
    """M * omega^2 in eV per square angstrom for a mode frequency given in kelvin."""
    w = KB * omega_kelvin / HBAR
    return mass * w * w / (EV / 1e-20)


def f(lam, mu):
    return tc([(lam, 1.0)], mu, 10.0)


if __name__ == "__main__":
    lams = [1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 7.0, 10.0]
    out = {"m_h_omega2_at_1000K": m_omega2(1000.0), "requirements": [], "isotherms": {}}
    for mu in (0.0, 0.10, 0.13, 0.16):
        for lam in lams:
            ratio = f(lam, mu)
            row = {"mu_star": mu, "lambda": lam, "tc_over_omega": ratio}
            for target in (300.0, 400.0):
                w = target / ratio
                row[f"omega_K_for_{int(target)}"] = w
                row[f"omega_meV_for_{int(target)}"] = w / K_PER_MEV
                row[f"eta_H_for_{int(target)}"] = lam * m_omega2(w)
            out["requirements"].append(row)

    # Large-lambda limit of the Hopfield requirement at mu* = 0: Tc = 0.1827 sqrt(eta / M).
    out["eta_H_asymptote_300"] = m_omega2(300.0 / 0.1827)
    out["eta_H_asymptote_400"] = m_omega2(400.0 / 0.1827)

    grid = np.concatenate([np.arange(0.5, 2.0, 0.05), np.arange(2.0, 4.51, 0.1)])
    ratios = [f(float(l), 0.13) for l in grid]
    for target in (39.0, 100.0, 200.0, 300.0):
        out["isotherms"][str(int(target))] = [
            {"lambda": float(l), "omega_K": target / r} for l, r in zip(grid, ratios) if r > 0
        ]
    out["isotherm_mu_star"] = 0.13

    (HERE / "requirements.json").write_text(json.dumps(out, indent=1))

    print(f"M_H omega^2 at 1000 K = {out['m_h_omega2_at_1000K']:.4f} eV/A^2")
    print(f"eta_H asymptote for 300 K (mu*=0): {out['eta_H_asymptote_300']:.2f} eV/A^2; for 400 K: {out['eta_H_asymptote_400']:.2f}")
    print(f"{'mu*':>5} {'lambda':>6} {'Tc/w':>7} {'w300 K':>8} {'w300 meV':>9} {'eta300':>7} {'eta400':>7}")
    for r in out["requirements"]:
        print(f"{r['mu_star']:5.2f} {r['lambda']:6.1f} {r['tc_over_omega']:7.4f} {r['omega_K_for_300']:8.0f} {r['omega_meV_for_300']:9.1f} {r['eta_H_for_300']:7.2f} {r['eta_H_for_400']:7.2f}")