Skip to content

Transient simulations

The groundinsight.simulation.transient module extends the package beyond the steady-state phasor solve into the time domain. A TransientStudy binds a network, an active fault and one or more user-defined source waveforms together and produces a ResultTransient with EPR and shield-current time series at the declared observation points.

Physical / modelling context

The frequency-domain solver answers "given a sinusoidal injection at frequency \(f\), what is the EPR?". Real fault currents are non-sinusoidal: they switch on at fault inception, may carry an exponentially decaying DC offset, and switch off again at clearing. Two complementary solver paths are implemented to capture that behaviour:

  • FFT solver (solver="fft") — samples the user waveform on a regular time grid, transforms to the frequency domain via NumPy's real-valued FFT, evaluates the existing nodal-admittance solve at every FFT bin, and transforms the bus voltages back via IFFT. It reuses BusType.impedance_formula and BranchType.self_impedance_formula and is therefore consistent with the stationary results bin-by-bin. Only current sources are accepted, and mutual coupling is not evaluated by this path.
  • State-space solver (solver="state_space") — assembles a modified-nodal-analysis ODE system \(\dot x = A x + B u\), \(y = C x + D u\) from the lumped RLC fields on BusType and BranchType (R_formula, L_formula, C_formula, R_self_formula, L_self_formula, C_self_formula, R_mutual_formula, M_mutual_formula) and integrates with scipy.signal.lsim. Voltage sources, Carson-style mutual coupling and pi-section branch capacitance are supported.

Source waveforms are produced by the small library in groundinsight.simulation.waveforms: step, sinusoidal_with_dc_offset (the textbook single-line-to-ground fault current with DC asymmetry) and damped_oscillation. Custom waveforms are any vectorised callable f(t) -> values.

Example

import groundinsight as gi
from groundinsight import waveforms

# Assume `net` is a built network with one current source 'infeed'
# and a fault 'fault1'.

study = gi.TransientStudy(network=net, fault_name="fault1")
study.set_source_waveform(
    "infeed",
    waveforms.sinusoidal_with_dc_offset(
        amplitude=1e3, frequency_hz=50.0,
        t_on=0.02, t_off=0.12,
        dc_amplitude=500.0, dc_decay_tau=0.05,
    ),
)
study.set_observation(buses=["bus_fault"], branches=["cable_1"])

result = study.solve(t_end=0.2, dt=1e-4, solver="fft")

# Plot the time series
gi.plot_epr_transient(result=result, title="EPR transient")
gi.plot_branch_current_transient(result=result, title="Shield current")

# Long-format DataFrame for further post-processing
df = result.to_polars()

Switching to the state-space solver only requires changing the solver argument and ensuring the network's bus and branch types carry the lumped RLC formulas required by the ODE form.

API reference

Transient study and result

transient

Transient Simulation Layer.

This module hosts the high-level :class:TransientStudy workflow and the matching :class:ResultTransient Pydantic model. The first solver path implemented here is FFT-based: a user-defined source waveform is sampled on a regular time grid, transformed to a frequency spectrum via NumPy's real-valued FFT, fed into the existing per-frequency network solve, and transformed back via IFFT. The state-space solver path is reserved for the next release; it will reuse the same study object and observation contract.

Design choices recorded in the Phase 3 discussion:

  • Default and currently only supported source mode for the FFT solver is the legacy current source (Source.source_type='current'). A current source's frequency-dependent values are replaced by the FFT spectrum of the user-supplied waveform; the rest of the network remains untouched. Voltage-mode sources will be supported by the state-space solver in Phase 4 (where the loop closure is naturally part of the ODE system).
  • Observation points (buses for EPR, branches for shield current) are passed explicitly to set_observation -- there is no "all" default, to keep memory and post-processing tractable on large networks.
  • Mutual coupling is not yet evaluated by the FFT solver; the phase current along each path would itself be time-dependent and is more natural to handle in the state-space formulation. The FFT solver ignores mutual impedance and emits no warning -- it is documented as a known limitation and the demo notebook shows how to interpret the resulting EPR.

