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


load 'Scarico_4bar_55V_bocca4.txt'
Q4_55V=Scarico_4bar_55V_bocca4(:,2);
Pserb_55V=Scarico_4bar_55V_bocca4(:,3)*0.1+Pa;
load 'Scarico_4bar_6V_bocca4.txt'
Q4_6V=Scarico_4bar_6V_bocca4(:,2);
Pserb_6V=Scarico_4bar_6V_bocca4(:,3)*0.1+Pa;

timeadd=(0:size(Q4_55V,1)-1)*1/fs; 
time=(0:size(Q4_6V,1)-1)*1/fs;


load 'Scarico_4bar_65V_bocca4.txt'
Q4_65V=Scarico_4bar_65V_bocca4(:,2);
Pserb_65V=Scarico_4bar_65V_bocca4(:,3)*0.1+Pa;
load 'Scarico_4bar_7V_bocca4.txt'
Q4_7V=Scarico_4bar_7V_bocca4(:,2);
Pserb_7V=Scarico_4bar_7V_bocca4(:,3)*0.1+Pa;
load 'Scarico_4bar_75V_bocca4.txt'
Q4_75V=Scarico_4bar_75V_bocca4(:,2);
Pserb_75V=Scarico_4bar_75V_bocca4(:,3)*0.1+Pa;
load 'Scarico_4bar_8V_bocca4.txt'
Q4_8V=Scarico_4bar_8V_bocca4(:,2);
Pserb_8V=Scarico_4bar_8V_bocca4(:,3)*0.1+Pa;
load 'Scarico_4bar_85V_bocca4.txt'
Q4_85V=Scarico_4bar_85V_bocca4(:,2);
Pserb_85V=Scarico_4bar_85V_bocca4(:,3)*0.1+Pa;
load 'Scarico_4bar_9V_bocca4.txt'
Q4_9V=Scarico_4bar_9V_bocca4(:,2);
Pserb_9V=Scarico_4bar_9V_bocca4(:,3)*0.1+Pa;
load 'Scarico_4bar_95V_bocca4.txt'
Q4_95V=Scarico_4bar_95V_bocca4(:,2);
Pserb_95V=Scarico_4bar_95V_bocca4(:,3)*0.1+Pa;
load 'Scarico_4bar_10V_bocca4.txt'
Q4_10V=Scarico_4bar_10V_bocca4(:,2);
Pserb_10V=Scarico_4bar_10V_bocca4(:,3)*0.1+Pa;

% load 'Scarico_4bar_5V_bocca4.txt'
% Q4_5V=Scarico_4bar_5V_bocca4(:,1);
% Pserb_5V=Scarico_4bar_5V_bocca4(:,3)*0.1+Pa;

%% Svuotamento nel tempo
figure
plot((timeadd),Pserb_55V,'b','MarkerSize',12,'linewidth',4),hold on
plot((time),Pserb_6V,':','MarkerSize',12,'linewidth',3),hold on
plot((time),Pserb_65V,':','MarkerSize',12,'linewidth',3),hold on
plot((time),Pserb_7V,':','MarkerSize',12,'linewidth',3),hold on
plot((time),Pserb_75V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time),Pserb_8V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time),Pserb_85V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time),Pserb_9V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time),Pserb_95V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time),Pserb_10V,':','MarkerSize',12,'linewidth',3),hold on 


ylabel('[P] MPa'),xlabel('[T] s','interpreter','latex')
xlim([0 20])
grid on,zoom on,title('Curve di svuotamento')
set(gca,'FontName','Times New Roman'),set(gca,'FontSize',32)  
legend('$ 5.5 V $ ','$ 6.0 V $ ','$ 6.5 V $ ','$ 7.0 V $ ','$ 7.5 V $ ','$ 8.0 V $ ','$ 8.5 V $ ','$ 9.0 V $ ','$ 9.5 V $ ','$ 10.0 V $','FontSize',32,'interpreter','latex')


%%
massimi=[max(abs(Q4_55V)),max(abs(Q4_6V)),max(abs(Q4_65V)),max(abs(Q4_7V)),max(abs(Q4_75V)),max(abs(Q4_8V)),max(abs(Q4_85V)),max(abs(Q4_9V)),max(abs(Q4_95V)),max(abs(Q4_10V))]';
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=1:10
    a(i,1)=A(i,1)/A(10,1);
