IEEE 9-bus fault at bus 7 (18D, pre-segmented)

def _ieee9_fault_bus7_rhs(t: float, y: np.ndarray) -> np.ndarray:
    """IEEE 9-bus 3-machine 6th-order subtransient model.

    Y_red switches at t=1.0 s (fault at bus 7) and t=1.15 s (line 5-7
    tripped).  State layout: 6 states per machine packed consecutively.
    """
    if t < _I9_FAULT_ON:
        Y_red = _Y9_PRE
    elif t < _I9_FAULT_OFF:
        Y_red = _Y9_FLT
    else:
        Y_red = _Y9_PST

    delta = np.array([y[0], y[6], y[12]])
    omega = np.array([y[1], y[7], y[13]])
    Eq_p = np.array([y[2], y[8], y[14]])
    Ed_p = np.array([y[3], y[9], y[15]])
    psi_ppd = np.array([y[4], y[10], y[16]])
    psi_ppq = np.array([y[5], y[11], y[17]])

    Vd_pp = -psi_ppq
    Vq_pp = psi_ppd

    sin_d = np.sin(delta)
    cos_d = np.cos(delta)
    E_D = Vd_pp * sin_d + Vq_pp * cos_d
    E_Q = -Vd_pp * cos_d + Vq_pp * sin_d

    I_complex = Y_red @ (E_D + 1j * E_Q)
    I_D = I_complex.real
    I_Q = I_complex.imag

    Id = I_D * sin_d - I_Q * cos_d
    Iq = I_D * cos_d + I_Q * sin_d
    Pe = Vd_pp * Id + Vq_pp * Iq

    dy = np.zeros(18)
    for i in range(3):
        j = i * 6
        dy[j] = _I9_ws * (omega[i] - 1.0)
        dy[j + 1] = (
            (_I9_Pm[i] - Pe[i] - _I9_D[i] * (omega[i] - 1.0))
            / (2.0 * _I9_H[i])
        )
        dy[j + 2] = (
            (_I9_Efd[i] - Eq_p[i] - (_I9_Xd[i] - _I9_Xdp[i]) * Id[i])
            / _I9_Td0p[i]
        )
        dy[j + 3] = -(Ed_p[i] + (_I9_Xq[i] - _I9_Xqp[i]) * Iq[i]) / _I9_Tq0p[i]
        dy[j + 4] = (
            (Eq_p[i] - psi_ppd[i] - (_I9_Xdp[i] - _I9_Xppd[i]) * Id[i])
            / _I9_Td0pp[i]
        )
        # d/q-symmetric with the psi''d equation above; the previous leading
        # minus sign made the psi''q self-term positive feedback (+1/T''q0),
        # giving 3 unstable real eigenvalues (+21.7/+12.7/+4.6 /s) at the
        # pre-fault equilibrium (model audit 2026-07-17, jobs 10592/10593).
        dy[j + 5] = (
            (Ed_p[i] - psi_ppq[i] + (_I9_Xqp[i] - _I9_Xppq[i]) * Iq[i])
            / _I9_Tq0pp[i]
        )
    return dy