ResultTransient

Bases: BaseModel

Container for the time-domain results of a transient simulation.

Attributes:

Name Type Description
time_s list of float

Time samples in seconds, equally spaced.

epr_t dict of str to list of float

Mapping of observed bus name to its EPR time series in volts.

i_branch_t dict of str to list of float

Mapping of observed branch name to its branch current time series in amperes.

source_t dict of str to list of float

The sampled source waveforms keyed by source name, in the natural unit of the source (amperes for current sources, volts for voltage sources).

fault str

Name of the fault that was active during the simulation.

solver str

Identifier of the solver that produced this result ("fft" or "state_space").

Notes

The list-typed fields use plain Python lists for JSON serialisability. Convert to numpy.ndarray at the call site if you need vectorised post-processing.

to_polars

to_polars() -> pl.DataFrame

Convert the result to a Polars DataFrame in long form.

Returns:

Type Description
DataFrame

A DataFrame with columns time_s, signal_kind ("epr" / "i_branch" / "source"), name and value.

Source code in src/groundinsight/simulation/transient.py
def to_polars(self) -> "pl.DataFrame":  # noqa: F821 -- forward ref
    """
    Convert the result to a Polars DataFrame in long form.

    Returns
    -------
    polars.DataFrame
        A DataFrame with columns ``time_s``, ``signal_kind``
        (``"epr"`` / ``"i_branch"`` / ``"source"``), ``name`` and
        ``value``.
    """
    import polars as pl

    rows = []
    for bus_name, series in self.epr_t.items():
        for t, v in zip(self.time_s, series):
            rows.append(
                {
                    "time_s": t,
                    "signal_kind": "epr",
                    "name": bus_name,
                    "value": v,
                }
            )
    for branch_name, series in self.i_branch_t.items():
        for t, v in zip(self.time_s, series):
            rows.append(
                {
                    "time_s": t,
                    "signal_kind": "i_branch",
                    "name": branch_name,
                    "value": v,
                }
            )
    for src_name, series in self.source_t.items():
        for t, v in zip(self.time_s, series):
            rows.append(
                {
                    "time_s": t,
                    "signal_kind": "source",
                    "name": src_name,
                    "value": v,
                }
            )
    return pl.DataFrame(rows)

TransientStudy

TransientStudy(network: Network, fault_name: str)

High-level entry point for transient simulations.

A study binds together a :class:Network, the active fault to study, a set of source waveforms (one per source contributing to the fault), and the list of observation points (buses and branches) that should be returned as time series. Calling :meth:solve produces a :class:ResultTransient.

Parameters:

Name Type Description Default
network Network

The network model. Buses, branches, sources and the active fault are taken from this object. The network is not mutated.

required
fault_name str

Name of the fault to activate for this study.

required

Examples:

>>> import groundinsight as gi
>>> from groundinsight.simulation import waveforms
>>> study = gi.TransientStudy(network, fault_name='F1')
>>> study.set_source_waveform(
...     'src',
...     waveforms.sinusoidal_with_dc_offset(
...         amplitude=1e3, frequency_hz=50.0,
...         t_on=0.02, t_off=0.12,
...         dc_amplitude=500.0, dc_decay_tau=0.05,
...     ),
... )
>>> study.set_observation(buses=['bus_fault'], branches=['line1'])
>>> result = study.solve(t_end=0.2, dt=1e-4)
Source code in src/groundinsight/simulation/transient.py
def __init__(self, network: Network, fault_name: str):
    if fault_name not in network.faults:
        raise ValueError(
            f"Fault '{fault_name}' does not exist in the network."
        )
    self.network = network
    self.fault_name = fault_name
    self._source_waveforms: Dict[str, WaveformFunc] = {}
    self._obs_buses: List[str] = []
    self._obs_branches: List[str] = []
    # Stand-in impedance for elements that are an ideal short circuit at
    # the 0 Hz bin, and the index of that bin. Both stay ``None`` unless
    # the FFT solver actually meets such an element. See
    # :meth:`_dc_substitute_for` and :meth:`_resolved_z`.
    self._dc_substitute: Optional[float] = None
    self._dc_bin: Optional[int] = None

