%% Getto secondario ON, tutto in HIGH SPEED NUMERICS 
close all
clear
clc

gamma=1.4; %aria
delta=(gamma-1)*0.5;
R=287; %cost. gas aria [J/(kg*K)]
c1=(gamma+1)/(2*(gamma-1));
c2=gamma/(gamma-1);
Gamma=sqrt(gamma)*(2/(gamma+1))^c1;
T0=300; %[K] temperatura totale aka Tc

L=0.10435; %[m] lunghezza divergente
rt=1; %[cm]
re=4.472; %[cm]
s=0.1*rt; %[cm] spessore parete ugello, sezione exit
Ae2=pi*s*(2*re+s); %[cm^2] sezione corona circolare

pa=101325;
p0_10=10*pa;
p0_12=12*pa;
p0_14=14*pa;
p0_16=16*pa;
p0_18=18*pa;
p0_20=20*pa;
p0_25=25*pa;
p0_30=30*pa;
p0_35=35*pa;
p0_40=40*pa;
p0_45=45*pa;
p0_50=50*pa;

pw_EUL=readmatrix("pw_wall_inv.txt");  % 120x80 Valentyn

% Portate limite, NPR=50 come riferimento

% Ho scelto NPR=50 perchè in queste condizioni il punto di separazione 
% ha superato il 95% della lunghezza L del divergente, perciò sono quasi
% condizioni full flowing ed ha senso prendere una portata di riferimento
% qui piuttosto che ad NPR=10 dove sarebbe troppo piccola


m_dot_50=3.631593; %[kg/s] misurata, vedi simulazione viscida, solo ugello
m_dot_lim=0.05*m_dot_50; %[kg/s] portata lim1 getto 2

Me2=[2,1.8,1.6,1]; % Mach d'uscita pari a 2 


pe2=pa; %[Pa] getto2 secondario adattato
p02=pe2*(1+delta*Me2.^2).^c2; %[Pa] pressione di camera getto2
npr2=p02./pa

m_dot2_1D=(Ae2*10^(-4))*(p02./sqrt(T0)).*(sqrt(gamma/R)).*Me2...
    .*(1+delta*Me2.^2).^(-c1); %[kg/s] prevista da approx 1D, vedi NASA

% Portate misurate con getto adattato Me2=[2,1.8,1.6,1]
m_dot2_mis=[0.2320778,0.20887,0.1856623,0.1160389]; %[kg/s] 

% Vediamo di quanto sforano il 5% della portata di riferimento
ratio_1D=m_dot2_1D./m_dot_50 % molto grossolana come stima
ratio_mis=m_dot2_mis./m_dot_50

% Getto sottoespanso (pe>pa)

p02_sot=p02(3);  % pressione  di camera del caso Me=1.6 (quindi stesso npr)
npr_sot=p02_sot/pa
pe_sot=pa+0.2*pa %[Pa] getto sottoespanso, +20% rispetto pa
k1=2/(gamma-1);
k2=(gamma-1)/gamma;
Me_sot=sqrt(k1*[(p02_sot/pe_sot)^k2-1])
m_dot_sot=0.2054167; %[kg/s] valore misurato da Fluent
ratio_sot=m_dot_sot/m_dot_50;


% Confronto pressioni

% Lettura pressioni di parete, HSN senza getto secondario

pw_data=leggi_pw('C:\Users\ricca\Documents\MATLAB\Lavori_tesi\TOP_TESI\pw_hsn');
pw6_10=pw_data{1};
pw6_12=pw_data{2};
pw6_14=pw_data{3};
pw6_16=pw_data{4};
pw6_18=pw_data{5};
pw6_20=pw_data{6};
pw6_22=pw_data{7};
pw6_25=pw_data{8};
pw6_30=pw_data{9};
pw6_35=pw_data{10};
pw6_40=pw_data{11};
pw6_45=pw_data{12};
pw6_50=pw_data{13};

% Lettura pressioni di parete, Me2=2 (Mach getto secondario = 2)

