"""f10_boyle -- second virial coefficient B2*(T*) of a Lennard-Jones fluid,
computed by numerical integration of the Mayer function:

    B2*(T*) = -2*pi * integral_0^inf [exp(-u*(r*)/T*) - 1] r*^2 dr*
    u*(r*) = 4*(r*^-12 - r*^-6)   (reduced LJ potential, eps=sigma=1)

The r*->0 singularity is handled by splitting the integral at r*=0.3: for
r* < 0.3, u*/T* is enormous and exp(-u*/T*) underflows to 0, so the Mayer
function is -1 there and that piece integrates analytically as
-2*pi*integral_0^0.3 (0-1) r*^2 dr* = +2*pi*(0.3^3)/3. For r* >= 0.3 we
integrate exp(-u*/T*)-1 numerically out to r*=50 (negligible tail beyond).

Boyle temperature T_B* is where B2*(T_B*) = 0; the standard LJ result is
T_B* ~= 3.42 (e.g. Hirschfelder, Curtiss & Bird). We solve for it and print
the value; the code asserts it lands close to that reference.

Right panel: initial slope dZ/dP at P->0 (equal to B2*/T* in reduced form)
changes sign at T_B*, shown as Z vs P* for a few reduced temperatures using
the low-density expansion Z = 1 + B2*(T*) rho*, P* = rho* T* Z.
"""
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:      # underflow guard, exp(-u/T) -> 0
        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 contribution (r < r_split): Mayer function is -1 there (hard-core-like)
    core = -2*np.pi * (-(r_split**3)/3.0)   # = +2*pi*r_split^3/3
    # outer contribution, numerical
    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

# ---- verify Boyle temperature -----------------------------------------
T_B = brentq(B2star, 2.5, 4.5, xtol=1e-8)
print(f"Computed LJ Boyle temperature T_B* = {T_B:.4f} (reference ~3.42)")
assert abs(T_B - 3.42) < 0.05, "Boyle temperature integration is off; check quadrature"

Tstars = np.linspace(0.55, 12.0, 260)
B2vals = np.array([B2star(T) for T in Tstars])

# The two panels are shown about 6 in wide on the slide, so the canvas is kept
# small and the type is enlarged: what matters for legibility is the ratio of
# font size to figure width, not the font size on its own.
plt.rcParams.update({
    "xtick.labelsize": 13,
    "ytick.labelsize": 13,
    "axes.labelsize": 14,
})
FS = 12.0        # in-plot annotations

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8.0, 4.3))
for ax in (ax1, ax2):
    ax.grid(True, zorder=0)
    ax.set_axisbelow(True)

ax1.plot(Tstars, B2vals, color=NAVY, lw=2.1)
ax1.axhline(0, color=GRAPHITE, lw=1.0, alpha=0.7)
ax1.plot([T_B], [0], marker="o", ms=7, color=AMBER, zorder=5, markeredgecolor="white", markeredgewidth=1.0)

# the three captions are put on three different levels (top strip, and two
# rows in the large empty area under the curve) so that no two of them share
# a baseline; each keeps a short leader to its own feature on the curve
ax1.annotate("repulsion dominates\n($B_2^*>0$)", xy=(9.6, B2star(9.6)), xytext=(11.9, 5.2),
             fontsize=FS, color=TEAL, ha="right", va="top", linespacing=1.3,
             arrowprops=dict(arrowstyle="-", color=TEAL, lw=1.0))
ax1.annotate("Boyle temperature\n$T_B^*=3.42$", xy=(T_B, 0), xytext=(4.4, -2.6),
             fontsize=FS, color="#8a6b16", ha="left", va="top", linespacing=1.3,
             arrowprops=dict(arrowstyle="-", color=AMBER, lw=1.1))
ax1.annotate("attraction dominates\n($B_2^*<0$)", xy=(1.3, B2star(1.3)), xytext=(2.1, -9.2),
             fontsize=FS, color=RUST, ha="left", va="top", linespacing=1.3,
             arrowprops=dict(arrowstyle="-", color=RUST, lw=1.0))

ax1.set_xlabel(r"Reduced temperature, $T^* = k_BT/\varepsilon$")
ax1.set_ylabel("Reduced second virial coeff.,\n" r"$B_2^*=B_2/N_A\sigma^3$", linespacing=1.4)
ax1.set_xlim(0.5, 12)
ax1.set_ylim(-16.5, 5.5)
ax1.set_xticks([2, 4, 6, 8, 10, 12])
ax1.set_yticks([-15, -10, -5, 0, 5])

# ---- right panel: initial slope of Z vs P changing sign at T_B --------
# fixed reduced-density grid for every temperature (not P-normalised), so
# curves naturally fan out; direct labels are placed at each curve's own
# endpoint with a minimum vertical pitch enforced to avoid overlap.
rho_star = np.linspace(0.0, 0.06, 60)
T_show = [1.5, 2.5, T_B, 5.0, 9.0]
colors_show = [RUST, "#c98a6e", AMBER, SKYBLUE, TEAL]

ends = []
for Tst, col in zip(T_show, colors_show):
    B2 = B2star(Tst)
    Z = 1.0 + B2*rho_star
    Pstar = rho_star*Tst*Z
    ax2.plot(Pstar, Z, color=col, lw=2.0)
    ends.append((Tst, col, Pstar[-1], Z[-1]))

ax2.axhline(1.0, color=GRAPHITE, lw=1.2, linestyle=":", zorder=2)
# the ideal-gas note goes to the far right, past the end of the hottest curve,
# where nothing else is drawn
ax2.text(0.50, 1.012, "$Z=1$ (ideal gas)", fontsize=FS, color=GRAPHITE, va="bottom", ha="left")

# End labels. Three curves are labelled straight off their own end point; the
# T*=T_B curve runs along Z=1 and the T*=5 curve is hidden under the T*=9 one,
# so those two get a short leader into clear space instead.
ends_by_T = {round(e[0], 3): e for e in ends}
plain = {1.5: (0.042, 0.006), 2.5: (0.012, 0.0), 9.0: (0.015, 0.0)}
for Tst, col, Px, Z in ends:
    key = round(Tst, 3)
    if key in plain:
        dx, dy = plain[key]
        ax2.text(Px + dx, Z + dy, f"$T^*={Tst:.2f}$", color=col, fontsize=FS,
                 va="center", ha="left")
    elif abs(Tst - 5.0) < 1e-6:
        ax2.annotate(f"$T^*={Tst:.2f}$", xy=(Px, Z), xytext=(0.10, 1.062),
                     color=col, fontsize=FS, va="top", ha="left",
                     arrowprops=dict(arrowstyle="-", color=col, lw=0.9, alpha=0.85))
    else:   # T_B: sits exactly on the Z = 1 line
        ax2.annotate(f"$T^*={Tst:.2f}$", xy=(Px, Z), xytext=(0.425, 0.968),
                     color=col, fontsize=FS, va="center", ha="left",
                     arrowprops=dict(arrowstyle="-", color=col, lw=0.9, alpha=0.85,
                                     shrinkA=14, shrinkB=3))

ax2.set_xlabel(r"Reduced pressure, $P^*=P\sigma^3/\varepsilon$")
ax2.set_ylabel("Compressibility factor, $Z$\n(low density)", linespacing=1.4)
ax2.set_xlim(0, 0.80)
ax2.set_ylim(0.835, 1.075)

fig.subplots_adjust(left=0.085, right=0.99, wspace=0.34, bottom=0.155, top=0.97)
save(fig, "f10_boyle")
