close all
clear all
clc

Ts=293;              %Temperatura ambiente durante la prova [K]
k=1.4;               % Rapporto calori specifici
Rg=287.053;          %Costante dell'aria [J/kg/K]
Pa=0.101325;          %Pressione ambiente durante la prova
rho=Pa*1e6/(Rg*Ts);
cd_ISO=0.87;
b=0.46;
Ps=0.2+Pa;%MPa
R0=0.72e-3;
cG=60000/rho;%da kg/s a l/min 
fs=100;
k=1.4;

load 'Apertura_2bar_0V_bocca2.txt'
Q2_0V=Apertura_2bar_0V_bocca2(:,1);
Pserb_0V=Apertura_2bar_0V_bocca2(:,3)*0.1+Pa;
load 'Apertura_2bar_05V_bocca2.txt'
Q2_05V=Apertura_2bar_05V_bocca2(:,1);
Pserb_05V=Apertura_2bar_05V_bocca2(:,3)*0.1+Pa;

time=(0:size(Q2_05V,1)-1)*1/fs;

load 'Apertura_2bar_1V_bocca2.txt'
Q2_1V=Apertura_2bar_1V_bocca2(:,1);
Pserb_1V=Apertura_2bar_1V_bocca2(:,3)*0.1+Pa;
load 'Apertura_2bar_15V_bocca2.txt'
Q2_15V=Apertura_2bar_15V_bocca2(:,1);
Pserb_15V=Apertura_2bar_15V_bocca2(:,3)*0.1+Pa;
load 'Apertura_2bar_2V_bocca2.txt'
Q2_2V=Apertura_2bar_2V_bocca2(:,1);
Pserb_2V=Apertura_2bar_2V_bocca2(:,3)*0.1+Pa;
load 'Apertura_2bar_25V_bocca2.txt'
Q2_25V=Apertura_2bar_25V_bocca2(:,1);
Pserb_25V=Apertura_2bar_25V_bocca2(:,3)*0.1+Pa;
load 'Apertura_2bar_3V_bocca2.txt'
Q2_3V=Apertura_2bar_3V_bocca2(:,1);
Pserb_3V=Apertura_2bar_3V_bocca2(:,3)*0.1+Pa;
load 'Apertura_2bar_35V_bocca2.txt'
Q2_35V=Apertura_2bar_35V_bocca2(:,1);
Pserb_35V=Apertura_2bar_35V_bocca2(:,3)*0.1+Pa;
load 'Apertura_2bar_4V_bocca2.txt'
Q2_4V=Apertura_2bar_4V_bocca2(:,1);
Pserb_4V=Apertura_2bar_4V_bocca2(:,3)*0.1+Pa;
load 'Apertura_2bar_45V_bocca2.txt'
Q2_45V=Apertura_2bar_45V_bocca2(:,1);
Pserb_45V=Apertura_2bar_45V_bocca2(:,3)*0.1+Pa;
timeadd=(0:size(Q2_45V,1)-1)*1/fs;

load 'Apertura_2bar_5V_bocca2.txt'
Q2_5V=Apertura_2bar_5V_bocca2(:,1);
Pserb_5V=Apertura_2bar_5V_bocca2(:,3)*0.1+Pa;

figure
plot((time),Pserb_0V,'b','MarkerSize',12,'linewidth',4),hold on
plot((time)*0.954+1.155,Pserb_05V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.8833+1.531,Pserb_1V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.810+1.955,Pserb_15V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.712+1.777,Pserb_2V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.567+2.1945,Pserb_25V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.43+2.588,Pserb_3V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.294+2.6,Pserb_35V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.166+2.87,Pserb_4V,':','MarkerSize',12,'linewidth',3),hold on
plot((timeadd)*0.055+3.04,Pserb_45V,':','MarkerSize',12,'linewidth',3),hold on 

ylabel('P [$MPa$]','interpreter','latex'),xlabel('t [$s$]','interpreter','latex')
xlim([2 6])
grid on,zoom on,  %title('Curve di riempimento')
set(gca,'FontName','Times New Roman'),set(gca,'FontSize',25)  
legend('$ 0 V $ ','$ 0.5 V $ ','$ 1.0 V $ ','$ 1.5 V $ ','$ 2.0 V $ ','$ 2.5 V $ ','$ 3.0 V $ ','$ 3.5 V $ ','$ 4.0 V $ ','$ 4.5 V $','FontSize',20,'interpreter','latex')






%% Sovrapposizione manuale

figure
plot((time),Pserb_0V,'b','MarkerSize',12,'linewidth',4),hold on
plot((time)*0.954+1.155,Pserb_05V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.8833+1.531,Pserb_1V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.810+1.955,Pserb_15V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.65+1.895,Pserb_2V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.54+2.239,Pserb_25V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.4+2.626,Pserb_3V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.25+2.68,Pserb_35V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.124+2.937,Pserb_4V,':','MarkerSize',12,'linewidth',3),hold on 
plot((timeadd)*0.037+3.071,Pserb_45V,':','MarkerSize',12,'linewidth',3),hold on 


ylabel('P [$MPa$]','interpreter','latex'),xlabel('t [$s$]','interpreter','latex')
xlim([2 6])
grid on,zoom on,  %title('Curve di riempimento')
set(gca,'FontName','Times New Roman'),set(gca,'FontSize',25)  
legend('$ 0 V $ ','$ 0.5 V $ ','$ 1.0 V $ ','$ 1.5 V $ ','$ 2.0 V $ ','$ 2.5 V $ ','$ 3.0 V $ ','$ 3.5 V $ ','$ 4.0 V $ ','$ 4.5 V $','FontSize',20,'interpreter','latex')