pw_data2=leggi_pw('C:\Users\ricca\Documents\MATLAB\Lavori_tesi\TOP_TESI\pw_me2');
pw6_10_me2=pw_data2{1};
pw6_12_me2=pw_data2{2};
pw6_14_me2=pw_data2{3};
pw6_16_me2=pw_data2{4};
pw6_18_me2=pw_data2{5};
pw6_20_me2=pw_data2{6};
pw6_25_me2=pw_data2{7};
pw6_30_me2=pw_data2{8};

% Lettura pressioni di parete, Me2=1.8 

pw_data18=leggi_pw('C:\Users\ricca\Documents\MATLAB\Lavori_tesi\TOP_TESI\pw_me18');
pw6_10_me18=pw_data18{1};
pw6_12_me18=pw_data18{2};
pw6_14_me18=pw_data18{3};
pw6_16_me18=pw_data18{4};
pw6_18_me18=pw_data18{5};
pw6_20_me18=pw_data18{6};
pw6_25_me18=pw_data18{7};
pw6_30_me18=pw_data18{8};

% Lettura pressioni di parete, Me2=1.6 

pw_data16=leggi_pw('C:\Users\ricca\Documents\MATLAB\Lavori_tesi\TOP_TESI\pw_me16');
pw6_10_me16=pw_data16{1};
pw6_12_me16=pw_data16{2};
pw6_14_me16=pw_data16{3};
pw6_16_me16=pw_data16{4};
pw6_18_me16=pw_data16{5};
pw6_20_me16=pw_data16{6};
pw6_25_me16=pw_data16{7};
pw6_30_me16=pw_data16{8};

% Lettura pressioni di parete, Me2=1 

pw_data1=leggi_pw('C:\Users\ricca\Documents\MATLAB\Lavori_tesi\TOP_TESI\pw_me1');
pw6_10_me1=pw_data1{1};
pw6_12_me1=pw_data1{2};
pw6_14_me1=pw_data1{3};
pw6_16_me1=pw_data1{4};
pw6_18_me1=pw_data1{5};
pw6_20_me1=pw_data1{6};
pw6_25_me1=pw_data1{7};
pw6_30_me1=pw_data1{8};

% Lettura pressioni di parete, Getto SOTTOESPANSO

pw_data_sot=leggi_pw('C:\Users\ricca\Documents\MATLAB\Lavori_tesi\TOP_TESI\pw6_sot');
pw6_10sot=pw_data_sot{1};
pw6_12sot=pw_data_sot{2};
pw6_14sot=pw_data_sot{3};
pw6_16sot=pw_data_sot{4};
pw6_18sot=pw_data_sot{5};
pw6_20sot=pw_data_sot{6};
pw6_25sot=pw_data_sot{7};
pw6_30sot=pw_data_sot{8};


% Confronto wall pressure arbitrario, 
% lo uso per simulazioni npr=14,20 a getto OFF

% figure
% % plot(pw_EUL(:,1)./L,pw_EUL(:,2)./p0_10,'k--','LineWidth',1.0)
% % hold on
% plot(pw6_20(:,1)./L,pw6_20(:,2)./10^5,'r-','LineWidth',2.0)
% hold on
% yline(pa/10^5,'k-.','LineWidth',1)
% grid on
% title('Pressione di parete','NPR=20')
% xlabel('x / L')
% ylabel('p_{wall} [10^5 Pa]')
% xlim([0.02,1.05])

%% Studio spostamento della separazione

% Posizioni punti di separazione

Pw=[pw6_10(:,2),pw6_12(:,2),pw6_14(:,2),pw6_16(:,2),pw6_18(:,2),...
    pw6_20(:,2),pw6_22(:,2),pw6_25(:,2),pw6_30(:,2),pw6_35(:,2),pw6_40(:,2),...
    pw6_45(:,2),pw6_50(:,2)];
Pw2=[pw6_10_me2(:,2),pw6_12_me2(:,2),pw6_14_me2(:,2),pw6_16_me2(:,2),...
    pw6_18_me2(:,2),pw6_20_me2(:,2),pw6_25_me2(:,2),pw6_30_me2(:,2)];
