 codex/rc-sinusoidal-turn-on.py

"""Original switched-sinusoid RC response and a specified settling threshold."""
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


def main():
    tau = 1e-3
    omega = 1000.0
    amplitude = 1.0
    ratio = omega * tau
    t = np.linspace(0, 12 * tau, 2401)
    steady = amplitude * (np.cos(omega * t) + ratio * np.sin(omega * t)) / (1 + ratio**2)
    transient = -amplitude * np.exp(-t / tau) / (1 + ratio**2)
    exact = steady + transient
    derivative = amplitude * (-omega * np.sin(omega * t)
                              + ratio * omega * np.cos(omega * t)
                              + np.exp(-t / tau) / tau) / (1 + ratio**2)
    assert exact[0] == 0
    assert np.allclose(tau * derivative + exact, amplitude * np.cos(omega * t))
    relative = np.exp(-t / tau) / np.sqrt(1 + ratio**2)
    tolerance = 0.01
    settling = tau * np.log(1 / (tolerance * np.sqrt(1 + ratio**2)))
    assert np.isclose(settling / tau, 4.258596595708118)
    assert np.allclose(abs(transient) / (amplitude / np.sqrt(1 + ratio**2)), relative)

    fig, (response, decay) = plt.subplots(1, 2, figsize=(9.6, 4.2), dpi=100,
                                        layout="constrained")
    fig.set_facecolor("white")
    response.plot(t * 1000, steady, "--", color="#286f99", label="Phasor steady state")
    response.plot(t * 1000, exact, color="#b66b31", label="Exact switch-on response")
    response.scatter([0, 0], [0.5, 0], color=["#286f99", "#b66b31"], s=22, zorder=4)
    response.set(ylabel="Capacitor potential difference (V)",
                 title="The initial condition requires a transient", ylim=(-0.85, 1.15))
    response.legend(fontsize=8.5, loc="upper right")
    decay.semilogy(t * 1000, 100 * relative, color="#286f99", label="Relative transient")
    decay.axhline(100 * tolerance, color="#444444", linestyle="--", label="1% allowance")
    decay.axvline(settling * 1000, color="#b66b31", linestyle=":",
                   label=f"Settled after {settling * 1000:.3f} ms")
    decay.set(ylabel="Transient / steady-state amplitude (%)", ylim=(1e-4, 120),
              title="An exponential error bound sets the wait")
    decay.legend(fontsize=8.5, loc="upper right")
    for ax in [response, decay]:
        ax.set(xlabel="Time after switch-on (ms)", xlim=(0, 12))
        ax.grid(alpha=0.15)
    fig.savefig(Path(__file__).with_suffix(".png").name, facecolor="white")
    plt.close(fig)


if __name__ == "__main__":
    main()