Frontal polymerization Crossply 3layer

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_crossply_3layer_rhs(t, y):
    T = np.clip(y[:_CP1_N_SPATIAL], _T_FLOOR, _T_CEIL)
    alpha = np.clip(y[_CP1_N_SPATIAL:2 * _CP1_N_SPATIAL], 0.0, 1.0)
    P = y[2 * _CP1_N_SPATIAL:]

    dadt = _kamal_sourour_rate_vec(T, alpha)

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

    for layer in range(_CP1_LAYERS):
        start = layer * _CP1_NODES_PER_LAYER
        end = start + _CP1_NODES_PER_LAYER
        diff_layer = _CP1_DIFF[layer]

        for i in range(start, end):
            local = i - start  # node index within layer

            if i == 0:
                T_left = _CP1_T_BOTTOM
            elif local == 0:
                # First node of a non-bottom layer: inter-ply interface
                prev_layer = layer - 1
                k_left = _CP1_K_THRU[prev_layer]
                k_right = _CP1_K_THRU[layer]
                R_left = _CP1_DX / (2.0 * k_left) + _R_CONTACT
                R_right = _CP1_DX / (2.0 * k_right)
                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 == _CP1_N_SPATIAL - 1:
                # Top node: convective BC  ->  k * dT/dx = h * (T_amb - T)
                # One-sided ghost: T_ghost = T[i] + (h*dx/k)*(T_amb - T[i])
                k_top = _CP1_K_THRU[layer]
                T_right = T[i] + (_H_CONV * _CP1_DX / k_top) * (_T_AMBIENT - T[i])
            elif local == _CP1_NODES_PER_LAYER - 1 and layer < _CP1_LAYERS - 1:
                # Last node of a non-top layer: inter-ply interface
                k_left = _CP1_K_THRU[layer]
                k_right = _CP1_K_THRU[layer + 1]
                R_left = _CP1_DX / (2.0 * k_left)
                R_right = _CP1_DX / (2.0 * k_right) + _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) * _CP1_INV_DX2
            dT[i] = diff_layer * lap + _SRC_COEFF * dadt[i]

    # Bottom node held at constant temperature
    dT[0] = 0.0

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

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