Pw18=[pw6_10_me18(:,2),pw6_12_me18(:,2),pw6_14_me18(:,2),pw6_16_me18(:,2),...
    pw6_18_me18(:,2),pw6_20_me18(:,2),pw6_25_me18(:,2),pw6_30_me18(:,2)];
Pw16=[pw6_10_me16(:,2),pw6_12_me16(:,2),pw6_14_me16(:,2),pw6_16_me16(:,2),...
    pw6_18_me16(:,2),pw6_20_me16(:,2),pw6_25_me16(:,2),pw6_30_me16(:,2)];
Pw1=[pw6_10_me1(:,2),pw6_12_me1(:,2),pw6_14_me1(:,2),pw6_16_me1(:,2),...
    pw6_18_me1(:,2),pw6_20_me1(:,2),pw6_25_me1(:,2),pw6_30_me1(:,2)];
Pw_sot=[pw6_10sot(:,2),pw6_12sot(:,2),pw6_14sot(:,2),pw6_16sot(:,2),...
    pw6_18sot(:,2),pw6_20sot(:,2),pw6_25sot(:,2),pw6_30sot(:,2)];

X=[pw6_10(:,1)];
[min_off,ind_off]=min(Pw);
[min2,ind2]=min(Pw2);
[min18,ind18]=min(Pw18);
[min16,ind16]=min(Pw16);
[min1,ind1]=min(Pw1);
[min_sot,ind_sot]=min(Pw_sot);
Xsep_off=X(ind_off);
Xsep2=X(ind2);
Xsep18=X(ind18);
Xsep16=X(ind16);
Xsep1=X(ind1);
Xsep_sot=X(ind_sot);

NPR_v=[10 12 14 16 18 20 22 25 30 35 40 45 50]';

% %Fitting di tre tratti, info getto OFF (OLD)
% 
% % help polyfit e polyval, osservo due tratti lineari FSS, RSS e poi 
% % un polinomiale, sarebbe la riapertura della bolla

% %Tratto 1, NPR=10,12,14,16
% % X1=Xsep_off(1:4);
% % xv1=linspace(X1(1),X1(end),100);
% % NPR1=[10,12,14,16];
% % [p1,S1]=polyfit(X1,NPR1,1);
% % [yv1,delta1]=polyval(p1,xv1,S1);
% 
% %Tratto 2, NPR=18,20,22,25
% % X2=Xsep_off(5:8);
% % xv2=linspace(X2(1),X2(end),100);
% % NPR2=[18,20,22,25];
% % [p2,S2]=polyfit(X2,NPR2,1);
% % [yv2,delta2]=polyval(p2,xv2,S2);
% 
% %Tratto3, NPR=30,35,40,45,50
% % X3=Xsep_off(9:13);
% % xv3=linspace(X3(1),X3(end),100);
% % NPR3=[30,35,40,45,50];
% % %[p3,S3]=polyfit(X3,NPR3,3); %cubica
% % [p3,S3]=polyfit(X3,NPR3,2); %parabola
% % [yv3,delta3]=polyval(p3,xv3,S3);
% 
% 
% % figure
% % subplot(2,1,1)
% % plot(Xsep_off./L,NPR_v,'k o','MarkerFaceColor','k','MarkerSize',5)
% % hold on
% % plot(xv1./L,yv1,'b-')
% % plot(xv2./L,yv2,'r-')
% % plot(xv3./L,yv3,'g-')
% % title('\fontname{Times New Roman}Posizione punto di separazione, NPR variabile',...
% % '\fontname{Times New Roman} Getto secondario spento','FontSize',12)
% % legend('NPR CFD ','$NPR=a_1\,(X_{sep}/L)+b_1$',...
% %     '$NPR=a_2\,(X_{sep}/L)+b_2$',...
% %     '$NPR=a_3\,(X_{sep}/L)^3+b_3\,(X_{sep}/L)^2+c_3\,(X_{sep}/L)+d_3$','Interpreter', 'latex')
% % 
% % ylim([5,55])
% % xlim([0,1.1])
% % xlabel('X_{sep} / L')
% % ylabel(' NPR ')
% % grid minor
% % 
% % m1=mean(delta1);
% % m2=mean(delta2);
% % m3=mean(delta3);
% % disp('Errori std medi')
% % M=[m1,m2,m3];
% % 
% % subplot(2,1,2)
% % plot(xv1./L,delta1,'b- ','LineWidth',1.5)
% % hold on
% % plot(xv1./L,m1*ones(1,100),'b-.')
% % 
% % plot(xv2./L,delta2,'r- ','LineWidth',1.5)
% % plot(xv2./L,m2*ones(1,100),'r-.')
% % 
% % plot(xv3./L,delta3,'g- ','LineWidth',1.5)
% % plot(xv3./L,m3*ones(1,100),'g-.')
% % xlim([0,1.1])
% % xlabel('X_{sep} / L')
% % ylabel('MAX |NPR-NPR_{cfd}| ')
% % grid minor
% % title('\fontname{Times New Roman}Errore standard correlazioni NPR = f (X_{sep})',...
% %     'FontSize',12)
% % legend('Err. correlazione FSS','','Err. correlazione RSS','', 'Err. correlazione FSS - end')

