%% Calculating water height in a channel
% Bernoulli equation 
%
% $\frac{u^2}{2g} + \frac{p}{\rho g} + (z_b + h - d)  = H$
%
% If the horizontal flow is steady and there are vertical components then
% 
% $p = \rho g d \rightarrow \frac{p}{\rho g}= d$ 
%
% if $d$ is $h$ then this simplifies to 
% 
% $\frac{u^2}{2g} + h + z_b  = H$
%
% we can equate energy $E + z_b = H$ and if we assume that $z_b=0$ then we
% get relationship between the Energy and the water volume
%
% $$E=\frac{U^2}{2g}+h$$
%
% We assume that the channel is rectangular 
% and that we know the energy and the amount of water that flows out
% 
% $$Q = Uwh$$  
%
% with width w=width, and h=height
%
% $$q_w = Q/w = Uh$$  
%
% were $q_ww$ is the outflow, using $U = qw/h$ we can estimate the the height 
% if we know the energy E and the flow rate $q_w$
%
% $$E = \frac{q_w^2}{2 g h^2} + h$$
% or
% $$ h^3 - Eh^2 + \frac{q_w^2}{2g} = 0$$
%
%% Simplified Bernoulli fluid equation
% this equation is in the form E=f(x)
% for the root finding further down in the script we will need to change it
% to g(x)=0., for example using 0=f(x)-E
bfluid = @(qw,h,g) qw.^2./(2.*g.* h .* h) + h
%% Plot of  energy x depth (h)
depth = -0.5:0.01:1.5;
g = 9.81
qw = 0.5
energy = bfluid(qw,depth,g);
hold off
h1 = plot(0.75 * ones(1,length(depth)),depth,'--');
hold on
h2 = plot(energy, depth);
set([h1 h2],'LineWidth',2)
title('Energy versus depths of channel');
xlabel('Energy');
ylabel('Depth');
axis([0,1.5,-0.5,1.5])
%% Finding the height for Energy=0.75 m and for $q_w = 0.5 m^2/s$
energy = 0.75;
results = NaN*ones(10,1);
startingpoints=linspace(-0.5,1.5,10);
for i=1:10
results(i) = fzero(@(x) bfluid(qw,x,g)-energy,startingpoints(i));
end
%% All solutions for h :
results = uniquetol(results) %results(diff(results)>1e-12)
% and potting them into the graph
plot([energy energy energy],results,'bo')
%% A step in the water, how does the height change, before the step
% we use the E=0.8 before the step with hight difference of 0.1
energy1 = 0.8;
results = NaN*ones(10,1);
startingpoints=linspace(-0.5,1.5,10);
for i=1:10
results(i) = fzero(@(x) bfluid(qw,x,g)-energy1,startingpoints(i));
end
%% All solutions for h (before step):
results2 = uniquetol(results)
%% A step in the water, how does the height change, during the step $\Delta z$
deltaz=0.1
energy2 = energy1-deltaz;
results = NaN*ones(10,1);
startingpoints=linspace(-0.5,1.5,10);
for i=1:10
results(i) = fzero(@(x) bfluid(qw,x,g)-energy2,startingpoints(i));
end
%% All solutions for h (after step):
results2 = uniquetol(results) %results(diff(results)>1e-12)
plot([energy2 energy2 energy2],results2,'bx')
