"""Frequency-cap rows, the Leavens check, Phi at large lambda, and the response of Tc to scaled force constants and scaled eta for the published spectra. Used in the section on what 300 K costs.

Run with the project environment: .venv/bin/python numerics/derived/pairing_spectra_response.py
The printed output is stored beside this file as pairing_spectra_response.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, allen_dynes, moments
from spectra import SPECTRA
from scipy.optimize import brentq, minimize_scalar
K=11.6045; MH=1.7904
# Esterlis cap: Tc/omega_E = c
for c in (0.08,0.10,0.12):
    lam=brentq(lambda l: tc([(l,1.0)],0.13,10.0)-c,0.5,3.0)
    wE=300/c  # K
    eta=lam*(wE/1000)**2*MH
    phi=c/(0.1827*math.sqrt(lam))
    # derivative at fixed S
    f=lambda x: tc([(x,1.0/math.sqrt(x))],0.13,10.0)
    d=(math.log(f(lam*1.02))-math.log(f(lam/1.02)))/(2*math.log(1.02))
    print(f"cap {c}: lam={lam:.3f} wE={wE:.0f} K ({wE/K:.0f} meV) eta={eta:.1f} Phi={phi:.3f} dlnTc/dlnlam={d:.2f}")
# Leavens: max over lambda of Tc/A, A = lam*omega/2 for Einstein
for mu in (0.0,0.13):
    r=minimize_scalar(lambda l: -tc([(l,1.0)],mu,10.0)/(l/2), bounds=(0.5,3.0), method='bounded')
    print('Leavens mu',mu,'lam*',r.x,'Tc/A max',-r.fun, 'A for 300K meV', 300/K/(-r.fun))
# Phi at lambda 3.7 and 10
for l in (3.7,4.7,10):
    print('Phi',l, tc([(l,1.0)],0.13,10.0)/(0.1827*math.sqrt(l)))
# real spectra: scale force constants by s at fixed S: omega -> omega*sqrt(s), lam -> lam/s
print()
for name,s in SPECTRA.items():
    modes=s['bins']
    lam,wl,w2=moments(modes)
    t0=tc(modes,0.13,10.0)
    S=sum(l*w*w for l,w in modes)
    out=[]
    for sc in (1.3,1/1.3):
        m2=[(l/sc, w*math.sqrt(sc)) for l,w in modes]
        t=tc(m2,0.13,10.0)
        out.append((t/t0-1)*100)
    # eta scaled at fixed frequency
    oute=[]
    for sc in (1.3,1/1.3):
        m2=[(l*sc,w) for l,w in modes]
        oute.append((tc(m2,0.13,10.0)/t0-1)*100)
    # soft modes below 30 meV
    lsoft=sum(l for l,w in modes if w<30); Ssoft=sum(l*w*w for l,w in modes if w<30)/S
    mhard=[(l,w) for l,w in modes if w>=30]
    tdel=tc(mhard,0.13,10.0*max(w for _,w in modes)/max(w for _,w in mhard)) if len(mhard) else 0
    lam_hard=sum(l for l,w in mhard)
    # soft-mode force-constant error
    outs=[]
    for sc in (1.3,1/1.3):
        m2=[((l/sc, w*math.sqrt(sc)) if w<30 else (l,w)) for l,w in modes]
        outs.append((tc(m2,0.13,10.0)/t0-1)*100)
    # derivative slope
    xs=[];ys=[]
    for sc in np.geomspace(1/1.3,1.3,9):
        m2=[(l/sc, w*math.sqrt(sc)) for l,w in modes]
        t=tc(m2,0.13,10.0)
        if t>0: xs.append(math.log(lam/sc)); ys.append(math.log(t))
    # PdH falls below the temperature floor of the solver when softened; no slope is reported for it
    slope=np.polyfit(xs,ys,1)[0] if len(xs)==9 else float('nan')
    f=lambda x: tc([(x,1.0/math.sqrt(x))],0.13,10.0)
    dsm=lambda L:(math.log(f(L*1.02))-math.log(f(L/1.02)))/(2*math.log(1.02))
    print(f"{name}: lam={lam:.2f} k-scale {out[0]:.1f}/{out[1]:.1f}%  eta-scale {oute[0]:.1f}/{oute[1]:.1f}%  soft lam={lsoft:.2f} ({lsoft/lam*100:.0f}%) S share {Ssoft*100:.1f}%  delete soft: {(tdel/t0-1)*100:.1f}%  soft k err {outs[0]:.2f}/{outs[1]:.2f}%  lam>30meV={lam_hard:.2f} slope={slope:.2f} single-mode slope at lam {dsm(lam):.2f} at lam_hard {dsm(lam_hard) if lam_hard>0.5 else float('nan'):.2f}")