% Inversione grafici (Xsep/L - NPR) (NEW)

% %Tratto 1, NPR=10,12,14,16
% NPR1=[10:2:16];
% NPR1v=linspace(NPR1(1),NPR1(end),100);
% X1=Xsep_off(1:4)./L;
% [p1,S1]=polyfit(NPR1,X1,1);
% [yv1,delta1]=polyval(p1,NPR1v,S1);
% 
% %Tratto 2, NPR=18,20,22,25
% NPR2=[18,20,22,25];
% NPR2v=linspace(NPR2(1),NPR2(end),100);
% X2=Xsep_off(5:8)./L;
% [p2,S2]=polyfit(NPR2,X2,1);
% [yv2,delta2]=polyval(p2,NPR2v,S2);
% 
% %Tratto3, NPR=30,35,40,45,50
% NPR3=[30,35,40,45,50];
% NPR3v=linspace(NPR3(1),NPR3(end),100);
% X3=Xsep_off(9:13)./L;
% %[p3,S3]=polyfit(X3,NPR3,3); %cubica
% [p3,S3]=polyfit(NPR3,X3,2); %parabola
% [yv3,delta3]=polyval(p3,NPR3v,S3);

% m1=mean(delta1);
% m2=mean(delta2);
% m3=mean(delta3);
% M=[m1,m2,m3];
% 
% figure
% plot(NPR_v,Xsep_off./L,'k *','MarkerFaceColor','k','MarkerSize',8)
% hold on
% plot(NPR1v,yv1,'b-','LineWidth',1)
% plot(NPR2v,yv2,'r-','LineWidth',1)
% plot(NPR3v,yv3,'g-','LineWidth',1)
% 
% plot(NPR1v,yv1+2*delta1,'m-.')  % 95% probabilità di beccare valore
% plot(NPR1v,yv1-2*delta1,'m-.')
% plot(NPR2v,yv2+2*delta2,'m-.')
% plot(NPR2v,yv2-2*delta2,'m-.')
% plot(NPR3v,yv3+2*delta3,'m-.')
% plot(NPR3v,yv3-2*delta3,'m-.')
% 
% %xline([10,16,18,25,30,50],'k--')
% 
% xlabel('NPR')
% ylabel('X_{sep} / L')
% title('\fontname{Times New Roman}Posizione punto di separazione, NPR variabile',...
% '\fontname{Times New Roman} Getto secondario spento','FontSize',12)
%  legend('NPR CFD ','$X_{sep}/L=f_1(NPR)\;$ [FSS]',...
%      '$X_{sep}/L=f_2(NPR)\;$ [RSS]',...
%      '$X_{sep}/L=f_3(NPR)\;$ [FSS-end]',...
%      '$95\%$ prediction interval','Interpreter','latex')
% xlim([8,52])
% ylim([0,1.1])
% grid minor