set_observation

set_observation(
    *,
    buses: Optional[List[str]] = None,
    branches: Optional[List[str]] = None
)

Declare which buses and branches should be returned as time series.

Parameters:

Name Type Description Default
buses list of str

Names of buses whose EPR u(t) should be returned. Defaults to no buses.

None
branches list of str

Names of branches whose shield current i(t) should be returned. Defaults to no branches.

None

Raises:

Type Description
ValueError

If a name is not present in the network.

Source code in src/groundinsight/simulation/transient.py
def set_observation(
    self,
    *,
    buses: Optional[List[str]] = None,
    branches: Optional[List[str]] = None,
):
    """
    Declare which buses and branches should be returned as time series.

    Parameters
    ----------
    buses : list of str, optional
        Names of buses whose EPR ``u(t)`` should be returned. Defaults
        to no buses.
    branches : list of str, optional
        Names of branches whose shield current ``i(t)`` should be
        returned. Defaults to no branches.

    Raises
    ------
    ValueError
        If a name is not present in the network.
    """
    for name in buses or []:
        if name not in self.network.buses:
            raise ValueError(f"Bus '{name}' not in network.")
    for name in branches or []:
        if name not in self.network.branches:
            raise ValueError(f"Branch '{name}' not in network.")
    self._obs_buses = list(buses or [])
    self._obs_branches = list(branches or [])

set_source_waveform

set_source_waveform(
    source_name: str, waveform: WaveformFunc
)

Bind a time-domain waveform to a network source.

For source_type='current' the waveform is the injected current in amperes. For source_type='voltage' the waveform is the Thevenin EMF in volts (the source_impedance is taken from the :class:Source definition). Voltage sources are accepted by the state-space solver but rejected by the FFT solver, which still only supports current sources — the relevant solver checks the type at solve() time.

Parameters:

Name Type Description Default
source_name str

Name of a source defined on the network.

required
waveform callable

Vectorised function f(t) -> values mapping a 1-D time array to source values of the same shape.

required

Raises:

Type Description
ValueError

If the source does not exist on the network.

Source code in src/groundinsight/simulation/transient.py
def set_source_waveform(self, source_name: str, waveform: WaveformFunc):
    """
    Bind a time-domain waveform to a network source.

    For ``source_type='current'`` the waveform is the injected current
    in amperes. For ``source_type='voltage'`` the waveform is the
    Thevenin EMF in volts (the ``source_impedance`` is taken from the
    :class:`Source` definition). Voltage sources are accepted by the
    state-space solver but rejected by the FFT solver, which still
    only supports current sources — the relevant solver checks the
    type at ``solve()`` time.

    Parameters
    ----------
    source_name : str
        Name of a source defined on the network.
    waveform : callable
        Vectorised function ``f(t) -> values`` mapping a 1-D time
        array to source values of the same shape.

    Raises
    ------
    ValueError
        If the source does not exist on the network.
    """
    src = self.network.sources.get(source_name)
    if src is None:
        raise ValueError(
            f"Source '{source_name}' does not exist in the network."
        )
    if src.source_type not in ("current", "voltage"):
        raise ValueError(
            f"Source '{source_name}' has unknown source_type "
            f"'{src.source_type}'."
        )
    self._source_waveforms[source_name] = waveform

solve

solve(
    *, t_end: float, dt: float, solver: str = "fft"
) -> ResultTransient

Run the transient solve and return the resulting time series.

Parameters:

Name Type Description Default
t_end float

End time of the simulation in seconds. t=0 is always the start.

required
dt float

Time-step in seconds. For the FFT solver this determines the Nyquist frequency f_max = 1/(2*dt) and therefore the highest harmonic resolved. For the state-space solver dt is the integration time-step (first-order hold on the source signal between samples).

required
solver (fft, state_space)

Solver to use. "fft" uses BusType.impedance_formula and BranchType.self_impedance_formula (frequency-domain impedances); "state_space" uses the lumped RLC fields on BusType / BranchType and integrates with :func:scipy.signal.lsim. Defaults to "fft".

'fft'

Returns:

Type Description
ResultTransient

Time-domain results at the observation points.

Raises:

Type Description
ValueError

If no source waveform is set, if the time parameters are invalid, or if the network is missing the lumped RLC fields required by the state-space solver.

NotImplementedError

If solver is not one of the supported identifiers.

Source code in src/groundinsight/simulation/transient.py
def solve(
    self,
    *,
    t_end: float,
    dt: float,
    solver: str = "fft",
) -> ResultTransient:
    """
    Run the transient solve and return the resulting time series.

    Parameters
    ----------
    t_end : float
        End time of the simulation in seconds. ``t=0`` is always the
        start.
    dt : float
        Time-step in seconds. For the FFT solver this determines the
        Nyquist frequency ``f_max = 1/(2*dt)`` and therefore the
        highest harmonic resolved. For the state-space solver ``dt``
        is the integration time-step (first-order hold on the source
        signal between samples).
    solver : {'fft', 'state_space'}, optional
        Solver to use. ``"fft"`` uses ``BusType.impedance_formula``
        and ``BranchType.self_impedance_formula`` (frequency-domain
        impedances); ``"state_space"`` uses the lumped RLC fields on
        ``BusType`` / ``BranchType`` and integrates with
        :func:`scipy.signal.lsim`. Defaults to ``"fft"``.

    Returns
    -------
    ResultTransient
        Time-domain results at the observation points.

    Raises
    ------
    ValueError
        If no source waveform is set, if the time parameters are
        invalid, or if the network is missing the lumped RLC fields
        required by the state-space solver.
    NotImplementedError
        If ``solver`` is not one of the supported identifiers.
    """
    if solver == "fft":
        return self._solve_fft(t_end=t_end, dt=dt)
    if solver == "state_space":
        return self._solve_state_space(t_end=t_end, dt=dt)
    raise NotImplementedError(
        f"Unknown transient solver: '{solver}'. "
        "Supported: 'fft', 'state_space'."
    )

Waveforms

waveforms

Waveform Library for Transient Simulations.

This module provides a small library of factory functions that return Callable[[np.ndarray], np.ndarray] waveforms suitable for the :class:groundinsight.simulation.transient.TransientStudy solver. Each factory captures its own parameters and returns a vectorised function f(t) -> values that can be called with an array of time samples.

The included waveforms cover the typical fault-current scenarios for low- and medium-voltage grounding studies:

  • :func:step -- a Heaviside-style on/off pulse, the simplest fault model (the source comes on at t_on and off again at t_off).
  • :func:sinusoidal_with_dc_offset -- a power-frequency current with an exponentially decaying DC component, the textbook representation of a single-line-to-ground fault current including the asymmetry caused by the inductive loop.
  • :func:damped_oscillation -- a damped sinusoid useful for switching transients and ringing studies.

Custom waveforms can be defined as any user function that accepts a 1-D np.ndarray of time samples and returns an array of the same shape.

damped_oscillation

damped_oscillation(
    amplitude: float,
    frequency_hz: float,
    decay_tau: float,
    *,
    phase_rad: float = 0.0,
    t_on: float = 0.0,
    t_off: Optional[float] = None
) -> Callable[[np.ndarray], np.ndarray]

Construct a damped sinusoid, useful for switching-transient studies.

