% The leaky bucket problem
%
% Statement: Water is pouring into cylindrical bucket with a small hole in the 
% bottom:
%       1. Calculate the height of water in the bucket over time
%       2. Find the time when the bucket will be full.
%
% T.O and A.P - 25 Sep 2015

%% ---------------- Define problem parameters and constants --------------

Qin  = 7;            % Inflow                        [ml/s]
r    = 2.5;          % Radius bucket                 [cm]
A    = pi*r^2;       % Bucket area                   [cm^2]
H    = 30;           % Bucket height                 [cm]
rh   = 0.1;          % Bottom hole radius            [cm]
Ah   = pi*rh^2;      % Area hole                     [cm^2]
g    = 981;          % Acceleration due to gravity   [cm/s^2]
rho  = 1;            % Density of water              [g/cm^3]
mu   = 0.01;         % Viscosity of water            [g/(cm s)]
dt   = 0.01;         % Time step                     [s]
total_time = 250;    % Total simulation time         [s]
step = 0:250/dt;     % Iteration steps               [step]
t    = step*dt';     % Time of simulation            [s]

% Look up Reynolds number and coefficient of discharge.
Re_table = [1    10   100  1000 10000 100000]; % Reynolds number
Re_table = log10(Re_table); % Take log10 to allow linear interpolation.
CD_table = [0.04 0.28 0.73 0.74 0.64  0.61  ];% Coefficient of discharge


%% ----------------------- Pre-allocate variables -------------------------

L = length(step);
h = nan(L,1);
Qout = nan(L,1);
Re = nan(L,1);
CD = nan(L,1);

%% ----------------------- Pre-allocate variables -------------------------

% Time t=0. No flow into the bucket yet.
h(1) = 0;      % Bucket is empty at time t=0.
Qout(1) = 0;  % No flow since the bucket is empty.

% Time t=1
h(2) = h(1) + dt*Qin/A;             % Assume no leaking in the first t step
Qout(2) = 0;

% Time t=2 to end. Notice that t=2 is in t(3).
for i = 3:L
    Re(i)       = (2*rh*rho)/mu*sqrt(g*h(i-1));  
    CD(i)       = interp1(Re_table,CD_table,log10(Re(i)));
    Qout(i)     = CD(i)*Ah*sqrt(2*g*h(i-1));       % [cm^3/s]
    h_out       = Qout(i)/A;
    h_in        = Qin/A;
    h_increment = dt*(h_in - h_out);
    h(i)        = min(h(i-1) + h_increment, H);
end

% Detect maximum height
row_full = find(h==H,1); % Check if bucket is full. 
time_full = t(row_full); % Find time at which the bucket was full

if ~isempty(row_full)
    display(['The bucket is filled after ',num2str(time_full),' seconds']);
else
    display('You need more time or higher Qin to fill this bucket');
end

%% Plotting
font_axis = 13;
figure
subplot(1,2,1)
plot(t,h)
ylim([0 max(h)+max(h)*0.1])
xlabel('Time (Seconds)','FontSize',font_axis)
ylabel('Height of water in the bucket (cm)','FontSize',font_axis)
set(gca,'FontSize',font_axis)

subplot(1,2,2)
scatter(log10(Re), CD,'o'); hold on
xlabel('Reynolds Number','FontSize',font_axis)
ylabel('Coefficient of discharge','FontSize',font_axis)

Re_interp = log10(0:10:1e5);
CD_interp = interp1(Re_table,CD_table,Re_interp);
plot(Re_interp, CD_interp,'--r')
xlabel('log_{10} Re','FontSize',font_axis)
ylabel('Coefficient of discharge','FontSize',font_axis)

legend('Re in simulation',...
       'Complete CD(Re) interpolation',...
       'Location', 'SE');
legend boxoff
set(gca,'FontSize',font_axis)
box on
set(gcf,'Position',[50 50 1000 500])