function [t,h]=pset7problem1
%solves the problem 4.B.1 in Beer's textbook
%input : none
%output: returns time t and height h of the tank

%initialize the parameters
params.rho = 1000; %SI units
params.mu = 0.001; %SI units
params.Dp = 0.01; %SI units
params.KL = 0.2;
params.Lp = 0.2; %SI units

h0 = 2; %initial value of h in m

options = odeset('Events', @odeevent,'RelTol',1e-4,'AbsTol',1e-14);
[t,h, TE, YE, IE] = ode23s(@odefun, [0 10000], h0,options,params);

plot(t,h); title('The height in the tank');
xlabel('Time (sec)'); ylabel('Height (m)');ylim([0.0 2.1]);

output = strcat('The tank gets empty at time:\t',num2str(TE),' s\n');
fprintf(output);

%calculates the volume in the tank for an array of given heights
V = calculateVolume(h);

figure;
plot(t,V);title('The volume of water in tank');
xlabel('Time(sec)'); ylabel('Volume (m^3)');;
return;

%*************************************************************************


function f = odefun(t,h,params)
%calculates the dhdt 

f = 0;

%solve the value of dhdt
option = optimset('Display','off');
if (h<=0)
    dhdt = 0;
elseif (h>=2.0001)
    dhdt = 0;
else
    dhdt = fsolve(@fsolvefun, 1.26e-3,option,h,params);
end

f = -1*abs(dhdt);
return

%*************************************************************************

function residue = fsolvefun(dhdt,h,params)
%unpack the parameters
rho = params.rho;

%calculate the drain velocity
vp = v_drain(h,dhdt, params);


%calculate the frictionfactor
f = ff(h, dhdt, params);

Lp = params.Lp;
KL = params.KL;
Dp = params.Dp;


residue = 1/2*dhdt^2 + 9.8 * h - vp^2/2 + 9.8*Lp - KL * (1/2*vp^2) - f * (Lp/Dp) * (1/2*vp^2);


return;


%*************************************************************************

function f = ff(h, dhdt, params)
%caclulates the friction factor

%calculate the drain velocity
vp = v_drain(h,dhdt, params);


%calculate the reynolds number
Re = Reynolds(vp, params);
option = optimset('Display','off');
if (Re < 2100)
    f = 64/Re;
elseif (Re >4000)
    f = fsolve(@solvef, 3.1e-2,option,Re,params);
    %fprintf('The Reynolds number are greater than 4000\n');
else
    f = 64/Re; %this will print a error message
    fprintf('The Reynolds number is between 2100 and 4000\n');
end

return;
 


%*************************************************************************

function vp = v_drain(h, dhdt, params)
%returns the drain velocity vp

%unpack the parameters
rho = params.rho;
mu = params.mu;
Dp = params.Dp;

d = (diameter(h));

vp = (d^2/Dp^2)*abs(dhdt);
return;

%*************************************************************************

function Re = Reynolds(vp, params)
%calculates the reynolds number

%unpack the parameters
rho = params.rho;
mu = params.mu;
Dp = params.Dp;

Re = rho * vp * Dp/mu;
return;


%*************************************************************************

function V = calculateVolume(h_given)
%for a given array of heights it returns an array of volumes.

for i = 1:length(h_given)
    h(length(h_given)-i+1) = h_given(i);
end

[h_out,V_calc] = ode45(@diffVdh,h,0,[]);

V = zeros(length(V_calc),1);
for i = 1:length(V_calc)
    V(i) = V_calc(length(V_calc)+1-i);
end

return;

%*************************************************************************

function f = diffVdh(h,V);

d = diameter(h);

f = pi*d^2/4;
return;

%*************************************************************************

function d = diameter(h)
%calculate the diameter at a given h

if (h < -0.000001)
    d = 35;
    fprintf('The height cannot be less than 0 \n height is:\t',num2str(h),' m\');

elseif(h <= 0.5)
    d = 35;
elseif(h > 0.5 && h<=0.75)
    d = 35 +  5 * (h-0.50)/0.25;
elseif(h > 0.75 && h<=1)
    d = 40 + 15 * (h-0.75)/0.25;
elseif(h > 1 && h<=1.25)
    d = 55 + 30 * (h-1.00)/0.25;
elseif(h > 1.25 && h<=1.50)
    d = 85 + 15 *(h-1.25)/0.25;
elseif(h>1.50 && h<=2.00001)
    d = 100;
else
    d = 100;%this will print an error message
    fprintf('The height cannot be greater than 2 m \n');
    fprintf('The height is:\t',num2str(h),' m\');
end

d = d/100; %convert to SI units

return

%*************************************************************************

function r = solvef(f,Re,params)

Dp = params.Dp;

r = 1/f^0.5 + 2*log10(0.045e-3/Dp/3.7 + 2.51/Re/f^0.5);
return;

%************************************************************************

function [value,isterminal,direction] = odeevent(t,x,params)
%is the events function used in ode23s which stops the integration as soon
%as h becomes less than 0.

isterminal(1) = 1;
direction(1) = 0;
value(1) = x;

return

%************************************************************************