"""f12_virial_B -- B2*(T*) of the Lennard-Jones fluid computed by numerical
integration of the Mayer function (same method as f10_boyle.py), plotted in
reduced form with the high-T repulsive plateau annotated.
"""
import sys
sys.path.insert(0, "/home/claude/eos_lecture/code")
from eosstyle import *
from scipy.integrate import quad
from scipy.optimize import brentq

def u_star(r):
    return 4.0*(r**-12 - r**-6)

def mayer_integrand(r, Tstar):
    expo = -u_star(r)/Tstar
    if expo < -700:
        f = -1.0
    else:
        f = np.exp(expo) - 1.0
    return f * r**2

def B2star(Tstar, r_split=0.3, r_max=50.0):
    core = -2*np.pi * (-(r_split**3)/3.0)
    outer, _ = quad(mayer_integrand, r_split, r_max, args=(Tstar,), limit=400,
                     points=[1.0], epsabs=1e-10, epsrel=1e-9)
    outer = -2*np.pi*outer
    return core + outer

T_B = brentq(B2star, 2.5, 4.5, xtol=1e-8)
print(f"f12: computed LJ Boyle temperature T_B* = {T_B:.4f} (reference ~3.42)")

Tstars = np.geomspace(0.5, 100.0, 300)
B2vals = np.array([B2star(T) for T in Tstars])

# NOTE on high-T behaviour (verified numerically, not assumed): B2* does
# NOT approach a true hard-sphere plateau. It rises through zero at T_B*,
# passes through a broad maximum around T* ~ 20-30 (~1.10, well below the
# rigid hard-sphere value 2*pi/3 = 2.094 because the LJ repulsive core is
# soft, not rigid), and then decays slowly back toward zero at still higher
# T* as the effective core shrinks further (checked up to T*=1e5, where
# B2* ~ 0.20). The figure is annotated for what the data actually show: a
# broad, only weakly-varying repulsive region across the plotted window,
# not an asymptotic excluded-volume constant.
imax = np.argmax(B2vals)
Tmax, B2max = Tstars[imax], B2vals[imax]
print(f"f12: B2* maximum in window = {B2max:.4f} at T*={Tmax:.1f} "
      f"(hard-sphere value would be {2*np.pi/3:.3f}; LJ core is soft, not rigid)")

fig, ax = newfig(7.2, 4.6)
ax.set_xscale("log")
ax.plot(Tstars, B2vals, color=NAVY, lw=2.1)
ax.axhline(0, color=GRAPHITE, lw=1.0, alpha=0.7)

Tb_marker = brentq(B2star, 2.5, 4.5, xtol=1e-8)
ax.plot([Tb_marker], [0], marker="o", ms=7, color=AMBER, zorder=5, markeredgecolor="white", markeredgewidth=1.0)
ax.annotate(f"$T_B^*={Tb_marker:.2f}$", xy=(Tb_marker, 0), xytext=(4.8, -3.4),
            fontsize=11, color="#8a6b16", ha="left", va="center",
            arrowprops=dict(arrowstyle="-", color=AMBER, lw=1.1))

ax.plot([Tmax], [B2max], marker="o", ms=6, color=TEAL, zorder=5, markeredgecolor="white", markeredgewidth=0.8)
# the three-line note is placed in the clear band above the curve; the axis
# top is raised so that the whole block fits inside the panel.
ax.annotate(f"broad repulsive region: $B_2^*>0$\nvaries only slowly with $T^*$\n(soft core, not a rigid-sphere plateau)",
            xy=(Tmax, B2max), xytext=(1.35, 7.2), fontsize=11, color=TEAL,
            ha="left", va="top",
            linespacing=1.3, arrowprops=dict(arrowstyle="-", color=TEAL, lw=1.0))

# moved off the steeply descending branch into the empty region below the
# curve, joined back to it with a short leader.
ax.annotate("attraction dominates\nat low $T^*$", xy=(0.78, B2star(0.78)),
            xytext=(1.6, -9.0), fontsize=11, color=RUST, ha="left", va="top",
            linespacing=1.25, arrowprops=dict(arrowstyle="-", color=RUST, lw=1.0))

ax.set_xlabel(r"Reduced temperature, $T^* = k_BT/\varepsilon$ (log scale)")
ax.set_ylabel(r"Reduced second virial coeff., $B_2^*=B_2/N_A\sigma^3$")
# the y label is a touch longer than the axes; nudge it down so that its
# superscript is not shaved off by the tight bounding box on save
ax.yaxis.label.set_y(0.46)
ax.set_xlim(0.5, 100)
ax.set_ylim(-16, 7.6)

save(fig, "f12_virial_B")
