clear;clc;

% Parameters and constants
Qin = 7;
nbuckets = 3;
r = 2.5;
A = pi*r^2;
H = 30;
rh = 0.1;
Ah = pi*rh^2;
g = 981;
rho = 1;
mu = 0.01;
dt = 50;
total_time = 250;
t_vec = [0:round(total_time/dt)]';
time = t_vec*dt;

% Pre-allocate
h = nan(length(t_vec),nbuckets);
Qout = nan(length(t_vec),nbuckets);

% Reynolds number table
Re_table = log10([1 10 100 1000 10000 100000]);
CD_table = [0.04 0.28 0.73 0.74 0.64 0.61];

% Initial conditions
h(1,:) = 0; % First bucket empty at t=0
Qout(1,:) = 0;

for t = 2:length(t_vec)
    h_in = Qin/A;   
    Re(t) = (2*rh*rho)/mu *sqrt(2*g*h(t-1));
    CD(t) = max(interp1(Re_table,CD_table,log10(Re(t))),0);
    Qout(t) = CD(t)*Ah*sqrt(2*g*h(t-1));
    h_out(t) = Qout(t)/A;
   
    h(t) = min(h(t-1)+dt*(h_in-h_out(t)),H);
end

% Plots
plot(time,h,'--b','LineWidth',2)
xlabel('Time (s)')
ylabel('Height (cm)')
ylim([0 H+5])
legend('Explicit method')
legend boxoff
hold all