% figure
% % plot(NPR1v,delta1,'k-.','LineWidth',1)
% % hold on
% % plot(NPR2v,delta2,'k-.','LineWidth',1)
% % plot(NPR3v,delta3,'k-.','LineWidth',1)
% area(NPR1v,m1*ones(1,100), 'FaceColor', 'b', 'FaceAlpha', 0.3)
% hold on
% area(NPR2v,m2*ones(1,100), 'FaceColor', 'r', 'FaceAlpha', 0.3)
% area(NPR3v,m3*ones(1,100), 'FaceColor', 'g', 'FaceAlpha', 0.3)
% plot(NPR1v,m1*ones(1,100),'b-')
% plot(NPR2v,m2*ones(1,100),'r-')
% plot(NPR3v,m3*ones(1,100),'g-')
% 
% title('Errore standard medio','Correlazioni X_{sep} / L = f (NPR)')
% xlabel('NPR')
% ylabel('$\frac{1}{2}\cdot$max$(|X_{sep}-X_{sep,CFD}|)/L$','Interpreter','latex')
% xlim([8,52])
% ylim([0,15]*10^(-3))
% grid minor
% xline([10,16,18,25,30,50],'k--')

%[y,delta] = polyval(p,x,S) uses the optional output structure S produced by polyfit to generate error estimates. delta 
%is an estimate of the standard error in predicting a future observation at x by p(x).
%Standard error for prediction, returned as a vector of the same size as the query points x. 
%Generally, an interval of y ± Δ corresponds to a roughly 68% prediction interval for future observations of large samples, and y ± 2Δ a roughly 95% prediction interval.
%If the coefficients in p are least-squares estimates computed by polyfit, and the errors in the data input to polyfit are independent, 
%normal, and have constant variance, then y ± Δ is at least a 50% prediction interval.


%% Grafico totale, getto secondario ON

figure
%plot(Xsep_off./L,NPR_v,'k o','MarkerFaceColor','k','MarkerSize',4)
%hold on
%plot(Xsep2./L,[10 12 14 16 18 20 25 30],'r v','MarkerFaceColor','r','MarkerSize',5)
%plot(Xsep18./L,[10 12 14 16 18 20 25 30],'g ^','MarkerFaceColor','g','MarkerSize',5)
plot(Xsep16./L,[10 12 14 16 18 20 25 30],'b +','MarkerFaceColor','b','MarkerSize',5)
%plot(Xsep1./L,[10 12 14 16 18 20 25 30],'m *','MarkerFaceColor','m','MarkerSize',5)
hold on
plot(Xsep_sot./L,[10 12 14 16 18 20 25 30],'r x ','MarkerSize',8)

% yline(18,'k--')
% yline(30,'k--')


%transizioni FSS-RSS
%plot(Xsep_off(4:5)./L,NPR_v(4:5),'k--','LineWidth',0.5)
%plot(Xsep2(1:2)./L,NPR_v(1:2),'r--','LineWidth',0.5)
%plot(Xsep18(2:3)./L,NPR_v(2:3),'g--','LineWidth',0.5)
plot(Xsep16(2:3)./L,NPR_v(2:3),'b-','LineWidth',0.5)
plot(Xsep_sot(2:3)./L,NPR_v(2:3),'r--','LineWidth',0.5)
%plot(Xsep1(3:4)./L,NPR_v(3:4),'m--','LineWidth',0.5)


%transizione RSS-FSS-end
%plot(Xsep_off(8:9)./L,NPR_v(8:9),'k:','LineWidth',0.5)

%transizione RSS-full flow
%plot(Xsep2(end-1:end)./L,[25,30],'k-.','LineWidth',0.5)


 % legend(' Getto2 OFF','$M_e=2$','$M_e=1.8$',...
 %      '$M_e=1.6$','$M_e=1$','FSS$\rightarrow$RSS (Getto2 OFF)','FSS$\rightarrow$RSS ($M_e=2$)','FSS$\rightarrow$RSS ($M_e=1.8$)',...
 %      'FSS$\rightarrow$RSS ($M_e=1.6$)','FSS$\rightarrow$RSS ($M_e=1.0$)',...
 %      'Interpreter', 'latex','Numcolumns',2,...
 %      'Location','best')

