r"""
Clausius's inequality and entropy (1854-1865)
=============================================

Rudolf Clausius (1854) showed that around any cycle the heat received,
each part divided by the temperature at which it crosses the boundary,
sums to zero or less: :math:`\oint \delta Q/T \le 0`
(:func:`~engineeringkit.thermodynamics.clausius_inequality`), with equality
only for a reversible cycle. Then :math:`\delta Q_{rev}/T` is the change of
a state function, which in 1865 he named *entropy*. The example builds a
Carnot cycle from an ideal gas, checks that its Clausius sum is zero, and
shows that drawing the same heat across a finite temperature difference
makes the sum negative.
"""

# %%
import matplotlib.pyplot as plt
import numpy as np

from engineeringkit.thermodynamics import (
    clausius_inequality,
    entropy_change_isothermal_ideal_gas,
    isentropic_pressure_volume,
    isentropic_temperature_ratio,
    reversible_isothermal_work_ideal_gas,
)

# %%
# A Carnot cycle of one mole of air
# ---------------------------------
#
# Isothermal expansion at 600 K from 1 to 3 litres, adiabatic expansion to
# 300 K, isothermal compression, adiabatic compression back. On the
# adiabats :math:`TV^{\gamma-1}` is constant.

n, gamma, T_h, T_c = 1.0, 1.4, 600.0, 300.0
V1, V2 = 1e-3, 3e-3
V3 = V2 * (T_h / T_c) ** (1 / (gamma - 1))
V4 = V1 * (T_h / T_c) ** (1 / (gamma - 1))
Q_h = reversible_isothermal_work_ideal_gas(n, T_h, V1, V2)
Q_c = reversible_isothermal_work_ideal_gas(n, T_c, V3, V4)  # negative: heat leaves
dS_h = entropy_change_isothermal_ideal_gas(n, T_h, V1, V2)
dS_c = entropy_change_isothermal_ideal_gas(n, T_c, V3, V4)
print(f"Q_h = {Q_h:.1f} J, Q_c = {Q_c:.1f} J; sum Q/T = {Q_h / T_h + Q_c / T_c:.2e} J/K")
print(f"entropy change on the hot isotherm {dS_h:.4f} J/K = Q_h/T_h; on the cold one {dS_c:.4f} J/K")
assert clausius_inequality([Q_h / T_h, Q_c / T_c])
assert np.isclose(Q_h / T_h, dS_h) and np.isclose(dS_h + dS_c, 0)
# The adiabats are consistent with the pressure-ratio form of the isentropic law.
p2, p3 = 8.314462618 * T_h / V2, 8.314462618 * T_c / V3
assert np.isclose(isentropic_pressure_volume(p2, V2, V3, gamma), p3)
assert np.isclose(isentropic_temperature_ratio(p3 / p2, gamma), T_c / T_h)

# %%
# The cycle in the T-S plane
# --------------------------
#
# In entropy coordinates the Carnot cycle is a rectangle, and its area is
# the net work.

S = [0, dS_h, dS_h, 0, 0]
T = [T_h, T_h, T_c, T_c, T_h]
fig, ax = plt.subplots(figsize=(5, 3.2))
ax.plot(S, T, "o-")
ax.fill(S, T, alpha=0.2)
ax.set(xlabel="entropy change [J/K]", ylabel="temperature [K]", title="Carnot cycle: area = net work")
fig.tight_layout()
assert np.isclose(dS_h * (T_h - T_c), Q_h + Q_c)

# %%
# An irreversible cycle
# ---------------------
#
# Now the same engine takes :math:`Q_h` from a furnace at 800 K: the heat
# crosses the boundary at 800 K, but the gas still receives it at 600 K.
# The Clausius sum, taken at the boundary temperature, becomes negative.

for T_furnace in (600.0, 700.0, 800.0, 1000.0):
    total = Q_h / T_furnace + Q_c / T_c
    print(f"furnace at {T_furnace:6.0f} K: sum dQ/T = {total:+.4f} J/K, holds: {clausius_inequality([Q_h / T_furnace, Q_c / T_c])}")
    assert clausius_inequality([Q_h / T_furnace, Q_c / T_c])
assert not clausius_inequality([Q_h / T_h, Q_c / 400.0])  # rejecting heat at 400 K while claiming 300 K work is impossible

plt.show()