Frontal polymerization Angleply 8layer

def _kamal_sourour_rate_vec(T: np.ndarray, alpha: np.ndarray) -> np.ndarray:
    """Vectorised Kamal-Sourour autocatalytic cure rate.

    Clamps inputs for numerical safety before evaluating the Arrhenius
    terms.  Returns dα/dt for each node.
    """
    T_safe = np.clip(T, _T_FLOOR, _T_CEIL)
    alpha_safe = np.clip(alpha, 0.0, 1.0)

    inv_RT = 1.0 / (_R_GAS * T_safe)
    arg1 = np.clip(_E1 * inv_RT, 0.0, _EXP_ARG_MAX)
    arg2 = np.clip(_E2 * inv_RT, 0.0, _EXP_ARG_MAX)

    k1 = _A1 * np.exp(-arg1)
    k2 = _A2 * np.exp(-arg2)

    return (k1 + k2 * np.power(alpha_safe, _M)) * np.power(1.0 - alpha_safe, _N_ORD)

def _fp_angleply_8layer_rhs(t, y):
    T = np.clip(y[:_AP_N_SPATIAL], _T_FLOOR, _T_CEIL)
    alpha = np.clip(y[_AP_N_SPATIAL:2 * _AP_N_SPATIAL], 0.0, 1.0)
    P = y[2 * _AP_N_SPATIAL:]

    dadt = _kamal_sourour_rate_vec(T, alpha)

    dT = np.empty(_AP_N_SPATIAL)
    dP = np.empty(_AP_N_SPATIAL)

    for layer in range(_AP_LAYERS):
        start = layer * _AP_NODES_PER_LAYER
        end = start + _AP_NODES_PER_LAYER

        for i in range(start, end):
            local = i - start

            if i == 0:
                T_left = _AP_T_BOTTOM
            elif local == 0:
                # Inter-ply interface (uniform k on both sides for ±45°)
                R_half = _AP_DX / (2.0 * _AP_K_EFF)
                R_left = R_half + _R_CONTACT
                R_right = R_half
                T_left = (T[i - 1] / R_left + T[i] / R_right) / (1.0 / R_left + 1.0 / R_right)
            else:
                T_left = T[i - 1]

            if i == _AP_N_SPATIAL - 1:
                # Top node: convective BC
                T_right = T[i] + (_H_CONV * _AP_DX / _AP_K_EFF) * (_T_AMBIENT - T[i])
            elif local == _AP_NODES_PER_LAYER - 1 and layer < _AP_LAYERS - 1:
                R_half = _AP_DX / (2.0 * _AP_K_EFF)
                R_left = R_half
                R_right = R_half + _R_CONTACT
                T_right = (T[i] / R_left + T[i + 1] / R_right) / (1.0 / R_left + 1.0 / R_right)
            else:
                T_right = T[i + 1]

            lap = (T_left - 2.0 * T[i] + T_right) * _AP_INV_DX2
            dT[i] = _AP_DIFF * lap + _SRC_COEFF * dadt[i]

    dT[0] = 0.0

    dP[:] = (
        (_RHO_RESIN * _V_GAS_SPECIFIC * dadt * _R_GAS_IDEAL * T) / _V_PORE
        - P * _PERM_LOSS
    )

    dy = np.empty(_AP_DIM)
    dy[:_AP_N_SPATIAL] = dT
    dy[_AP_N_SPATIAL:2 * _AP_N_SPATIAL] = dadt
    dy[2 * _AP_N_SPATIAL:] = dP
    return dy