# kring.py - arbeid en rendement van kringprocessen: numeriek W = ∮ p dV, en de formule
import numpy as np

R, g = 8.314, 1.4                       # J/(mol K); lucht: cp/cv


def arbeid(takken):
    """Som van de trapeziumregel over alle takken: elke tak is (V, p) in m³ en Pa, in volgorde doorlopen."""
    return sum(np.trapezoid(p, V) for V, p in takken)


def isotherm(n, T, Va, Vb):
    V = np.geomspace(Va, Vb, 2000)
    return V, n * R * T / V


def adiabaat(p0, V0, Vb):
    V = np.geomspace(V0, Vb, 2000)
    return V, p0 * (V0 / V) ** g


# Carnot: 0,2005 mol tussen 600 en 300 K, V van 1 naar 2 liter op de hete isotherm
n, Th, Tc, V1, V2 = 0.2005, 600.0, 300.0, 1e-3, 2e-3
V3, V4 = V2 * (Th / Tc) ** (1 / (g - 1)), V1 * (Th / Tc) ** (1 / (g - 1))
t1 = isotherm(n, Th, V1, V2)
t2 = adiabaat(n * R * Th / V2, V2, V3)
t3 = isotherm(n, Tc, V3, V4)
t4 = adiabaat(n * R * Tc / V4, V4, V1)
W = arbeid([t1, t2, t3, t4])
Qh = n * R * Th * np.log(V2 / V1)
print(f"Carnot:  W = {W:6.1f} J, Qh = {Qh:6.1f} J, rendement {W / Qh:.4f} (formule {1 - Tc / Th:.4f})")

# Otto: r = 10, 1 bar en 300 K in, 2500 K na de verbranding
r, Vmax, T1, T3 = 10.0, 0.5556e-3, 300.0, 2500.0
n = 1e5 * Vmax / (R * T1)
V2 = Vmax / r
T2 = T1 * r ** (g - 1)
p2, p3 = n * R * T2 / V2, n * R * T3 / V2
W = arbeid([adiabaat(1e5, Vmax, V2), adiabaat(p3, V2, Vmax)])      # de isochoren leveren geen arbeid
Qin = n * 2.5 * R * (T3 - T2)
print(f"Otto:    W = {W:6.1f} J, Qin = {Qin:6.1f} J, rendement {W / Qin:.4f} (formule {1 - r ** (1 - g):.4f})")

# Stirling: twee isothermen en twee isochoren, 0,2005 mol tussen 600 en 300 K, V van 1 naar 2 liter
n, V1, V2 = 0.2005, 1e-3, 2e-3
W = arbeid([isotherm(n, Th, V1, V2), isotherm(n, Tc, V2, V1)])
Qiso = n * R * Th * np.log(V2 / V1)
Qv = n * 2.5 * R * (Th - Tc)                  # opwarmen bij constant volume
print(f"Stirling zonder regenerator: rendement {W / (Qiso + Qv):.4f}")
print(f"Stirling met regenerator:    rendement {W / Qiso:.4f}")