 codex/rc-pulse-train-baseline.py

"""Original ideal RC pulse train; time is expressed in units of RC."""
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


def main():
    high, quiet, threshold = 0.4, 0.6, 0.5
    period = high + quiet
    alpha, beta = np.exp(-high), np.exp(-quiet)
    q = alpha * beta
    peak_star = (1 - alpha) / (1 - q)
    base_star = beta * peak_star
    up = np.log((1 - base_star) / (1 - threshold))
    down = high + np.log(peak_star / threshold)
    assert np.isclose(base_star, 0.28623051789026865)
    assert np.isclose(peak_star, 0.5215460078933705)
    assert np.isclose(down - up, 0.08623743675739143)
    bases = base_star * (1 - q ** np.arange(6))
    peaks = 1 - alpha + alpha * bases
    assert np.all(peaks[:3] < threshold) and peaks[3] > threshold
    assert np.allclose(bases[1:], beta * peaks[:-1])
    assert np.isclose(1 - (1 - base_star) * np.exp(-up), threshold)
    assert np.isclose(peak_star * np.exp(-(down - high)), threshold)

    time = np.linspace(0, 5 * period, 10001)
    cycle = np.floor(time / period).astype(int)
    phase = time - cycle * period
    baseline = base_star * (1 - q**cycle)
    peak = 1 - alpha + alpha * baseline
    output = np.where(phase <= high, 1 - (1 - baseline) * np.exp(-phase),
                      peak * np.exp(-(phase - high)))
    steady_phase = np.linspace(0, period, 2001)
    steady = np.where(steady_phase <= high,
                      1 - (1 - base_star) * np.exp(-steady_phase),
                      peak_star * np.exp(-(steady_phase - high)))
    fig, (startup, settled) = plt.subplots(1, 2, figsize=(10, 4.3), dpi=100,
                                         layout="constrained")
    fig.set_facecolor("white")
    startup.plot(time, (phase < high).astype(float), color="#999999", linewidth=1,
                 label="Applied pulse train")
    startup.plot(time, output, color="#286f99", linewidth=1.8, label="Capacitor output")
    startup.axhline(threshold, color="#9f4b27", linestyle="--", label="Threshold")
    startup.scatter([3 * period + high], [peaks[3]], color="#286f99", s=30, zorder=5)
    startup.annotate("Fourth pulse first crosses", (3 * period + high, peaks[3]),
                     xytext=(1.3, 0.74), fontsize=9,
                     arrowprops={"arrowstyle": "->", "color": "#286f99"})
    startup.set(xlim=(0, 5), ylim=(-0.03, 1.06), xlabel="Time / τ",
                 ylabel="Electric potential difference / V₀",
                 title="Repeated short pulses change their own baseline")
    startup.legend(loc="upper left", fontsize=8)
    settled.axvspan(0, high, color="#286f99", alpha=0.06)
    settled.plot(steady_phase, steady, color="#286f99", linewidth=2,
                 label="Repeatable response")
    settled.axhline(threshold, color="#9f4b27", linestyle="--", label="Threshold")
    settled.scatter([up, down], [threshold, threshold], color="#9f4b27", s=30, zorder=5)
    settled.annotate("Above threshold for 0.0862τ", ((up + down) / 2, threshold),
                     xytext=(0.42, 0.45), fontsize=9,
                     arrowprops={"arrowstyle": "->", "color": "#9f4b27"})
    settled.text(0.03, 0.25, "Baseline = 0.28623 V₀", fontsize=9)
    settled.set(xlim=(0, period), ylim=(0.22, 0.58), xlabel="Time within one period / τ",
                 ylabel="Electric potential difference / V₀",
                 title="A nonzero low interval still permits residual charge")
    settled.legend(loc="upper right", fontsize=8)
    for ax in (startup, settled):
        ax.grid(alpha=0.14)
    fig.savefig(Path(__file__).with_suffix(".png").name, facecolor="white")
    plt.close(fig)


if __name__ == "__main__":
    main()