"""Lab to Code — How many trays does a distillation column need? (benzene + toluene, Fenske equation)

Relative volatility: alpha = Psat_benzene / Psat_toluene (Raoult's law, ideal mixture). Psat from the Antoine
equation (NIST WebBook, P in bar, T in K): log10(Psat) = A - B / (T + C), the same constants as our Raoult video.
Fenske: minimum number of stages at total reflux
    N_min = ln[(xD / (1 - xD)) * ((1 - xB) / xB)] / ln(alpha)
"""
from math import log, sqrt

ANT = {"benzene": (4.01814, 1203.835, -53.226), "toluene": (4.07827, 1343.943, -53.773)}
P = 1.01325  # bar


def psat(name, T):
    A, B, C = ANT[name]
    return 10 ** (A - B / (T + C))


def boiling_T(x, lo=340.0, hi=390.0):        # bisection, as in the Raoult video
    for _ in range(50):
        mid = (lo + hi) / 2
        total = x * psat("benzene", mid) + (1 - x) * psat("toluene", mid)
        lo, hi = (lo, mid) if total == max(total, P) else (mid, hi)
    return mid


def fenske(xD, xB, a):
    return log((xD / (1 - xD)) * ((1 - xB) / xB)) / log(a)


Tt, Tb = boiling_T(0.95), boiling_T(0.05)
a = sqrt(psat("benzene", Tt) / psat("toluene", Tt) * psat("benzene", Tb) / psat("toluene", Tb))
print(f"top {Tt - 273.15:.1f} C, bottom {Tb - 273.15:.1f} C, alpha {a:.2f}")
for xD in (0.95, 0.99):
    print(f"{xD:.0%} pure: {fenske(xD, 1 - xD, a):.1f} stages min")