clearvars
close all

run('param_exp_general.m');
load('res_Q6Q9Q14_array.mat');

S_totale = 0.045*0.045;

%% Calcul de l'integral de la porosité dans la canopée
phi_PHem40_04 =@(z) 1 - (pi * (R_hemi^2 - (z + R_hemi-h_04).^2))/S_totale;
phi_PHem40_08 =@(z) 1 - (pi * (R_hemi^2 - (z + R_hemi-h_08).^2))/S_totale;
phi_PHem40_12 =@(z) 1 - (pi * (R_hemi^2 - (z + R_hemi-h_12).^2))/S_totale;
phi_PHem40_16 =@(z) 1 - (pi * (R_hemi^2 - (z + R_hemi-h_16).^2))/S_totale;

int_Phi_PHem40_04 = integral(phi_PHem40_04,0,h_04);
int_Phi_PHem40_08 = integral(phi_PHem40_08,0,h_08);
int_Phi_PHem40_12 = integral(phi_PHem40_12,0,h_12);
int_Phi_PHem40_16 = integral(phi_PHem40_16,0,h_16);

for pp = 1:100
    
    hh(pp) = pp/5000;
    phi =@(z) 1 - (pi * (R_hemi^2 - (z + R_hemi-hh(pp)).^2))/S_totale;
    int_Phi(pp) = integral(phi,0,hh(pp));
    
    clear phi
end
PP = 1:100;



