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.84


void_mat_l84=[];

[ivend_l84,jvend_l84]=size(array_sim(11).z(:,1:end));


void_mat_l84(1:ivend_l84,1:jvend_l84)=NaN;


mat_zcurv_l84=[void_mat_l84(1:66,:); array_sim(11).z(:,1:end)];

mat_ucurv_l84=[void_mat_l84(1:66,:); array_sim(11).U(:,1:end)];

mat_wcurv_l84=[void_mat_l84(1:66,:); array_sim(11).W(:,1:end)];

mat_kcurv_l84=[void_mat_l84(1:66,:); array_sim(11).k(:,1:end)];

A_l84=[];
B_l84=[];

for j_l84=1:jvend_l84
    
    b1_l84=array_sim(11).z(end,j_l84);

    A_l84=array_sim(11).z(end,j_l84);
    
    for i_l84=1:round(b1_l84/0.015384615)
     
        if A_l84(i_l84)>0.007856769230769%0.019975384615385
        
            A_l84(i_l84+1)=A_l84(i_l84)-0.015384615;
   
            B_l84(i_l84)=A_l84(i_l84);
        end
    
    
    end
    
    mat_zcurv_l84(:,j_l84)=[void_mat_l84(1:abs((length(A_l84(2:(end-1)))-66)),j_l84);mat_zcurv_l84(67:end,j_l84);A_l84(2:(end-1))'];
    mat_ucurv_l84(:,j_l84)=[void_mat_l84(1:abs((length(A_l84(2:(end-1)))-66)),j_l84);array_sim(11).U(:,j_l84);void_mat_l84(1:length(A_l84(2:(end-1))))'];
    mat_wcurv_l84(:,j_l84)=[void_mat_l84(1:abs((length(A_l84(2:(end-1)))-66)),j_l84);array_sim(11).W(:,j_l84);void_mat_l84(1:length(A_l84(2:(end-1))))'];
    mat_kcurv_l84(:,j_l84)=[void_mat_l84(1:abs((length(A_l84(2:(end-1)))-66)),j_l84);array_sim(11).k(:,j_l84);void_mat_l84(1:length(A_l84(2:(end-1))))'];

    clear A_l84 b1_l84
end


add_matx_l84=array_sim(11).x(1:66,1:end);
mat_xcurv_l84=[add_matx_l84; array_sim(11).x(:,1:end)];

xh1_l84 = array_sim(11).xst(1:end);
zh1_l84 = array_sim(11).zst(1:end);

iend_l84=length(xh1_l84);

Jxh1_l84=zeros(iend_l84,1);

ucurv_st_l84=zeros(length(mat_xcurv_l84(:,1)),iend_l84);

kcurv_st_l84=zeros(length(mat_xcurv_l84(:,1)),iend_l84);

for i_l84=1:10:5201
    
    

Jxh1_l84(i_l84)=find(xh1_l84(i_l84)-3e-3<mat_xcurv_l84(1,:) & mat_xcurv_l84(1,:)<xh1_l84(i_l84)+3e-3);

ucurv_st_l84(:,i_l84)=mat_ucurv_l84(:,Jxh1_l84(i_l84));

kcurv_st_l84(:,i_l84)=mat_kcurv_l84(:,Jxh1_l84(i_l84));

end

% initialisation modèle Pope turbulent flow (2000)

xcurv_pope=mat_xcurv_l84(1,Jxh1_l84(1):end).*H;

zcurv_pope=mat_zcurv_l84(:,Jxh1_l84(1):end).*H;

u_pope_l84=mat_ucurv_l84(:,Jxh1_l84(1):end).*U_deb;

k_pope_l84=mat_kcurv_l84(:,Jxh1_l84(1):end).*(U_deb^2);

                                                                                                        
zm_l84_pope=zeros(length(u_pope_l84(1,:)),1);
z_09a_l84_pope=zeros(length(u_pope_l84(1,:)),1);
z_01a_l84_pope=zeros(length(u_pope_l84(1,:)),1);
e_coord_l84_pope=zeros(length(u_pope_l84(:,1)),length(u_pope_l84(1,:)));

Uscale_l84_pope=zeros(length(u_pope_l84(:,1)),length(u_pope_l84(1,:)));
kscale_l84_pope=zeros(length(u_pope_l84(:,1)),length(u_pope_l84(1,:)));

delta_l84=zeros(length(u_pope_l84(1,:)),1);


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

Ut_l84=u_pope_l84(157,1); 



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

Uc_l84=0.5*(Ut_l84+Ub_l84);

Us_l84=Ut_l84-Ub_l84;

U_09a_l84=Ub_l84+0.9*Us_l84;

U_01a_l84=Ub_l84+0.1*Us_l84;




for iv_l84=1:length(u_pope_l84(1,:))



idelta_l84=find(U_01a_l84<=u_pope_l84(:,iv_l84) & u_pope_l84(:,iv_l84)<=U_09a_l84);


z_09a_l84_pope(iv_l84)=zcurv_pope(min(idelta_l84),iv_l84);

z_01a_l84_pope(iv_l84)=zcurv_pope(max(idelta_l84),iv_l84);

delta_l84(iv_l84)=z_09a_l84_pope(iv_l84)-z_01a_l84_pope(iv_l84);



% methode Pope turbulent flow (2000)

zm_l84_pope(iv_l84)=0.5*(z_09a_l84_pope(iv_l84)+z_01a_l84_pope(iv_l84));

e_coord_l84_pope(:,iv_l84)=(zcurv_pope(:,iv_l84)-zm_l84_pope(iv_l84))./delta_l84(iv_l84);

Uscale_l84_pope(:,iv_l84)=(u_pope_l84(:,iv_l84)-Uc_l84)./Us_l84;

kscale_l84_pope(:,iv_l84)=k_pope_l84(:,iv_l84)./(Us_l84.^2);




clear idelta_l84
end


%% tracé

y_lim = [0 0.08];
x_lim = [-1 1.2];



fig=figure('rend','painters','pos',[10 10 1500 600]);

% exemple set(axes,'Position', [x y width, height])
set(axes,'Position',[0.1 0.2 0.7 0.7])

box on
hold on

plot(e_coord_l84_pope(1:5:end,302),kscale_l84_pope(1:5:end,302),'kp');
plot(e_coord_l84_pope(1:5:end,353),kscale_l84_pope(1:5:end,353),'kx');
plot(e_coord_l84_pope(1:5:end,403),kscale_l84_pope(1:5:end,403),'kh');
plot(e_coord_l84_pope(1:5:end,453),kscale_l84_pope(1:5:end,453),'kv');
plot(e_coord_l84_pope(1:5:end,503),kscale_l84_pope(1:5:end,503),'k^');
plot(e_coord_l84_pope(1:5:end,554),kscale_l84_pope(1:5:end,554),'ks');

plot(e_coord_l84_pope(1:5:end,604),kscale_l84_pope(1:5:end,604),'k>');
plot(e_coord_l84_pope(1:5:end,654),kscale_l84_pope(1:5:end,654),'k.');


plot(e_coord_l84_pope(1:5:end,755),kscale_l84_pope(1:5:end,755),'k*');
plot(e_coord_l84_pope(1:5:end,838),kscale_l84_pope(1:5:end,838),'k+');
plot(e_coord_l84_pope(1:5:end,905),kscale_l84_pope(1:5:end,905),'ko');
plot(e_coord_l84_pope(1:5:end,1006),kscale_l84_pope(1:5:end,1006),'kd');

plot([0 0],[-1 1],'k--');


    hy = ylabel('$k/U_{s}^2$','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 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08])
    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.6$'],...
['$x/H_{1}=7.2$'],...
['$x/H_{1}=7.8$'],...
['$x/H_{1}=9$'],...
['$x/H_{1}=10$'],...
['$x/H_{1}=11$'],...
['$x/H_{1}=12$']);

set(leg,'Interpreter','Latex','Fontsize',20);
pos_letter_x = 20;
pos_letter_y = 1.02;
    letter = '$\lambda/H_{1}=12.9$';
    text(pos_letter_x,pos_letter_y,letter,'Interpreter','latex','FontSize',30)    

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