r"""
Swing curves and the equal-area criterion (1926)
================================================

As transmission lines grew longer in the 1920s, generators at their two
ends began to fall out of step after faults. Westinghouse engineers
R. D. Evans and C. F. Wagner computed *swing curves*, the rotor angle
:math:`\delta` against time, by integrating the swing equation
:math:`(2H/\omega_s)\,\ddot\delta = P_m - P_{max}\sin\delta` step by step
(:func:`~engineeringkit.electrical_engineering.swing_response`). For one
machine against an infinite bus there is a shortcut: the kinetic energy
gained while the fault holds the electrical power down (the accelerating
area) must be given back after clearing (the decelerating area) before
:math:`\delta` passes the unstable equilibrium. Equating the two areas
gives the critical clearing angle
(:func:`~engineeringkit.electrical_engineering.equal_area_clearing_angle`)
and, for a bolted fault, the critical clearing time
(:func:`~engineeringkit.electrical_engineering.critical_clearing_time`).
This example checks the criterion against the integrated swing curves.
"""

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

from engineeringkit.electrical_engineering import critical_clearing_time, equal_area_clearing_angle, swing_response

# %%
# The critical clearing angle and time
# ------------------------------------
#
# A 60 Hz generator with :math:`H = 5` s delivers 0.9 pu over a double line
# (:math:`P_{max} = 2.0` pu). A fault on one line drops the transfer to
# 0.5 pu, and clearing it trips that line, leaving 1.5 pu.

H, Pm, P1, P2, P3 = 5.0, 0.9, 2.0, 0.5, 1.5
d0, dcr = equal_area_clearing_angle(Pm, P1, P2, P3)
dmax = np.pi - np.arcsin(Pm / P3)
print(f"delta0 = {np.degrees(d0):.1f} deg, critical clearing angle {np.degrees(dcr):.1f} deg")
print(f"unstable equilibrium {np.degrees(dmax):.1f} deg")

# The two areas are equal at the critical angle.
accel = Pm * (dcr - d0) + P2 * (np.cos(dcr) - np.cos(d0))
decel = P3 * (np.cos(dcr) - np.cos(dmax)) - Pm * (dmax - dcr)
assert np.isclose(accel, decel)

# %%
# Finding the critical time from the swing curves
# -----------------------------------------------
#
# With power still flowing during the fault there is no closed form for the
# time, so bisect on the clearing time until the first swing just reaches
# the unstable equilibrium; at that clearing time the angle at clearing is
# the critical angle.

low, high = 0.05, 1.0
for _ in range(40):
    mid = 0.5 * (low + high)
    low, high = (mid, high) if swing_response(H, Pm, P1, P2, P3, mid, 3.0).stable else (low, mid)
r = swing_response(H, Pm, P1, P2, P3, low, 3.0, n_points=6001)
angle_at_clearing = r.delta[np.searchsorted(r.t, r.clearing_time)]
print(f"critical clearing time {low:.4f} s, angle at clearing {np.degrees(angle_at_clearing):.2f} deg")
assert np.isclose(angle_at_clearing, dcr, atol=np.radians(0.1))

# %%
# A bolted three-phase fault
# --------------------------
#
# With :math:`P_2 = 0` the angle grows as a parabola and the critical time
# has a closed form.

tcr = critical_clearing_time(H, Pm, P1, frequency=60.0, p_max_post=P3)
assert swing_response(H, Pm, P1, 0.0, P3, 0.98 * tcr, 3.0).stable
assert not swing_response(H, Pm, P1, 0.0, P3, 1.02 * tcr, 3.0).stable
print(f"bolted fault: critical clearing time {tcr:.3f} s ({60 * tcr:.1f} cycles)")

fig, axes = plt.subplots(1, 2, figsize=(11, 3.9))
delta = np.linspace(0, np.pi, 300)
ax = axes[0]
for p, label in ((P1, "prefault"), (P2, "fault on"), (P3, "post-fault")):
    ax.plot(np.degrees(delta), p * np.sin(delta), label=label)
ax.axhline(Pm, color="k", lw=0.8)
fill = np.linspace(d0, dcr, 100)
ax.fill_between(np.degrees(fill), P2 * np.sin(fill), Pm, alpha=0.3, label="accelerating area")
fill = np.linspace(dcr, dmax, 100)
ax.fill_between(np.degrees(fill), Pm, P3 * np.sin(fill), alpha=0.3, label="decelerating area")
ax.set(xlabel="rotor angle [deg]", ylabel="electrical power [pu]", title="Equal areas at the critical angle")
ax.legend(fontsize="x-small")
for tc in (0.9 * low, 1.1 * low):
    s = swing_response(H, Pm, P1, P2, P3, tc, 2.0)
    axes[1].plot(s.t, np.degrees(s.delta), label=f"cleared at {tc:.3f} s ({'stable' if s.stable else 'unstable'})")
axes[1].axhline(np.degrees(dmax), color="gray", ls=":", lw=0.8)
axes[1].set(xlabel="time [s]", ylabel="rotor angle [deg]", ylim=(0, 360), title="Swing curves")
axes[1].legend(fontsize="small")
fig.tight_layout()

plt.show()