"""f13_pair_potential -- Lennard-Jones u(r) and the Mayer function
f(r) = exp(-u/kT) - 1 at two temperatures, twin axes; mark sigma and
epsilon. u*(r*) = 4*(r*^-12 - r*^-6), reduced units (eps=sigma=1).
"""
import sys
sys.path.insert(0, "/home/claude/eos_lecture/code")
from eosstyle import *

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

r = np.linspace(0.85, 3.0, 800)
u = u_star(r)

fig, ax1 = newfig(7.0, 4.6)
ax2 = ax1.twinx()

l1, = ax1.plot(r, u, color=NAVY, lw=2.2, zorder=4)
ax1.axhline(0, color=GRAPHITE, lw=1.0, alpha=0.6)

# mark sigma (u=0 crossing at r*=1) and epsilon (well depth at r*=2^(1/6))
r_min = 2.0**(1.0/6.0)
u_min = u_star(r_min)
ax1.plot([1.0], [0.0], marker="o", ms=6.5, color=AMBER, zorder=6, markeredgecolor="white", markeredgewidth=0.9)
ax1.annotate(r"$\sigma$: $u(\sigma)=0$", xy=(1.0, 0.0), xytext=(1.36, 1.80),
             fontsize=11, color=AMBER, va="center",
             arrowprops=dict(arrowstyle="-", color=AMBER, lw=1.1))
ax1.plot([r_min], [u_min], marker="o", ms=6.5, color=RUST, zorder=6, markeredgecolor="white", markeredgewidth=0.9)
ax1.annotate(r"$-\varepsilon$: well minimum at $r=2^{1/6}\sigma$", xy=(r_min, u_min),
             xytext=(1.52, -0.58), fontsize=11, color=RUST, va="center",
             arrowprops=dict(arrowstyle="-", color=RUST, lw=1.1))

ax1.text(2.55, 0.35, "$u(r)$", color=NAVY, fontsize=12, weight="bold")

Tstars = [0.8, 3.42]
colors_f = [TEAL, SKYBLUE]
for Tst, col in zip(Tstars, colors_f):
    f = np.exp(-u/Tst) - 1.0
    l, = ax2.plot(r, f, color=col, lw=2.0, linestyle="--", zorder=3)
    # direct label near the trough / tail of each Mayer function curve
    if Tst < 1.0:
        idx = np.argmin(f)
        ax2.annotate(f"$f(r)$, $T^*={Tst:.2f}$", xy=(r[idx], f[idx]), xytext=(2.08, -0.95),
                     fontsize=11, color=col, va="center",
                     arrowprops=dict(arrowstyle="-", color=col, lw=0.9))
    else:
        # this label used to sit on the u(r) = 0 line and overran the right
        # spine; it now lives in the empty upper-right block with a leader
        # back to the descending flank of its own Mayer curve.
        idx = np.argmin(np.abs(r-1.30))
        ax2.annotate(f"$f(r)$, $T^*={Tst:.2f}$", xy=(r[idx], f[idx]), xytext=(2.00, 1.05),
                     fontsize=11, color=col, ha="left", va="center",
                     arrowprops=dict(arrowstyle="-", color=col, lw=0.9))

ax2.axhline(0, color=SLATE, lw=0.8, linestyle=":", alpha=0.7)

ax1.set_xlabel(r"Separation, $r/\sigma$")
ax1.set_ylabel(r"Pair potential, $u(r)/\varepsilon$", color=NAVY)
ax2.set_ylabel(r"Mayer function, $f(r)=e^{-u(r)/k_BT}-1$", color=GRAPHITE)
ax1.tick_params(axis="y", colors=NAVY)
ax2.tick_params(axis="y", colors=GRAPHITE)
ax1.set_xlim(0.85, 3.0)
ax1.set_ylim(-1.4, 2.5)
ax2.set_ylim(-1.15, 2.05)
ax2.grid(False)

save(fig, "f13_pair_potential")