%%
massimi=[max(Q2_0V),max(Q2_05V),max(Q2_1V),max(Q2_15V),max(Q2_2V),max(Q2_25V),max(Q2_3V),max(Q2_35V),max(Q2_4V),max(Q2_45V),max(Q2_5V)]';
for i=1:length(massimi)
A(i,1)=(massimi(i,1)/cG)/(cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6); %m^2
end

V=linspace(0,5,11);
A_man=[A(1,1);A(1,1)*0.954;A(1,1)*0.8833;A(1,1)*0.810;A(1,1)*0.65;A(1,1)*0.54;A(1,1)*0.4;A(1,1)*0.25;A(1,1)*0.124;A(1,1)*0.037;A(11,1)];


myfittype = fittype('m*V+q',...
    'dependent',{'A_man'},'independent',{'V'},...
    'coefficients',{'m','q'})
myfit = fit(V',A_man,myfittype)

figure 
plot(V,A_man,'bo','MarkerSize',12,'linewidth',3), hold on 
plot(myfit,V,A_man,'b')
ylabel('$A_{2}$ [$m^2$]','interpreter','latex'), xlabel(' V [$V$] ','interpreter','latex')
grid on, zoom on,  %title('Sezione passaggio uscita 2 man')
set(gca,'FontName','Times New Roman'),set(gca,'FontSize',25)
%% Confronto Qth/QISO

%V=0.0
for i=1:length(Pserb_0V)
    A_eff=A(1,1);
%Subsonico
if Pserb_0V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_0V(i)/Ps-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff;
end

G_ISOt(i)=g_ISO;
G_ISO_0V=G_ISOt';
end

%V=0.5
for i=1:length(Pserb_05V)
    A_eff=A(2,1);
%Subsonico
if Pserb_05V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_05V(i)/Ps-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff;
end

G_ISOt(i)=g_ISO;
G_ISO_05V=G_ISOt';
end


%V=1
for i=1:length(Pserb_1V)
    A_eff=A(3,1);
%Subsonico
if Pserb_1V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_1V(i)/Ps-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff;
end

G_ISOt(i)=g_ISO;
G_ISO_1V=G_ISOt';
end

%V=1.5
for i=1:length(Pserb_15V)
    A_eff=A(4,1);
%Subsonico
if Pserb_15V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_15V(i)/Ps-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff;
end

G_ISOt(i)=g_ISO;
G_ISO_15V=G_ISOt';
end

%V=2.0
for i=1:length(Pserb_2V)
    A_eff=A(5,1);
%Subsonico
if Pserb_2V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_2V(i)/Ps-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff;
end

G_ISOt(i)=g_ISO;
G_ISO_2V=G_ISOt';
end

%V=2.5
for i=1:length(Pserb_25V)
    A_eff=A(6,1);
%Subsonico
if Pserb_25V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_25V(i)/Ps-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff;
end

G_ISOt(i)=g_ISO;
G_ISO_25V=G_ISOt';
end

%V=3.0
for i=1:length(Pserb_3V)
    A_eff=A(7,1);
%Subsonico
if Pserb_3V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_3V(i)/Ps-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff;
end

G_ISOt(i)=g_ISO;
G_ISO_3V=G_ISOt';
end

%V=3.5
for i=1:length(Pserb_35V)
    A_eff=A(8,1);
%Subsonico
if Pserb_35V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_35V(i)/Ps-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff;
end

G_ISOt(i)=g_ISO;
G_ISO_35V=G_ISOt';
end

%V=4
for i=1:length(Pserb_4V)
    A_eff=A(9,1);
%Subsonico
if Pserb_4V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_4V(i)/Ps-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff;
end

G_ISOt(i)=g_ISO;
G_ISO_4V=G_ISOt';
end

%V=4.5
for i=1:length(Pserb_45V)
    A_eff=A(10,1);
%Subsonico
if Pserb_45V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_45V(i)/Ps-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff;
end

G_ISOt(i)=g_ISO;
G_ISO_45V=G_ISOt';
end

figure
plot(Pserb_0V./Ps,G_ISO_0V.*cG,'k','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_0V./Ps,Q2_0V,'ko','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_05V./Ps,G_ISO_05V.*cG,'r','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_05V./Ps,Q2_05V,'ro','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_15V./Ps,G_ISO_15V.*cG,'b','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_15V./Ps,Q2_15V,'bo','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_2V./Ps,G_ISO_2V.*cG,'y','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_2V./Ps,Q2_2V,'yo','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_25V./Ps,G_ISO_25V.*cG,'c','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_25V./Ps,Q2_25V,'co','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_3V./Ps,G_ISO_3V.*cG,'g','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_3V./Ps,Q2_3V,'go','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_35V./Ps,G_ISO_35V.*cG,'k','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_35V./Ps,Q2_35V,'ko','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_4V./Ps,G_ISO_4V.*cG,'m','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_4V./Ps,Q2_4V,'mo','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_45V./Ps,G_ISO_45V.*cG,'r','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_45V./Ps,Q2_45V,'ro','MarkerSize',3,'linewidth',1.3),hold on

ylabel('Q [$l/min$]','interpreter','latex'),xlabel('$P_{2}/P_{1}$ ','interpreter','latex')
grid on,zoom on, %title('Curve portata')
set(gca,'FontName','Times New Roman'),set(gca,'FontSize',25)  
legend('ISO ','Sperimetale','FontSize',20,'interpreter','latex')
