Multi-Node Thermal Network (13-state)

def _mission_profile(t: float) -> tuple[float, float]:
    """Return (motor_current_per_motor_A, airspeed_ms) for each phase."""
    if t < 30.0:
        return 4.0, 0.5        # spin-up / hover
    elif t < 120.0:
        return 3.0, 8.0        # cruise
    elif t < 180.0:
        return 5.5, 2.0        # aggressive hover / filming
    elif t < 300.0:
        return 3.0, 8.0        # cruise return
    elif t < 360.0:
        return 3.5, 1.0        # hover descent
    else:
        return 1.0, 0.0        # idle

def rhs_thermal_envelope(t, y):
    """13-state multi-node thermal network for a quadrotor."""
    dy = np.zeros(13)

    T_batt = y[0]
    T_motors = y[1:5]
    T_escs = y[5:9]
    T_frame = y[9]
    T_payload = y[10]

    I_motor, airspeed = _mission_profile(t)
    I_pack = I_motor * 4

    # Convective enhancement from airspeed (forced convection scaling)
    conv_factor = 1.0 + 2.0 * math.sqrt(max(0, airspeed))

    # Battery heat generation (I^2 R)
    P_batt = I_pack**2 * _BATT_R_INTERNAL

    # Motor heat generation (copper losses, per motor)
    P_motors = np.array([I_motor**2 * _MOTOR_R_PHASE] * 4)
    # Add iron losses (proportional to speed, approximated)
    P_motors += 0.3 * I_motor

    # ESC heat generation (conduction + switching, per ESC)
    P_escs = np.array([I_motor**2 * _ESC_R_DS_ON * 2 + 0.1 * I_motor] * 4)

    # Battery thermal dynamics
    Q_batt_to_frame = (T_batt - T_frame) / _R_BATT_FRAME
    Q_batt_conv = _H_BATT_BASE * conv_factor * (T_batt - _T_AMB)
    dy[0] = (P_batt - Q_batt_to_frame - Q_batt_conv) / _C_BATT

    # Motor thermal dynamics (4 motors)
    for i in range(4):
        Q_motor_frame = (T_motors[i] - T_frame) / _R_MOTOR_FRAME
        Q_motor_conv = _H_MOTOR_BASE * conv_factor * (T_motors[i] - _T_AMB)
        dy[1 + i] = (P_motors[i] - Q_motor_frame - Q_motor_conv) / _C_MOTOR

    # ESC thermal dynamics (4 ESCs)
    for i in range(4):
        Q_esc_frame = (T_escs[i] - T_frame) / _R_ESC_FRAME
        Q_esc_conv = _H_ESC_BASE * conv_factor * (T_escs[i] - _T_AMB)
        dy[5 + i] = (P_escs[i] - Q_esc_frame - Q_esc_conv) / _C_ESC

    # Frame thermal dynamics (receives heat from all components)
    Q_in_frame = (
        (T_batt - T_frame) / _R_BATT_FRAME +
        sum((T_motors[i] - T_frame) / _R_MOTOR_FRAME for i in range(4)) +
        sum((T_escs[i] - T_frame) / _R_ESC_FRAME for i in range(4))
    )
    Q_frame_air = (T_frame - _T_AMB) / _R_FRAME_AIR
    Q_frame_payload = (T_frame - T_payload) / _R_PAYLOAD_FRAME
    dy[9] = (Q_in_frame - Q_frame_air - Q_frame_payload) / _C_FRAME

    # Payload thermal dynamics
    dy[10] = (T_frame - T_payload) / (_R_PAYLOAD_FRAME * _C_PAYLOAD)

    # Cumulative heat tracking
    dy[11] = P_batt
    dy[12] = P_batt + float(np.sum(P_motors)) + float(np.sum(P_escs))

    return dy