clearvars
close all


load('array_simb1.mat');
% load('array_exp.mat');




U_deb=0.3567405;

H=0.065;%en m



%% calcul modèle couche de mélange enchainement colline lambda=0.55


void_mat_l55=[];

[ivend_l55,jvend_l55]=size(array_sim(7).z(:,1:end));


void_mat_l55(1:ivend_l55,1:jvend_l55)=NaN;


mat_zcurv_l55=[void_mat_l55(1:66,:); array_sim(7).z(:,1:end)];

mat_ucurv_l55=[void_mat_l55(1:66,:); array_sim(7).U(:,1:end)];

mat_wcurv_l55=[void_mat_l55(1:66,:); array_sim(7).W(:,1:end)];

mat_kcurv_l55=[void_mat_l55(1:66,:); array_sim(7).k(:,1:end)];

A_l55=[];
B_l55=[];

for j_l55=1:jvend_l55
    
    b1_l55=array_sim(7).z(end,j_l55);

    A_l55=array_sim(7).z(end,j_l55);
    
    for i_l55=1:round(b1_l55/0.015384615)
     
        if A_l55(i_l55)>0.007856769230769%0.019975384615385
        
            A_l55(i_l55+1)=A_l55(i_l55)-0.015384615;
   
            B_l55(i_l55)=A_l55(i_l55);
        end
    
    
    end
    
    mat_zcurv_l55(:,j_l55)=[void_mat_l55(1:abs((length(A_l55(2:(end-1)))-66)),j_l55);mat_zcurv_l55(67:end,j_l55);A_l55(2:(end-1))'];
    mat_ucurv_l55(:,j_l55)=[void_mat_l55(1:abs((length(A_l55(2:(end-1)))-66)),j_l55);array_sim(7).U(:,j_l55);void_mat_l55(1:length(A_l55(2:(end-1))))'];
    mat_wcurv_l55(:,j_l55)=[void_mat_l55(1:abs((length(A_l55(2:(end-1)))-66)),j_l55);array_sim(7).W(:,j_l55);void_mat_l55(1:length(A_l55(2:(end-1))))'];
    mat_kcurv_l55(:,j_l55)=[void_mat_l55(1:abs((length(A_l55(2:(end-1)))-66)),j_l55);array_sim(7).k(:,j_l55);void_mat_l55(1:length(A_l55(2:(end-1))))'];

    clear A_l55 b1_l55
end


add_matx_l55=array_sim(7).x(1:66,1:end);
mat_xcurv_l55=[add_matx_l55; array_sim(7).x(:,1:end)];






% projection des profils sur la ligne de courant

xh1_l55 = array_sim(7).xst(1:end);
zh1_l55 = array_sim(7).zst(1:end);

iend_l55=length(xh1_l55);

Jxh1_l55=zeros(iend_l55,1);

ucurv_st_l55=zeros(length(mat_xcurv_l55(:,1)),iend_l55);
kcurv_st_l55=zeros(length(mat_xcurv_l55(:,1)),iend_l55);

for i_l55=1:10:5841
    
    

Jxh1_l55(i_l55)=find(xh1_l55(i_l55)-1e-3<mat_xcurv_l55(1,:) & mat_xcurv_l55(1,:)<xh1_l55(i_l55)+1e-3);

ucurv_st_l55(:,i_l55)=mat_ucurv_l55(:,Jxh1_l55(i_l55));

kcurv_st_l55(:,i_l55)=mat_kcurv_l55(:,Jxh1_l55(i_l55));

end


alpha_l55=zeros(iend_l55,1);

cos_alpha_l55=zeros(iend_l55,1);

u_proj_st_l55=zeros(length(mat_xcurv_l55(:,1)),iend_l55);

k_proj_st_l55=zeros(length(mat_xcurv_l55(:,1)),iend_l55);

z_proj_st_l55=zeros(length(mat_xcurv_l55(:,1)),iend_l55);


for ii_l55=11:10:5841

dxproj_l55=(mat_xcurv_l55(1,Jxh1_l55(ii_l55)+1)-mat_xcurv_l55(1,Jxh1_l55(ii_l55)-1));

dzproj_l55=(zh1_l55(ii_l55+1)-zh1_l55(ii_l55-1));

tan_alpha_l55=dzproj_l55/dxproj_l55;

alpha_l55(ii_l55)=atan(tan_alpha_l55);

cos_alpha_l55(ii_l55)=cos(alpha_l55(ii_l55));

u_proj_st_l55(:,ii_l55)=cos_alpha_l55(ii_l55).*ucurv_st_l55(:,ii_l55);

k_proj_st_l55(:,ii_l55)=cos_alpha_l55(ii_l55).*kcurv_st_l55(:,ii_l55);

z_proj_st_l55(:,ii_l55)=cos_alpha_l55(ii_l55).*mat_zcurv_l55(:,Jxh1_l55(ii_l55));

clear dzproj_l55 dxproj_l55 tan_alpha_l55
end 


% initialisation modèle suivre ligne de courant

x_proj_l55=mat_xcurv_l55(:,Jxh1_l55(11:10:5841)).*H;
z_proj_l55=z_proj_st_l55(:,11:10:5841).*H;

u_proj_l55=u_proj_st_l55(:,11:10:5841).*U_deb;
k_proj_l55=k_proj_st_l55(:,11:10:5841).*(U_deb^2);

z_st_l55=zh1_l55(:,11:10:5841).*H;                                                                                                                                  
e_coord_l55=zeros(length(u_proj_l55(:,1)),length(u_proj_l55(1,:)));
zm_l55=zeros(length(u_proj_l55(1,:)),1);

Uscale_l55=zeros(length(u_proj_l55(:,1)),length(u_proj_l55(1,:)));
kscale_l55=zeros(length(u_proj_l55(:,1)),length(u_proj_l55(1,:)));





% initialisation modèle Pope turbulent flow (2000)

xcurv_pope=mat_xcurv_l55(1,Jxh1_l55(1):end).*H;

zcurv_pope=mat_zcurv_l55(:,Jxh1_l55(1):end).*H;

u_pope_l55=mat_ucurv_l55(:,Jxh1_l55(1):end).*U_deb;

k_pope_l55=mat_kcurv_l55(:,Jxh1_l55(1):end).*(U_deb^2);

                                                                                                        
zm_l55_pope=zeros(length(u_pope_l55(1,:)),1);
z_09a_l55_pope=zeros(length(u_pope_l55(1,:)),1);
z_01a_l55_pope=zeros(length(u_pope_l55(1,:)),1);
e_coord_l55_pope=zeros(length(u_pope_l55(:,1)),length(u_pope_l55(1,:)));

Uscale_l55_pope=zeros(length(u_pope_l55(:,1)),length(u_pope_l55(1,:)));
kscale_l55_pope=zeros(length(u_pope_l55(:,1)),length(u_pope_l55(1,:)));



% initialisation commune pour les deux modèles


z_09a_l55=zeros(length(u_proj_l55(1,:)),1);
z_01a_l55=zeros(length(u_proj_l55(1,:)),1);
delta_l55=zeros(length(u_proj_l55(1,:)),1);


% vitesse initial de la couche de mélange pour les deux modèles
Ub_l55=0;

Ut_l55=u_pope_l55(155,1); 



% vitesse caractéristique de la couche de mélange pour les deux modèles

Uc_l55=0.5*(Ut_l55+Ub_l55);

Us_l55=Ut_l55-Ub_l55;

U_09a_l55=Ub_l55+0.9*Us_l55;

U_01a_l55=Ub_l55+0.1*Us_l55;








for iv_l55=1:length(u_proj_l55(1,:))



idelta_l55=find(U_01a_l55<=u_proj_l55(:,iv_l55) & u_proj_l55(:,iv_l55)<=U_09a_l55);


z_09a_l55(iv_l55)=z_proj_l55(min(idelta_l55),iv_l55);

z_01a_l55(iv_l55)=z_proj_l55(max(idelta_l55),iv_l55);

delta_l55(iv_l55)=z_09a_l55(iv_l55)-z_01a_l55(iv_l55);

%% methode suivre la ligne de courant

zm_l55(iv_l55)=z_st_l55(iv_l55);

e_coord_l55(:,iv_l55)=(z_proj_l55(:,iv_l55)-zm_l55(iv_l55))./delta_l55(iv_l55);

Uscale_l55(:,iv_l55)=(u_proj_l55(:,iv_l55)-Uc_l55)./Us_l55;

kscale_l55(:,iv_l55)=k_proj_l55(:,iv_l55)./(Us_l55^2);


%% methode Pope turbulent flow (2000)

zm_l55_pope(iv_l55)=0.5*(z_09a_l55(iv_l55)+z_01a_l55(iv_l55));

e_coord_l55_pope(:,iv_l55)=(zcurv_pope(:,iv_l55)-zm_l55_pope(iv_l55))./delta_l55(iv_l55);

Uscale_l55_pope(:,iv_l55)=(u_pope_l55(:,iv_l55)-Uc_l55)./Us_l55;

kscale_l55_pope(:,iv_l55)=k_pope_l55(:,iv_l55)./(Us_l55.^2);




clear idelta_l55
end



%% tracé
y_lim = [-0.9 0.8];
x_lim = [-1 1.2];

F = figure('rend','painters','pos',[10 10 1200 1200],'DefaultAxesFontName','cmr10','defaulttextinterpreter','latex')

% fonction erreur et paramètre de fitting
a_uscale=-0.46;
b_uscale=1.65;

Fit_Uscale=(a_uscale+erf(e_coord_l55(1:5:end,306)))/b_uscale;


% plot de la méthode suivre la ligne de courant

Ax1=subplot(2,1,1);

% exemple set(Ax,'Position', [x y width, height])
set(Ax1,'Position',[0.1 0.62 0.65 0.35])

box on
hold on

plot(e_coord_l55(1:5:end,306),Uscale_l55(1:5:end,306),'kp');
plot(e_coord_l55(1:5:end,360),Uscale_l55(1:5:end,360),'kx');
plot(e_coord_l55(1:5:end,414),Uscale_l55(1:5:end,414),'kh');
plot(e_coord_l55(1:5:end,468),Uscale_l55(1:5:end,468),'kv');
plot(e_coord_l55(1:5:end,522),Uscale_l55(1:5:end,522),'k^');
plot(e_coord_l55(1:5:end,584),Uscale_l55(1:5:end,584),'ks');

plot(e_coord_l55(1:5:end,306),Fit_Uscale,'r-');



    hy = ylabel('$F=(\overline{u}_{st}-U_{c})/U_{s}$','FontSize',30);
    hx = xlabel('$\zeta$','FontSize',30);
    set(hx,'interpreter','Latex','FontSize',15)
    set(hy,'interpreter','Latex','FontSize',15)

    xlim([x_lim(1) x_lim(2)])
    ylim([y_lim(1) y_lim(2)])

    set(gca, 'xTick', [-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1])
    set(gca, 'yTick', [-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8])
    set(gca, 'FontSize', 20);
    set(gca,'TickLabelInterpreter', 'Latex');


leg=legend(['$x/H_{1}=3.6$'],...
['$x/H_{1}=4.2$'],...
['$x/H_{1}=4.8$'],...
['$x/H_{1}=5.4$'],...
['$x/H_{1}=6$'],...
['$x/H_{1}=6.7$'],...
['$error \; function$']);
set(leg,'Interpreter','Latex','Fontsize',20);


% position légende
% exemple pos_leg =[x y width, height];
pos_leg = [0.773, 0.759, .16, .18];
set(leg, 'Position', pos_leg)

% position lettre a
pos_letter_x = 0.06;
pos_letter_y = -1.4;
    letter = '(a)';
    text(pos_letter_x,pos_letter_y,letter,'Interpreter','latex','FontSize',20)



% plot de la méthode Pope turbulent flow (2000)


% fonction erreur et paramètre de fitting
a_uscale_pope=0;
b_uscale_pope=1.4;


Fit_Uscale_pope=(a_uscale_pope+erf(e_coord_l55_pope(1:5:end,306)))/b_uscale_pope;


Ax2=subplot(2,1,2);

% exemple set(Ax,'Position', [x y width, height])
set(Ax2,'Position',[0.1 0.13 0.65 0.35])

box on
hold on

plot(e_coord_l55_pope(1:5:end,307),Uscale_l55_pope(1:5:end,307),'kp');
plot(e_coord_l55_pope(1:5:end,361),Uscale_l55_pope(1:5:end,361),'kx');
plot(e_coord_l55_pope(1:5:end,415),Uscale_l55_pope(1:5:end,415),'kh');
plot(e_coord_l55_pope(1:5:end,469),Uscale_l55_pope(1:5:end,469),'kv');
plot(e_coord_l55_pope(1:5:end,523),Uscale_l55_pope(1:5:end,523),'k^');
plot(e_coord_l55_pope(1:5:end,573),Uscale_l55_pope(1:5:end,573),'ks');

plot(e_coord_l55_pope(1:5:end,307),Fit_Uscale_pope,'r-');



    hy = ylabel('$F=(\overline{u}-U_{c})/U_{s}$','FontSize',30);
    hx = xlabel('$\zeta$','FontSize',30);
    set(hx,'interpreter','Latex','FontSize',15)
    set(hy,'interpreter','Latex','FontSize',15)

    xlim([x_lim(1) x_lim(2)])
    ylim([y_lim(1) y_lim(2)])

    set(gca, 'xTick', [-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1])
    set(gca, 'yTick', [-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8])
    set(gca, 'FontSize', 20);
    set(gca,'TickLabelInterpreter', 'Latex');

% position lettre b

pos_letter_x2 = 0.06;
pos_letter_y2 = -1.4;
    letter2 = '(b)';
    text(pos_letter_x2,pos_letter_y2,letter2,'Interpreter','latex','FontSize',20)

