clc
close all

Q=0.35      %m3
dx=50;
n_max=200;      %antal dx
I0=0.2*10^-3;
b=4;        %Bundbredde
a=1;        %anlæg
alfa=1.2;
M=30;
%-----------Naturlig dybde----------
Q=0.35;
i_gaet=0.5; %startværdi
%Iteration af manningformlen vha intvalhalveringsmetoden-----------------------------
y_ovre=100;      %m
y_nedre=0.001;  %m
while abs(y_ovre-y_nedre)>0.0001
y=(y_ovre+y_nedre)/2;
A=b*y+y^2*a;
P=b+2*y*sqrt(1+a^2);
R=A/P;
i_gaet=(Q/A/(M*R^(2/3)))^2;
if i_gaet<I0
    y_ovre=y;
else
    y_nedre=y;
end
end
naturlig_dybde=y

%------Sænknings/stuvningskurve--------------

for y_start=[0.01 0.6 0.9 1.5]; % forskellige nedre randbetingelser y=0 svarer til kote -0.5
stop=1e-4;
g=9.81;
fejl=10;
Res=zeros(n_max,2);
yn=y_start;
y_ny=0;
for i=1:n_max;
y_op=1.3*yn;
y_ned=0.8*yn;
    while y_op-y_ned>stop
        y_ny=0.5*y_op+0.5*y_ned;
        A_ny=y_ny*(b+a*y_ny);
        P_ny=b+2*y_ny*sqrt(1+a^2);
        R_ny=A_ny/P_ny;
        An=yn*(b+a*yn);
        Pn=b+2*yn*sqrt(1+a^2);
        Rn=An/Pn;
        delta_H=(Q/(0.5*A_ny+0.5*An)/(M*(0.5*R_ny+0.5*Rn)^(2/3)))^2*dx;
        fejl=y_ny-(yn-I0*dx+alfa*Q^2/(2*g)*(1/An^2-1/A_ny^2)+delta_H);
        if fejl>=0
        y_op=y_ny;
        elseif fejl<0
        y_ned=y_ny;
        end
    end
yn=y_ny;
    %Res(i,1)=i*dx;
    Froude=(Q/A_ny)/sqrt(g*A_ny/(b+2*y_ny*a));
    Res(i,1)=Froude;
    Res(i,2)=y_ny;
end
hold on
x_akse=dx:dx:n_max*dx;
y_bund(:,1)=x_akse.*I0-0.5;
plot(x_akse,y_bund,'+black')
Xlabel('afstand fra fjord i meter')
Ylabel('vandstand kote')
Res(:,2)=Res(:,2)+y_bund;
plot(x_akse,x_akse*I0-0.5+naturlig_dybde,'r--') %plotter naturlig dybde
plot(x_akse,Res(:,2),'LineWidth',2);                          %plotter sænkningskurve

end

froude_max=max(Res(:,1))
froude_min=min(Res(:,1))