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.7+Pa;%MPa
R0=0.72e-3;
cG=60000/rho;%da kg/s a l/min 
fs=100;
k=1.4;

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

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

load 'Apertura_7bar_1V_bocca2.txt'
Q2_1V=Apertura_7bar_1V_bocca2(:,1);
Pserb_1V=Apertura_7bar_1V_bocca2(:,3)*0.1+Pa;
load 'Apertura_7bar_15V_bocca2.txt'
Q2_15V=Apertura_7bar_15V_bocca2(:,1);
Pserb_15V=Apertura_7bar_15V_bocca2(:,3)*0.1+Pa;
load 'Apertura_7bar_2V_bocca2.txt'
Q2_2V=Apertura_7bar_2V_bocca2(:,1);
Pserb_2V=Apertura_7bar_2V_bocca2(:,3)*0.1+Pa;
load 'Apertura_7bar_25V_bocca2.txt'
Q2_25V=Apertura_7bar_25V_bocca2(:,1);
Pserb_25V=Apertura_7bar_25V_bocca2(:,3)*0.1+Pa;
load 'Apertura_7bar_3V_bocca2.txt'
Q2_3V=Apertura_7bar_3V_bocca2(:,1);
Pserb_3V=Apertura_7bar_3V_bocca2(:,3)*0.1+Pa;
load 'Apertura_7bar_35V_bocca2.txt'
Q2_35V=Apertura_7bar_35V_bocca2(:,1);
Pserb_35V=Apertura_7bar_35V_bocca2(:,3)*0.1+Pa;
load 'Apertura_7bar_4V_bocca2.txt'
Q2_4V=Apertura_7bar_4V_bocca2(:,1);
Pserb_4V=Apertura_7bar_4V_bocca2(:,3)*0.1+Pa;
load 'Apertura_7bar_45V_bocca2.txt'
Q2_45V=Apertura_7bar_45V_bocca2(:,1);
Pserb_45V=Apertura_7bar_45V_bocca2(:,3)*0.1+Pa;
timeadd=(0:size(Q2_45V,1)-1)*1/fs;

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


massimi=[max(abs(Q2_0V)),max(abs(Q2_05V)),max(abs(Q2_1V)),max(abs(Q2_15V)),max(abs(Q2_2V)),max(abs(Q2_25V)),max(abs(Q2_3V)),max(abs(Q2_35V)),max(abs(Q2_4V)),max(abs(Q2_45V)),max(abs(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

for i=2:11
    a(1,1)=1;
    a(i,1)=A(i,1)/A(1,1);
end


% Sovrapposizione Aree

figure
plot((time)*a(1,1),Pserb_0V,'b','MarkerSize',12,'linewidth',4),hold on
plot((time)*a(2,1)-0.538,Pserb_05V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(3,1)+0.329,Pserb_1V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(4,1)+0.704,Pserb_15V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(5,1)+0.891,Pserb_2V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*a(6,1)+1.51,Pserb_25V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*a(7,1)+2.178,Pserb_3V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*a(8,1)+2.4743,Pserb_35V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*a(9,1)+2.8949,Pserb_4V,':','MarkerSize',12,'linewidth',3),hold on 
plot((timeadd)*a(10,1)+3.1391,Pserb_45V,':','MarkerSize',12,'linewidth',3),hold on 
% % plot((timeadd)*a(11,1),Pserb_5V,'b','MarkerSize',12,'linewidth',2),hold on

ylabel('P [$MPa$]','interpreter','latex'),xlabel('t [$s$]','interpreter','latex')
% xlim([0 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 ALBENDER
figure
plot((time)*a(1,1),Pserb_0V,'b','MarkerSize',12,'linewidth',4),hold on
plot((time)*a(2,1)-0.538,Pserb_05V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(3,1)+0.329,Pserb_1V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(4,1)+0.704,Pserb_15V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.65+1.032,Pserb_2V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.51+1.657,Pserb_25V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.373+2.272,Pserb_3V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.245+2.5629,Pserb_35V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.126+2.9521,Pserb_4V,':','MarkerSize',12,'linewidth',3),hold on 
plot((timeadd)*0.04+3.2052,Pserb_45V,':','MarkerSize',12,'linewidth',3),hold on 
% % plot((timeadd)*a(11,1),Pserb_5V,'b','MarkerSize',12,'linewidth',2),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')

%%

V=linspace(0,5,11);
A_man=[A(1,1)*a(1,1);A(1,1)*a(2,1);A(1,1)*a(3,1);A(1,1)*a(4,1);A(1,1)*0.65;A(1,1)*0.51;A(1,1)*0.373;A(1,1)*0.245;A(1,1)*0.126;A(1,1)*0.04;0];

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

figure 
plot(V,A_man,'bo','MarkerSize',12,'linewidth',3), hold on 
plot(myfit1,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')