%% Détermination du delta
% Q14
fit_Q14_Treynolds_FR = polyfit(z(90:140)',myarray(1).tau_DA(90:140),1);
yfit_Q14_FR = fit_Q14_Treynolds_FR(1)*z+fit_Q14_Treynolds_FR(2);
delta_Q14_FR = -fit_Q14_Treynolds_FR(2)/fit_Q14_Treynolds_FR(1) - myarray(1).zsmoy_P1;

fit_Treynolds_Q14_PHem40_04 = polyfit(z(59:99)',myarray(2).tau_DA(59:99),1);
yfit_Q14_PHem40_04 = fit_Treynolds_Q14_PHem40_04(1)*z+fit_Treynolds_Q14_PHem40_04(2);
delta_Q14_PHem40_04 = -fit_Treynolds_Q14_PHem40_04(2)/fit_Treynolds_Q14_PHem40_04(1) - myarray(2).zsmoy_P1 ;

fit_Treynolds_Q14_PHem40_08 = polyfit(z(76:130)',myarray(3).tau_DA(76:130),1);
yfit_Q14_PHem40_08 = fit_Treynolds_Q14_PHem40_08(1)*z+fit_Treynolds_Q14_PHem40_08(2);
delta_Q14_PHem40_08 = -fit_Treynolds_Q14_PHem40_08(2)/fit_Treynolds_Q14_PHem40_08(1) - myarray(3).zsmoy_P1;

fit_Treynolds_Q14_PHem40_12 = polyfit(z(107:151)',myarray(4).tau_DA(107:151),1);
yfit_Q14_PHem40_12 = fit_Treynolds_Q14_PHem40_12(1)*z+fit_Treynolds_Q14_PHem40_12(2);
delta_Q14_PHem40_12 = -fit_Treynolds_Q14_PHem40_12(2)/fit_Treynolds_Q14_PHem40_12(1) - myarray(4).zsmoy_P1;

fit_Treynolds_Q14_PHem40_16 = polyfit(z(136:178)',myarray(5).tau_DA(136:178),1);
yfit_Q14_PHem40_16 = fit_Treynolds_Q14_PHem40_16(1)*z+fit_Treynolds_Q14_PHem40_16(2);
delta_Q14_PHem40_16 = -fit_Treynolds_Q14_PHem40_16(2)/fit_Treynolds_Q14_PHem40_16(1) - myarray(5).zsmoy_P1;

%Q9
fit_Treynolds_Q9_FR = polyfit(myarray(6).z_P1(70:140)',myarray(6).tau_DA(70:140),1);
yfit_Q9_FR = fit_Treynolds_Q9_FR(1)*myarray(6).z_P1+fit_Treynolds_Q9_FR(2);
delta_Q9_FR = -fit_Treynolds_Q9_FR(2)/fit_Treynolds_Q9_FR(1) - myarray(1).zsmoy_P1;

fit_Treynolds_Q9_PHem40_04 = polyfit(myarray(7).z_P1(59:99)',myarray(7).tau_DA(59:99),1);
yfit_Q9_PHem40_04 = fit_Treynolds_Q9_PHem40_04(1)*myarray(7).z_P1+fit_Treynolds_Q9_PHem40_04(2);
delta_Q9_PHem40_04 = -fit_Treynolds_Q9_PHem40_04(2)/fit_Treynolds_Q9_PHem40_04(1) - myarray(7).zsmoy_P1;

fit_Treynolds_Q9_PHem40_08 = polyfit(myarray(8).z_P1(76:130)',myarray(8).tau_DA(76:130),1);
yfit_Q9_PHem40_08 = fit_Treynolds_Q9_PHem40_08(1)*myarray(8).z_P1+fit_Treynolds_Q9_PHem40_08(2);
delta_Q9_PHem40_08 = -fit_Treynolds_Q9_PHem40_08(2)/fit_Treynolds_Q9_PHem40_08(1) - myarray(8).zsmoy_P1;

fit_Treynolds_Q9_PHem40_12 = polyfit(myarray(9).z_P1(107:151)',myarray(9).tau_DA(107:151),1);
yfit_Q9_PHem40_12 = fit_Treynolds_Q9_PHem40_12(1)*myarray(9).z_P1+fit_Treynolds_Q9_PHem40_12(2);
delta_Q9_PHem40_12 = -fit_Treynolds_Q9_PHem40_12(2)/fit_Treynolds_Q9_PHem40_12(1) - myarray(9).zsmoy_P1;

fit_Treynolds_Q9_PHem40_16 = polyfit(myarray(10).z_P1(107:151)',myarray(10).tau_DA(107:151),1);
yfit_Q9_PHem40_16 = fit_Treynolds_Q9_PHem40_16(1)*myarray(10).z_P1+fit_Treynolds_Q9_PHem40_16(2);
delta_Q9_PHem40_16 = -fit_Treynolds_Q9_PHem40_16(2)/fit_Treynolds_Q9_PHem40_16(1) - myarray(10).zsmoy_P1;

% Q6
fit_Q6_Treynolds_FR = polyfit(z(90:140)',myarray(11).tau_DA(90:140),1);
yfit_Q6_FR = fit_Q6_Treynolds_FR(1)*z+fit_Q6_Treynolds_FR(2);
delta_Q6_FR = -fit_Q6_Treynolds_FR(2)/fit_Q6_Treynolds_FR(1) - myarray(1).zsmoy_P1;

fit_Treynolds_Q6_PHem40_04 = polyfit(z(59:99)',myarray(12).tau_DA(59:99),1);
yfit_Q6_PHem40_04 = fit_Treynolds_Q6_PHem40_04(1)*z+fit_Treynolds_Q6_PHem40_04(2);
delta_Q6_PHem40_04 = -fit_Treynolds_Q6_PHem40_04(2)/fit_Treynolds_Q6_PHem40_04(1) - myarray(2).zsmoy_P1;

fit_Treynolds_Q6_PHem40_08 = polyfit(z(76:130)',myarray(13).tau_DA(76:130),1);
yfit_Q6_PHem40_08 = fit_Treynolds_Q6_PHem40_08(1)*z+fit_Treynolds_Q6_PHem40_08(2);
delta_Q6_PHem40_08 = -fit_Treynolds_Q6_PHem40_08(2)/fit_Treynolds_Q6_PHem40_08(1) - myarray(3).zsmoy_P1;

fit_Treynolds_Q6_PHem40_12 = polyfit(z(107:151)',myarray(14).tau_DA(107:151),1);
yfit_Q6_PHem40_12 = fit_Treynolds_Q6_PHem40_12(1)*z+fit_Treynolds_Q6_PHem40_12(2);
delta_Q6_PHem40_12 = -fit_Treynolds_Q6_PHem40_12(2)/fit_Treynolds_Q6_PHem40_12(1) - myarray(4).zsmoy_P1;

fit_Treynolds_Q6_PHem40_16 = polyfit(z(107:151)',myarray(15).tau_DA(107:151),1);
yfit_Q6_PHem40_16 = fit_Treynolds_Q6_PHem40_16(1)*z+fit_Treynolds_Q6_PHem40_16(2);
delta_Q6_PHem40_16 = -fit_Treynolds_Q6_PHem40_16(2)/fit_Treynolds_Q6_PHem40_16(1) - myarray(5).zsmoy_P1;


%% Calcul de la force totale

tautot_Q14_FR = myarray(1).ustar_fit^2;
tautot_Q14_PHem40_04 = myarray(2).ustar_fit^2 * (1 + int_Phi_PHem40_04/delta_Q14_PHem40_04);
tautot_Q14_PHem40_08 = myarray(3).ustar_fit^2 * (1 + int_Phi_PHem40_08/delta_Q14_PHem40_08);
tautot_Q14_PHem40_12 = myarray(4).ustar_fit^2 * (1 + int_Phi_PHem40_12/delta_Q14_PHem40_12);
tautot_Q14_PHem40_16 = myarray(5).ustar_fit^2 * (1 + int_Phi_PHem40_16/delta_Q14_PHem40_16);

tautot_Q9_FR = myarray(6).ustar_fit^2;
tautot_Q9_PHem40_04 = myarray(7).ustar_fit^2 * (1 + int_Phi_PHem40_04/delta_Q9_PHem40_04);
tautot_Q9_PHem40_08 = myarray(8).ustar_fit^2 * (1 + int_Phi_PHem40_08/delta_Q9_PHem40_08);
tautot_Q9_PHem40_12 = myarray(9).ustar_fit^2 * (1 + int_Phi_PHem40_12/delta_Q9_PHem40_12);
tautot_Q9_PHem40_16 = myarray(10).ustar_fit^2 * (1 + int_Phi_PHem40_16/delta_Q9_PHem40_16);

tautot_Q6_FR = myarray(11).ustar_fit^2;
tautot_Q6_PHem40_04 = myarray(12).ustar_fit^2 * (1 + int_Phi_PHem40_04/delta_Q6_PHem40_04);
tautot_Q6_PHem40_08 = myarray(13).ustar_fit^2 * (1 + int_Phi_PHem40_08/delta_Q6_PHem40_08);
tautot_Q6_PHem40_12 = myarray(14).ustar_fit^2 * (1 + int_Phi_PHem40_12/delta_Q6_PHem40_12);
tautot_Q6_PHem40_16 = myarray(15).ustar_fit^2 * (1 + int_Phi_PHem40_16/delta_Q6_PHem40_16);

tautot_Q9 = [tautot_Q9_FR tautot_Q9_PHem40_04 tautot_Q9_PHem40_08 tautot_Q9_PHem40_12 tautot_Q9_PHem40_16];
tautot_Q6 = [tautot_Q6_FR tautot_Q6_PHem40_04 tautot_Q6_PHem40_08 tautot_Q6_PHem40_12 tautot_Q6_PHem40_16];
tautot_Q14 = [tautot_Q14_FR tautot_Q14_PHem40_04 tautot_Q14_PHem40_08 tautot_Q14_PHem40_12 tautot_Q14_PHem40_16];


%% Tracé
deltap = [delta_Q14_FR delta_Q14_PHem40_04 delta_Q14_PHem40_08 delta_Q14_PHem40_12 delta_Q14_PHem40_16];

ustar_fit_Q14 = [myarray(1).ustar_fit myarray(2).ustar_fit myarray(3).ustar_fit myarray(4).ustar_fit myarray(5).ustar_fit];
ustar_fit_Q9 = [myarray(6).ustar_fit myarray(7).ustar_fit myarray(8).ustar_fit myarray(9).ustar_fit myarray(10).ustar_fit];
ustar_fit_Q6 = [myarray(11).ustar_fit myarray(12).ustar_fit myarray(13).ustar_fit myarray(14).ustar_fit myarray(15).ustar_fit];

P = [0 20 40 60 80];

figure('rend','painters','pos',[10 10 1050 400])
ha = tight_subplot(1,2,[.1 .1],[.18 .03],[.07 .07])  % (nombre graphe vertical, nombre graphe horizontal,[],[distance bord bas, distance bord haut],[distance bord gauche distance bord droit])

axes(ha(1))
plot(P,tautot_Q14./tautot_Q14(1),'ks','MarkerFaceColor','k')
hold on
plot(P,tautot_Q9./tautot_Q9(1),'ko')
plot(P,tautot_Q6./tautot_Q6(1),'k^','MarkerFaceColor','k')
% plot([-5 100],[1 1],'k--','LineWidth',0.5)
% plot(P,ustar_maxTR_Q6,'k^:','MarkerFaceColor','k')
% plot(P,ustar_fit_Q6,'ks:','MarkerFaceColor','k')
hl = legend('$Re_{D} = 27700$','$Re_{D} = 18600$','$Re_{D} = 13300$','Location','southeast');
hy = ylabel('$\tau_{b}/\tau_{P0}$','FontSize',30);
hx = xlabel('$P (\%)$','FontSize',30);
set(gca, 'FontSize', 15);
set(hl,'interpreter','Latex','FontSize',14)
set(hx,'interpreter','Latex','FontSize',20)
set(hy,'interpreter','Latex','FontSize',20)
set(gca,'TickLabelInterpreter', 'Latex');
xlim([0 100])
ylim([1 4])
box on
letter = '(a)';
text(90,3.75,letter,'Interpreter','latex','FontSize',18)
box on


axes(ha(2))
plot(P,tautot_Q14./ustar_fit_Q14.^2,'ks','MarkerFaceColor','k')
hold on
plot(P,tautot_Q9./ustar_fit_Q9.^2,'ko')
plot(P,tautot_Q6./ustar_fit_Q6.^2,'k^','MarkerFaceColor','k')
% plot([-5 100],[1 1],'k--','LineWidth',0.5)
% plot(P,ustar_maxTR_Q6,'k^:','MarkerFaceColor','k')
% plot(P,ustar_fit_Q6,'ks:','MarkerFaceColor','k')
hl = legend('$Re_{D} = 27700$','$Re_{D} = 18600$','$Re_{D} = 13300$','Location','southeast');
hy = ylabel('$\tau_{b}/\rho_{f}u_{\star}^2$','FontSize',30);
hx = xlabel('$P (\%)$','FontSize',30);
set(gca, 'FontSize', 15);
set(hl,'interpreter','Latex','FontSize',14)
set(hx,'interpreter','Latex','FontSize',20)
set(hy,'interpreter','Latex','FontSize',20)
set(gca,'TickLabelInterpreter', 'Latex');
xlim([0 100])
ylim([1 1.6])
box on
letter = '(b)';
text(90,1.55,letter,'Interpreter','latex','FontSize',18)
box on


axes('Position',[.605 .60 .16 .35])
box on
plot(PP,int_Phi./hh,'Color',[0,0.45,0.74],'LineWidth',1)
hy = ylabel('$\int_{0}^{k}\phi(z)dz/k$','FontSize',20);
set(gca, 'FontSize', 10);
set(hy,'interpreter','Latex','FontSize',13)
set(gca,'TickLabelInterpreter', 'Latex');
xlim([0 100])
ylim([0.5 1])

% 
% F = figure
% plot(PP,int_Phi./hh,'Color',[0,0.45,0.74],'LineWidth',3)
% hold on
% hy = ylabel('$\int_{0}^{k}\phi(z)dz/k$','FontSize',30);
% set(gca, 'FontSize', 25);
% set(hy,'interpreter','Latex','FontSize',30)
% set(gca,'TickLabelInterpreter', 'Latex');
% xlim([0 100])
% ylim([0.5 1])
% box on
% tightfig(F)