clear;clc;

% Parameters and constants
Qin = 7;
nbuckets = 1;
r = 2.5;
A = pi*r^2;
H = 30;
rh = 0.1;
Ah = pi*rh^2;
g = 981;
rho = 1;
mu = 0.01;
dt = 10;
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];

% Equivalent input height (this is constant and could be removed
% outside the loop)
h_in = Qin/A;
    
% Initial conditions
h(1) = 0; % First bucket empty at t=0
Qout(1) = 0;

for t = 2:length(t_vec) 
    
    % Instead of using h(t-1) to estimate h(t), 
    h2_fun = @(h2) h2 - h(t-1) - dt*(h_in -...
                                    (interp1(Re_table,...
                                    CD_table,...
                                    max(log10((2*rh*rho)/mu *sqrt(2*g*h2)),1))...
                                    *Ah...
                                    *sqrt(2*g*h2))/A);
    x0 = h(t-1)+h_in*dt;
    h(t) = fzero(h2_fun,x0);  % Predict (used to estimate outflow
    h(t) = min(h(t),H);
    
end

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