xline(1,'k--')

legend('Getto2 adattato ($M_e=2.0$)','Getto2 sottoespanso ($p_e=1.2\cdot p_a$)',...
    '$FSS\rightarrow RSS$ (Getto2 adattato)', '$FSS\rightarrow RSS$ (Getto2 sottoespanso)',...
    'Condizioni full-flowing','Interpreter','latex')

ylim([5,35])
xlim([0,1.1])

title('\fontname{Times New Roman}Effetto del getto secondario sul punto di separazione',...
    '$p_e=p_a=101325$ Pa','FontSize',12)
xlabel('X_{sep} / L')
ylabel(' NPR ')
grid minor

% cftool per aprire interfaccia del Curve Fitting Toolbox

%% Profili velocità assiale, sezione uscita

% Lettura profili velocità assiale U(x=L,y)

Ue_data=leggi_Ue('C:\Users\ricca\Documents\MATLAB\Lavori_tesi\TOP_TESI\Profili_exit');

c=10^2;

Ue14_me1=Ue_data{1};
Ue14_me2=Ue_data{4};
Ue14_me16=Ue_data{2};
Ue14_me18=Ue_data{3};
Ue14_off=Ue_data{5};

Ue30_me1=Ue_data{6};
Ue30_me2=Ue_data{9};
Ue30_me16=Ue_data{7};
Ue30_me18=Ue_data{8};
Ue30_off=Ue_data{10};

figure
plot(Ue14_off(:,1),Ue14_off(:,2)*c,'k-.','LineWidth',1)
hold on
plot(Ue14_me1(:,1),Ue14_me1(:,2)*c,'m-','LineWidth',1)
plot(Ue14_me16(:,1),Ue14_me16(:,2)*c,'b-','LineWidth',1)
plot(Ue14_me18(:,1),Ue14_me18(:,2)*c,'g-','LineWidth',1)
plot(Ue14_me2(:,1),Ue14_me2(:,2)*c,'r-','LineWidth',1)
xline(0,'k--')
legend('Getto OFF','$M_e=1$','$M_e=1.6$',...
    '$M_e=1.8$','$M_e=2.0$','Interpreter','latex')
ylim([0,4.6])
xlim([-120,700])

xlabel('U [m/s]')
ylabel('y [cm]')
grid minor
title('\fontname{Times New Roman}Profilo velocità assiale in sezione d''uscita',...
    '\fontname{Times New Roman}NPR = 14', 'FontSize',12)


figure
plot(Ue30_off(:,1),Ue30_off(:,2)*c,'k-.','LineWidth',1)
hold on
plot(Ue30_me1(:,1),Ue30_me1(:,2)*c,'m-','LineWidth',1)
plot(Ue30_me16(:,1),Ue30_me16(:,2)*c,'b-','LineWidth',1)
plot(Ue30_me18(:,1),Ue30_me18(:,2)*c,'g-','LineWidth',1)
plot(Ue30_me2(:,1),Ue30_me2(:,2)*c,'r-','LineWidth',1)
xline(0,'k--')

xlim([-120,700])
ylim([0,4.6])
title('\fontname{Times New Roman}Profilo velocità assiale in sezione d''uscita',...
    '\fontname{Times New Roman}NPR = 30', 'FontSize',12)
legend('Getto OFF','$M_e=1$','$M_e=1.6$',...
    '$M_e=1.8$','$M_e=2.0$','','Interpreter','latex')
xlabel('U [m/s]')
ylabel('y [cm]')
grid minor
%% Getto 2 sottoexp confronto
%NPR2=4.25
X_diff=Xsep16-Xsep_sot;
figure
plot([10 12 14 16 18 20 25 30],X_diff,'k *')


