def warmtestap(c, z, tout, qm_req, qh_req, qwp, qnood,
               gains_m, gains_h, extra_H, hwoon, dt):
    """RK4 integreert toestanden EN dezelfde fluxen voor de balans."""
    capaciteiten = np.array([c["Ca_J_K"], c["Cb_J_K"],
        c["Cwoon_J_K"], c["Cwp_J_K"]])
    def f(t):
        a, b, h, w = t
        qm = min(qm_req, c["Gconvector_W_K"]*max(0, w-a))
        qh = min(qh_req, c["Gvloer_W_K"]*max(0, w-h))
        qbw = c["Hwp_naar_woon_W_K"]*(w-h)
        g = c["Gmodule_W_K"]*(a-b)
        lm = (c["Hmodule_W_K"]+extra_H)*(a-tout)
        lh = hwoon*(h-tout)
        vermogen = np.array([qm+gains_m-lm-g, g,
            qh+qbw+qnood+gains_h-lh, qwp-qm-qh-qbw])
        flux = np.array([qm, qh, lm, lh, qbw])
        return vermogen/capaciteiten, flux
    voor = z.copy()
    som = np.zeros(5)
    for _ in range(round(c["besturing_s"]/dt)):
        k1, f1 = f(z)
        k2, f2 = f(z+dt*k1/2)
        k3, f3 = f(z+dt*k2/2)
        k4, f4 = f(z+dt*k3)
        z = z+dt*(k1+2*k2+2*k3+k4)/6
        som += dt*(f1+2*f2+2*f3+f4)/6
    delta = float(capaciteiten@(z-voor))
    netto = (qwp+qnood+gains_m+gains_h)*c["besturing_s"]-som[2]-som[3]
    assert abs(delta-netto) < 1e-5
    return z, som, delta-netto