% Relation HQ
%
% Auteurs : Thibault Prévost.
% NOMA(s)   :0525 1500.
% Date    : 06/05/2021

clc;
clear all;
%close all;

DebitOpenfoam=[0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 3];
HauteurOpenFoam=[1.8335 1.9457 2.0529 2.1484 2.2379 2.3272 2.4081 2.4826 2.5445 2.6085 2.6745 2.7215];

%Dimensionnement
H0=0.5;
p=3*H0;

%valeurs théoriques
mu0=0.496;
%creager y/H0 = 0.47(x/H0)^(1.8)
L=1;
g=9.81;

Q0=mu0*L*sqrt(2*g)*H0^(3/2);

Hauteur=0.2:0.1:1.7;
Debit=length(Hauteur);

mu=mu0.*(Hauteur./H0).^0.17;

for i=1:length(Hauteur)
Debit(i)= mu(i) * L * sqrt(2*g) * Hauteur(i)^(3/2);
end

Hauteur=Hauteur+p;


%% Calcul valeurs expérimentales de mu

for i=1:length(DebitOpenfoam)

    muExpe(i)= DebitOpenfoam(i)/(L * sqrt(2*g) * (HauteurOpenFoam(i)-p)^(3/2));

end

%calcul du mu0 expériementale
mu0Expe=0;

for i=2:length(DebitOpenfoam)
    if HauteurOpenFoam(i)>2 && mu0Expe==0
        delta=(2-HauteurOpenFoam(i-1))/(HauteurOpenFoam(i)-HauteurOpenFoam(i-1));
        DebitDimensionnement=DebitOpenfoam(i-1)+delta*(DebitOpenfoam(i)-DebitOpenfoam(i-1));
        mu0Expe = DebitDimensionnement/(L * sqrt(2*g) * (2-p)^(3/2))
    else
    end    
end

Q0exp=mu0Expe*L*sqrt(2*g)*H0^(3/2);

%Calcul de la courbe avec la formule théorique (mu/mu0) = (H/H0)^0.17
muexpExtrapol=mu0Expe.*((HauteurOpenFoam-p)./0.5).^0.17;

%Calcul du nouvel exposant pour la formule théorique
expo= (log(muExpe./mu0Expe)./log((HauteurOpenFoam-p)./H0));
expoMoyen=mean(expo)

Exposigma=std(expo);

minExp=expoMoyen-2*Exposigma;
maxExp=expoMoyen+2*Exposigma;

%calcul de la nouvelle courbe de \mu
NouveauMu=mu0Expe.*((Hauteur-p)./H0).^expoMoyen;

figure
plot(DebitOpenfoam,HauteurOpenFoam-p,'ro')
hold on
plot(Debit,Hauteur-p)
xlabel('Debit [m^3 /s]')
ylabel('Hauteur H [m]')
legend('Relation Hauteur-Débit calculée avec le programme OpenFoam','Relation Hauteur-Débit théorique','Location','southeast')
hold off


figure
plot((HauteurOpenFoam-p)./H0,muExpe./mu0Expe,'ro')
hold on
plot((Hauteur-p)./H0,mu./mu0)
hold on
plot((Hauteur-p)./H0,NouveauMu./mu0Expe)
ylabel('\mu/\mu_0')
xlabel('Hauteur H/H_0')
legend('Valeurs de \mu_{expérimentale} calculé avec le programme OpenFoam','Relation \mu-Hauteur théorique','Relation \mu_{expérimentale}-Hauteur déduite des simulations','Location','southeast')
hold off