"""Single-mode Eliashberg table at mu*=0.13 (cutoff 10 omega_E), the floor on eta_H, lambda at 317 meV, and the functional derivative per unit of the Hopfield sum. Used in the section on what 300 K costs.

Run with the project environment: .venv/bin/python numerics/derived/pairing_single_mode.py
The printed output is stored beside this file as pairing_single_mode.out.txt.
"""
import pathlib
import sys, math, json
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
import numpy as np
from eliashberg import tc, functional_derivative, allen_dynes, moments
from scipy.optimize import brentq
K=11.6045
MH=1.7904  # eV/A^2 at 1000 K (script constant)
# (a) Einstein table at mu*=0.13
for lam in (1,1.5,2,2.5,3,4,5,10):
    f=tc([(lam,1.0)],0.13,10.0)
    wE=300/f/K
    eta=lam*(wE*K/1000)**2*MH
    print(f"lam={lam} f={f:.4f} wE={wE:.1f} meV {wE*K:.0f} K eta={eta:.2f} Phi={f/(0.1827*math.sqrt(lam)):.3f}")
print('floor', (300/136.54)**2, (300/ (182.7/math.sqrt(MH)))**2, 'prefactor', 182.7/math.sqrt(MH))
# lambda for 317 meV
g=lambda lam: tc([(lam,1.0)],0.13,10.0)-300/(317*K)
print('lambda at 317 meV', brentq(g,0.8,1.5))
# (b) functional derivative fine grid
t0, fd = functional_derivative([(2.0,1.0)],0.13,np.arange(0.04,3.0,0.01))
pk=max(fd,key=lambda d:d['dtc_dalpha2f'])
print('peak', pk['omega_over_tc'], 'tc',t0, 'meV for 300K', pk['omega_over_tc']*300/K)
mx=pk['dtc_dalpha2f']
within=[d['omega_over_tc'] for d in fd if d['dtc_dalpha2f']>=0.99*mx]
print('within 1%:', min(within), max(within))
# per unit S: derivative/(2 Omega)
def perS(x):
    d=min(fd,key=lambda d:abs(d['omega_over_tc']-x))
    return d['dtc_dalpha2f']/(2*d['omega_over_tc']*t0), d['omega_over_tc']
ref=perS(0.25)
for x in (0.25,0.5,1,2,4,6,7,10,13):
    v,xx=perS(x)
    print(f"per-S at {xx:.2f} kTc: {v/ref[0]:.3f}")