The waveform is::

x(t) = A * exp(-(t - t_on)/tau) * sin(omega*(t - t_on) + phi)

inside the on-window [t_on, t_off), and zero outside.

Parameters:

Name Type Description Default
amplitude float

Initial peak amplitude.

required
frequency_hz float

Oscillation frequency in Hz.

required
decay_tau float

Decay time constant in seconds.

required
phase_rad float

Phase offset of the sinusoid, in radians. Defaults to 0.0.

0.0
t_on float

Onset time, in seconds. Defaults to 0.0.

0.0
t_off float

Cut-off time, in seconds. None lets the oscillation decay naturally to the end of the time grid.

None

Returns:

Type Description
callable

A vectorised waveform.

Examples:

>>> w = damped_oscillation(
...     amplitude=1e3, frequency_hz=500.0, decay_tau=2e-3,
...     t_on=0.01,
... )
Source code in src/groundinsight/simulation/waveforms.py
def damped_oscillation(
    amplitude: float,
    frequency_hz: float,
    decay_tau: float,
    *,
    phase_rad: float = 0.0,
    t_on: float = 0.0,
    t_off: Optional[float] = None,
) -> Callable[[np.ndarray], np.ndarray]:
    """
    Construct a damped sinusoid, useful for switching-transient studies.

    The waveform is::

        x(t) = A * exp(-(t - t_on)/tau) * sin(omega*(t - t_on) + phi)

    inside the on-window ``[t_on, t_off)``, and zero outside.

    Parameters
    ----------
    amplitude : float
        Initial peak amplitude.
    frequency_hz : float
        Oscillation frequency in Hz.
    decay_tau : float
        Decay time constant in seconds.
    phase_rad : float, optional
        Phase offset of the sinusoid, in radians. Defaults to ``0.0``.
    t_on : float, optional
        Onset time, in seconds. Defaults to ``0.0``.
    t_off : float, optional
        Cut-off time, in seconds. ``None`` lets the oscillation decay
        naturally to the end of the time grid.

    Returns
    -------
    callable
        A vectorised waveform.

    Examples
    --------
    >>> w = damped_oscillation(
    ...     amplitude=1e3, frequency_hz=500.0, decay_tau=2e-3,
    ...     t_on=0.01,
    ... )
    """
    if not np.isfinite(frequency_hz) or frequency_hz <= 0:
        raise ValueError(
            "damped_oscillation: frequency_hz must be a finite positive "
            f"number (got {frequency_hz!r})."
        )
    if not np.isfinite(decay_tau) or decay_tau <= 0:
        raise ValueError(
            "damped_oscillation: decay_tau must be a finite positive "
            f"number (got {decay_tau!r}). decay_tau == 0 divides by zero "
            "and decay_tau < 0 produces an exponentially growing "
            "waveform — almost certainly not what was intended."
        )
    _validate_window(t_on, t_off, "damped_oscillation")
    omega = 2.0 * np.pi * frequency_hz

    def _wave(t: np.ndarray) -> np.ndarray:
        t = np.asarray(t, dtype=float)
        on = t >= t_on
        if t_off is not None:
            on = on & (t < t_off)
        tau_local = np.maximum(t - t_on, 0.0)
        envelope = amplitude * np.exp(-tau_local / decay_tau)
        signal = envelope * np.sin(omega * tau_local + phase_rad)
        return np.where(on, signal, 0.0)

    return _wave

sinusoidal_with_dc_offset

sinusoidal_with_dc_offset(
    amplitude: float,
    frequency_hz: float,
    *,
    phase_rad: float = 0.0,
    t_on: float = 0.0,
    t_off: Optional[float] = None,
    dc_amplitude: float = 0.0,
    dc_decay_tau: Optional[float] = None
) -> Callable[[np.ndarray], np.ndarray]

Construct a windowed sinusoid with optional exponentially decaying DC.