%% Mach sulla centerline (mesh 6, HSN, ROE 2 ordine)
% GETTO OFF QUI

% M_data=leggi_pw('C:\Users\ricca\Documents\MATLAB\Lavori_tesi\TOP_TESI\Mach_centerline');
% M10=M_data{1};
% M12=M_data{2};
% M14=M_data{3};
% M16=M_data{4};
% M18=M_data{5};
% M20=M_data{6};
% M22=M_data{7};
% M25=M_data{8};
% M30=M_data{9};
% M35=M_data{10};
% M40=M_data{11};
% M45=M_data{12};
% M50=M_data{13};
% 
% % Troviamo posizione disco di mach su centerline con x / ho max(Mach)
% M_mat1=[M10,M12(:,2),M14(:,2),M16(:,2),M18(:,2),...
%     M20(:,2),M22(:,2),M25(:,2),M30(:,2),M35(:,2)];
% M_mat2=[M40,M45(:,2),M50(:,2)];
% [max1,ind1]=max(M_mat1(:,2:end));
% [max2,ind2]=max(M_mat2(:,2:end));
% X1=M10(:,1);
% X2=M40(:,1);
% X_disk=[X1(ind1)',X2(ind2)'];
% NPR_disk=[10:2:22,25:5:50];
% 
% % Distanza assiale fra disco di mach e punto separazione
% dist=(X_disk-Xsep_off')*100;
% 
% % Disco mach position
% p1=polyfit(10:2:16,X_disk(1:4),1);
% p2=polyfit(NPR_disk(5:end),X_disk(5:end),2);
% X_diskFSS=polyval(p1,10:0.01:16);
% X_diskRSS=polyval(p2,18:0.1:50);
% 
% % Distanza MD - X_sep
%  p3=polyfit(10:2:16,dist(1:4),1);
%  p4=polyfit([18,20,22,25],dist(5:8),1);
%  p5=polyfit(30:5:50,dist(9:13),2);
%  dist3=polyval(p3,10:0.01:16);
%  dist4=polyval(p4,18:0.01:25);
%  dist5=polyval(p5,30:0.1:50);


% figure
% plot(NPR_disk,X_disk./L,'k *','MarkerSize',8)
% hold on
% plot(10:0.01:16,X_diskFSS./L,'g-.','LineWidth',1)
% plot(18:0.1:50,X_diskRSS./L,'r-.','LineWidth',1)
% legend('CFD')
% grid on
% title('Spostamento disco di Mach lungo l''asse dell''ugello',...
%     'Getto secondario spento')
% xlim([8,52])
% ylim([0.35,1.25])


% figure
% plot(NPR_disk,dist,'k *','MarkerSize',8)
% hold on
% plot(10:0.01:16,dist3,'g-.','LineWidth',1)
% plot(18:0.01:25,dist4,'r-.','LineWidth',1)
% plot(30:0.1:50,dist5,'b-.','LineWidth',1)
% grid on
% title('Distanza assiale fra disco di Mach e punto di separazione','Getto secondario spento')
% xlabel('NPR')
% ylabel('$X_{MD}-X_{sep}\;$[cm]','Interpreter','latex')
% xlim([8,52])
% ylim([0.5,3.2])

% figure
% plot(M10(:,1)./L,M10(:,2),'m-','LineWidth',1.1)
% hold on
% plot(M12(:,1)./L,M12(:,2),'LineWidth',1.1)
% plot(M14(:,1)./L,M14(:,2),'LineWidth',1.1)
% plot(M16(:,1)./L,M16(:,2),'LineWidth',1.1)
% plot(M18(:,1)./L,M18(:,2),'LineWidth',1.1)
% plot(M20(:,1)./L,M20(:,2),'b-','LineWidth',1.1)
% grid on
% title('Profili di Mach lungo asse ugello','NPR\in[10,20]')
% xlabel('x / L')
% ylabel('M')
% legend('NPR=10','NPR=12','NPR=14','NPR=16',...
%     'NPR=18','NPR=20','Interpreter','latex')

