function [s, L, E] = step_rk4(s0, dt, P, n, Zr, sstar, sw, ETmax, ETsw, Ks, b)
% STEP_RK4 Aggiornamento di un timestep dello stato del suolo, RK4, a bilancio
% di massa garantito: P - E - L = (s_new - s0)*n*Zr  ESATTAMENTE, per ogni step.
%
% Input
% s0     [-]    = soil moisture iniziale (relative saturation)
% dt     [h]    = timestep
% P      [mm]   = pioggia caduta durante il timestep (infiltrazione totale, no runoff)
% n      [-]    = porosita'
% Zr     [mm]   = profondita' radici
% sstar  [-]    = soglia per ET massima
% sw     [-]    = wilting point
% ETmax  [mm/h] = ET massima (a s>=sstar)
% ETsw   [mm/h] = ET a wilting point (a s<=sw)
% Ks     [mm/h] = conducibilita' idraulica satura
% b      [-]    = esponente di Clapp-Hornberger
%
% Output
% s      [-]    = soil moisture aggiornato
% L      [mm]   = leaching nel timestep (gia' corretto per eventuale clamp)
% E      [mm]   = evapotraspirazione nel timestep (gia' corretta per eventuale clamp)

nZr = n * Zr;
I_rate = P / dt;   % [mm/h] tasso di infiltrazione costante nel passo

%% -------- STEP 1 --------
E1 = compute_ET(s0, sw, sstar, ETsw, ETmax);
L1 = compute_L(s0, sw, Ks, b);

%% -------- STEP 2 --------
s2 = s0 + 0.5*(I_rate - E1 - L1)*dt/nZr;
E2 = compute_ET(s2, sw, sstar, ETsw, ETmax);
L2 = compute_L(s2, sw, Ks, b);

%% -------- STEP 3 --------
s3 = s0 + 0.5*(I_rate - E2 - L2)*dt/nZr;
E3 = compute_ET(s3, sw, sstar, ETsw, ETmax);
L3 = compute_L(s3, sw, Ks, b);

%% -------- STEP 4 --------
s4 = s0 + (I_rate - E3 - L3)*dt/nZr;
E4 = compute_ET(s4, sw, sstar, ETsw, ETmax);
L4 = compute_L(s4, sw, Ks, b);

%% -------- RK4: medie pesate dei tassi, poi integrazione su dt --------
Ebar = (E1 + 2*E2 + 2*E3 + E4)/6;   % [mm/h]
Lbar = (L1 + 2*L2 + 2*L3 + L4)/6;   % [mm/h]

E = Ebar * dt;   % [mm]
L = Lbar * dt;   % [mm]

%% -------- Aggiornamento soil moisture (bilancio "grezzo") --------
s_raw = s0 + (P - E - L) / nZr;

%% -------- Clamp CONSERVATIVO: ogni correzione su s viene rimessa nei flussi --------
if s_raw > 1
    % eccesso non immagazzinabile: scaricato come leaching aggiuntivo
    excess_mm = (s_raw - 1) * nZr;   % [mm], > 0
    L = L + excess_mm;
    s_raw = 1;
elseif s_raw < 0
    % deficit: l'ET "potenziale" calcolata non era sostenibile dall'acqua disponibile
    deficit_mm = (0 - s_raw) * nZr;  % [mm], > 0
    E = E - deficit_mm;
    s_raw = 0;
end

s = s_raw;

%% -------- Check di coerenza locale (rimuovibile, utile in debug) --------
% closure = P - E - L - (s - s0)*nZr;  % deve essere ~0 a meno di eps numerico
% assert(abs(closure) < 1e-9, 'step_rk4: bilancio locale non chiuso (%.3e mm)', closure);

end

%% ======================================================
% Evapotranspiration function
% ======================================================
function E = compute_ET(s, sw, sstar, ETsw, ETmax)
if s <= sw
    E = ETsw;
elseif s >= sstar
    E = ETmax;
else
    rel = (s - sw) / (sstar - sw);
    E = ETsw + rel * (ETmax - ETsw);
end
end

%% ======================================================
% Leaching function (Clapp-Hornberger: K = Ks * Se^(2b+3))
% ======================================================
function L = compute_L(s, sw, Ks, b)
if s <= sw
    L = 0;
else
    Se = (s - sw) / (1 - sw);
    Se = max(0, min(1, Se));
    L = Ks * Se^(2*b + 3);
end
end