Models the typical single-line-to-ground fault current::

i(t) = A * sin(omega*(t - t_on) + phi) + I_dc * exp(-(t - t_on)/tau)

multiplied by a rectangular window between t_on and t_off. The DC component captures the asymmetric peak that arises when the fault occurs at a non-zero-crossing instant of the source voltage and the loop has finite inductance.

Parameters:

Name Type Description Default
amplitude float

Peak value of the AC component.

required
frequency_hz float

Power frequency in Hz (e.g. 50 or 60).

required
phase_rad float

Initial phase of the AC component, in radians. Defaults to 0.0.

0.0
t_on float

Fault initiation time, in seconds. Defaults to 0.0.

0.0
t_off float

Fault clearing time, in seconds. None keeps the fault on until the end of the time grid.

None
dc_amplitude float

Initial value of the DC component. Defaults to 0.0 (no DC offset).

0.0
dc_decay_tau float

Time constant L/R of the loop governing the DC decay, in seconds. None disables the decay (DC stays constant during the on-window). Required when dc_amplitude != 0 and a finite decay is desired.

None

Returns:

Type Description
callable

A vectorised waveform f(t) -> values.

Examples:

>>> w = sinusoidal_with_dc_offset(
...     amplitude=20e3, frequency_hz=50.0,
...     t_on=0.02, t_off=0.12,
...     dc_amplitude=10e3, dc_decay_tau=0.05,
... )
Source code in src/groundinsight/simulation/waveforms.py
def sinusoidal_with_dc_offset(
    amplitude: float,
    frequency_hz: float,
    *,
    phase_rad: float = 0.0,
    t_on: float = 0.0,
    t_off: Optional[float] = None,
    dc_amplitude: float = 0.0,
    dc_decay_tau: Optional[float] = None,
) -> Callable[[np.ndarray], np.ndarray]:
    """
    Construct a windowed sinusoid with optional exponentially decaying DC.

    Models the typical single-line-to-ground fault current::

        i(t) = A * sin(omega*(t - t_on) + phi) + I_dc * exp(-(t - t_on)/tau)

    multiplied by a rectangular window between ``t_on`` and ``t_off``.
    The DC component captures the asymmetric peak that arises when the
    fault occurs at a non-zero-crossing instant of the source voltage
    and the loop has finite inductance.

    Parameters
    ----------
    amplitude : float
        Peak value of the AC component.
    frequency_hz : float
        Power frequency in Hz (e.g. 50 or 60).
    phase_rad : float, optional
        Initial phase of the AC component, in radians. Defaults to ``0.0``.
    t_on : float, optional
        Fault initiation time, in seconds. Defaults to ``0.0``.
    t_off : float, optional
        Fault clearing time, in seconds. ``None`` keeps the fault on
        until the end of the time grid.
    dc_amplitude : float, optional
        Initial value of the DC component. Defaults to ``0.0`` (no DC
        offset).
    dc_decay_tau : float, optional
        Time constant ``L/R`` of the loop governing the DC decay, in
        seconds. ``None`` disables the decay (DC stays constant during
        the on-window). Required when ``dc_amplitude != 0`` and a finite
        decay is desired.

    Returns
    -------
    callable
        A vectorised waveform ``f(t) -> values``.

    Examples
    --------
    >>> w = sinusoidal_with_dc_offset(
    ...     amplitude=20e3, frequency_hz=50.0,
    ...     t_on=0.02, t_off=0.12,
    ...     dc_amplitude=10e3, dc_decay_tau=0.05,
    ... )
    """
    if not np.isfinite(frequency_hz) or frequency_hz <= 0:
        raise ValueError(
            "sinusoidal_with_dc_offset: frequency_hz must be a finite "
            f"positive number (got {frequency_hz!r}). Negative frequency "
            "is silently equivalent to flipping the phase and "
            "frequency_hz == 0 collapses to a constant offset that "
            "masks user confusion with the dc_amplitude term."
        )
    if dc_decay_tau is not None and (
        not np.isfinite(dc_decay_tau) or dc_decay_tau <= 0
    ):
        raise ValueError(
            "sinusoidal_with_dc_offset: dc_decay_tau must be a finite "
            f"positive number when set (got {dc_decay_tau!r}). Pass "
            "dc_decay_tau=None to disable the decay altogether."
        )
    _validate_window(t_on, t_off, "sinusoidal_with_dc_offset")
    omega = 2.0 * np.pi * frequency_hz

    def _wave(t: np.ndarray) -> np.ndarray:
        t = np.asarray(t, dtype=float)
        on = t >= t_on
        if t_off is not None:
            on = on & (t < t_off)

        tau_local = t - t_on
        ac = amplitude * np.sin(omega * tau_local + phase_rad)
        if dc_amplitude == 0.0:
            dc = np.zeros_like(t)
        elif dc_decay_tau is None:
            dc = np.full_like(t, dc_amplitude)
        else:
            # Exponential decay starting at t_on, only meaningful in the
            # on-window. Outside the window the rectangular factor below
            # forces it to zero anyway.
            dc = dc_amplitude * np.exp(
                -np.maximum(tau_local, 0.0) / dc_decay_tau
            )

        return np.where(on, ac + dc, 0.0)

    return _wave

