%% initiatilisation
close all 
clear all

% paramètres généraux
g=9.81
q = 10 ; % débit par unité de largeur
Ch = 50 ; % Chézy

% paramètres bief 1
i1 = 0.05 ; % pente
h0 = 2 ; %hauteur initiale en x = 0
x0=0;
xA=10;

% paramètres bief 2
i2 = 0.002; % pente
p = 0.5 ; % pelle en B
xB=1000+xA;

%% calcul des conditions hydrauliques dans le bief 1
hn1 = (q/Ch/sqrt(i1))^(2/3)
Frn1 = q/hn1^1.5/sqrt(g) % Froude pour le régime asymptotique quand h= h_n
hc1 = (q^2/9.81)^(1/3)
Fr1 = q/2^1.5/sqrt(g) % Froude initial quand h= h_0

% calcul des conditions hydrauliques dans le bief 2
hn2 = (q/Ch/sqrt(i2))^(2/3)
Frn2 = q/hn2^1.5/sqrt(g) % Froude pour le régime asymptotique quand h= h_n
hc2 = (q^2/9.81)^(1/3)

%%solution de la courbe de remous (équation de Bessin) 
% sur l'intervalle [0, xA]
% 
[x1,y1] = ode45( @(x, y) besse (x, y, i1, hn1, hc1) , [x0 xA], h0); % Résolution de l'équation de Bresse

hA=y1(end) % hauteur d'eau en A déduite de l'équation de Besse précédente

%% solution de la courbe de remous (équation de Bessin) 
% sur l'intervalle [xA, 100]
[x2,y2] = ode45( @(x, y) besse (x, y, i2, hn2, hc2) , [xA 100], hA); % Résolution de l'équation de Bresse



subplot(121)
plot(x1,y1,'b')
title('tracé de la courbe de remous bief 1')
subplot(122)
plot(x2,y2,'b')
title('tracé de la courbe de remous bief 2 : branche supercritique')


% calcul de la condition à la limite à l'aval (point B)

% Atenntion aux options : un compromis entre précision et temps de calcul
% options = odeset('RelTol',1e-6,'AbsTol',1e-6);
% [x2,y2] = ode45( @(x, y) besse (x, y, i2, hn2, hc2) , [xA 100], hA, options); % Résolution
% 
% figure
% plot(x1,y1)
% hold on
% plot(x2,y2, 'r')

%% Attention : l'équation admet une singularité en h = h_c
% traiter cela avec l'option Events
options = odeset('RelTol',1e-4,'AbsTol',1e-4,'Events',@critique);
[x2,y2] = ode45( @(x, y) besse (x, y, i2, hn2, hc2) , [xA 100], hA, options); % Résolution

figure
plot(x1,y1, 'b')
hold on
plot(x2,y2, 'b')
title('courbe de remous $h(x)$ dans le bief 2 (branche supercritique)','Interpreter','latex')
xlabel('$x$ [m]','Interpreter','latex');
ylabel('$h$ [m]','Interpreter','latex');

%% calcul de la solution subcritique dans le bief 2
% calcul de la hauteur en B
func = @(h) (q/sqrt(g))^(2/3)*3/2+p-h-q^2/2/h^2/g;
x0 = 3; % starting point. Attention si vous prenez une valeur trop basse vous aurez la solution en régime supercritique
hB = fzero(func,x0);
fprintf("La hauteur d'eau au-dessus du seuil est h = %1.2f m.\n",hB)

% résolution de l'équation (dans le sens inverse) de l'équation de Besse
options = odeset('RelTol',1e-5,'AbsTol',1e-5,'Events',@critique);
[x3,y3] = ode45( @(x, y) besse (x, y, i2, hn2, hc2) , [xB xA], hB, options); % on intègre depuis xB donc on met [xB yB]

figure
plot(x3,y3, 'r')
title('courbe de remous $h(x)$ dans le bief 2 (branche subrcritique)','Interpreter','latex')
xlabel('$x$ [m]','Interpreter','latex');
ylabel('$h$ [m]','Interpreter','latex');

%% tracé de la courbe conjuguée et calcul du point d'intersection
conjugaison = @(h) h*(sqrt(8 *q^2/g/h^3+1)-1)*0.5;

y2conjugue = arrayfun(conjugaison,y2);
figure
plot(x2,y2, 'b')
hold on
plot(x2,y2conjugue, 'b--')
plot(x3,y3, 'r')
xlabel('$x$ [m]','Interpreter','latex');
ylabel('$h$ [m]','Interpreter','latex');
title('courbe de remous $h(x)$','Interpreter','latex')

 
% position du ressaut : c'est quand conjuguée(y2) = y3)
xr = fzero(@(x) interp1(x3,y3,x) - interp1(x2,y2conjugue,x),40) ;
% tracé du point
scatter(xr, interp1(x3,y3,xr),'or','filled')
hold off
fprintf("Le ressaut se produit en x = %1.1f m avec h1 = %1.2f m et  h2 = %1.2f m.\n"...
,xr, interp1(x2,y2,xr), interp1(x3,y3,xr))



% Mes fonctions
function [value,isterminal,direction] = critique(x,y)
global hc1
value = (y - 0.99*hc1); % condition d'arrêt
isterminal = 1;  % si la condition est vérifiée, alors arrêt du calcul
direction = 0;   
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ydot = besse (x, y, slope, hn, hc) 
 
  ydot = slope*(1-(hn/y)^3)/(1-(hc/y)^3);
 
end