end

%% Sovrapposto con rapporto aree
figure
plot((timeadd)*a(1,1)+2.8459,Pserb_55V+0.02,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(2,1)+2.604,Pserb_6V+0.003,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(3,1)+2.1998,Pserb_65V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(4,1)+1.8,Pserb_7V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(5,1)+1.632,Pserb_75V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(6,1)+1.075,Pserb_8V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(7,1)+0.801,Pserb_85V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(8,1)+0.549,Pserb_9V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(9,1)+0.273,Pserb_95V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(10,1),Pserb_10V,'m','MarkerSize',12,'linewidth',4),hold on %
% plot((timeadd),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 svuotamento')
set(gca,'FontName','Times New Roman'),set(gca,'FontSize',25)  
legend('$ 5.5 V $ ','$ 6.0 V $ ','$ 6.5 V $ ','$ 7.0 V $ ','$ 7.5 V $ ','$ 8.0 V $ ','$ 8.5 V $ ','$ 9.0 V $ ','$ 9.5 V $ ','$ 10.0 V $','FontSize',32,'interpreter','latex')



V=linspace(5.5,10,10);

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

figure
plot(V,A,'bo','MarkerSize',12,'linewidth',3), hold on 
plot(myfit,V,A,'b')
ylabel('[$A_{4}$] $m^2$','interpreter','latex'), xlabel(' V [$V$] ','interpreter','latex')
grid on, zoom on,title('Sezione passaggio uscita 4')
set(gca,'FontName','Times New Roman'),set(gca,'FontSize',32)

%% Svuotamento sovrapposto a mano ALBENDER
figure
plot((timeadd)*0.042+2.91256,Pserb_55V+0.02,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.118+2.7162,Pserb_6V+0.003,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.227+2.3676,Pserb_65V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.35+1.988,Pserb_7V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.5+1.745,Pserb_75V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*0.628+1.189,Pserb_8V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*0.756+0.893,Pserb_85V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*0.84+0.616,Pserb_9V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(9,1)+0.273,Pserb_95V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(10,1),Pserb_10V,'m','MarkerSize',12,'linewidth',4),hold on %
% plot((timeadd),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 svuotamento ALBENDER')
set(gca,'FontName','Times New Roman'),set(gca,'FontSize',25)  
legend('$ 5.5 V $ ','$ 6.0 V $ ','$ 6.5 V $ ','$ 7.0 V $ ','$ 7.5 V $ ','$ 8.0 V $ ','$ 8.5 V $ ','$ 9.0 V $ ','$ 9.5 V $ ','$ 10.0 V $','FontSize',32,'interpreter','latex')

A_man=[A(10,1)*0.042;A(10,1)*0.118;A(10,1)*0.227;A(10,1)*0.35;A(10,1)*0.5;A(10,1)*0.628;A(10,1)*0.756;A(10,1)*0.84;A(10,1)*a(9,1);A(10,1)*a(10,1)];

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

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



%% Confronto Qth/QISO

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

figure 
plot(Pserb_55V./Pa,-G_ISO_55V.*cG,'k','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_55V./Pa,Q4_55V,'ko','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_6V./Pa,-G_ISO_6V.*cG,'r','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_6V./Pa,Q4_6V,'ro','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_7V./Pa,-G_ISO_7V.*cG,'b','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_7V./Pa,Q4_7V,'bo','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_75V./Pa,-G_ISO_75V.*cG,'y','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_75V./Pa,Q4_75V,'yo','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_8V./Pa,-G_ISO_8V.*cG,'c','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_8V./Pa,Q4_8V,'co','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_85V./Pa,-G_ISO_85V.*cG,'g','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_85V./Pa,Q4_85V,'go','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_9V./Pa,-G_ISO_9V.*cG,'k','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_9V./Pa,Q4_9V,'ko','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_95V./Pa,-G_ISO_95V.*cG,'m','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_95V./Pa,Q4_95V,'mo','MarkerSize',5,'linewidth',1.3),hold on
plot(Pserb_10V./Pa,-G_ISO_10V.*cG,'r','MarkerSize',5,'linewidth',2),hold on
plot(Pserb_10V./Pa,Q4_10V,'ro','MarkerSize',3,'linewidth',1.3),hold on


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