step

step(
    amplitude: float,
    *,
    t_on: float = 0.0,
    t_off: Optional[float] = None
) -> Callable[[np.ndarray], np.ndarray]

Construct a rectangular pulse waveform.

The output is zero before t_on, equals amplitude between t_on and t_off and is zero again after t_off. With t_off=None the pulse extends to the end of the time grid (a classic Heaviside step).

Parameters:

Name Type Description Default
amplitude float

Plateau value of the pulse.

required
t_on float

Time at which the pulse switches on, in seconds. Defaults to 0.0.

0.0
t_off float

Time at which the pulse switches off, in seconds. None means "stays on indefinitely". Defaults to None.

None

Returns:

Type Description
callable

A vectorised waveform function f(t) -> values.

Examples:

>>> import numpy as np
>>> w = step(amplitude=100.0, t_on=0.02, t_off=0.12)
>>> y = w(np.linspace(0.0, 0.2, 5))
Source code in src/groundinsight/simulation/waveforms.py
def step(
    amplitude: float,
    *,
    t_on: float = 0.0,
    t_off: Optional[float] = None,
) -> Callable[[np.ndarray], np.ndarray]:
    """
    Construct a rectangular pulse waveform.

    The output is zero before ``t_on``, equals ``amplitude`` between
    ``t_on`` and ``t_off`` and is zero again after ``t_off``. With
    ``t_off=None`` the pulse extends to the end of the time grid
    (a classic Heaviside step).

    Parameters
    ----------
    amplitude : float
        Plateau value of the pulse.
    t_on : float, optional
        Time at which the pulse switches on, in seconds. Defaults to
        ``0.0``.
    t_off : float, optional
        Time at which the pulse switches off, in seconds. ``None`` means
        "stays on indefinitely". Defaults to ``None``.

    Returns
    -------
    callable
        A vectorised waveform function ``f(t) -> values``.

    Examples
    --------
    >>> import numpy as np
    >>> w = step(amplitude=100.0, t_on=0.02, t_off=0.12)
    >>> y = w(np.linspace(0.0, 0.2, 5))
    """
    _validate_window(t_on, t_off, "step")

    def _wave(t: np.ndarray) -> np.ndarray:
        t = np.asarray(t, dtype=float)
        on = t >= t_on
        if t_off is not None:
            on = on & (t < t_off)
        return np.where(on, amplitude, 0.0)

    return _wave

The matching matplotlib helpers plot_epr_transient and plot_branch_current_transient are documented on the Plotting page.