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


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

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

load 'Scarico_2bar_1V_bocca2.txt'
Q2_1V=Scarico_2bar_1V_bocca2(:,1);
Pserb_1V=Scarico_2bar_1V_bocca2(:,3)*0.1+Pa;
load 'Scarico_2bar_15V_bocca2.txt'
Q2_15V=Scarico_2bar_15V_bocca2(:,1);
Pserb_15V=Scarico_2bar_15V_bocca2(:,3)*0.1+Pa;
load 'Scarico_2bar_2V_bocca2.txt'
Q2_2V=Scarico_2bar_2V_bocca2(:,1);
Pserb_2V=Scarico_2bar_2V_bocca2(:,3)*0.1+Pa;
load 'Scarico_2bar_25V_bocca2.txt'
Q2_25V=Scarico_2bar_25V_bocca2(:,1);
Pserb_25V=Scarico_2bar_25V_bocca2(:,3)*0.1+Pa;
load 'Scarico_2bar_3V_bocca2.txt'
Q2_3V=Scarico_2bar_3V_bocca2(:,1);
Pserb_3V=Scarico_2bar_3V_bocca2(:,3)*0.1+Pa;
load 'Scarico_2bar_35V_bocca2.txt'
Q2_35V=Scarico_2bar_35V_bocca2(:,1);
Pserb_35V=Scarico_2bar_35V_bocca2(:,3)*0.1+Pa;
load 'Scarico_2bar_4V_bocca2.txt'
Q2_4V=Scarico_2bar_4V_bocca2(:,1);
Pserb_4V=Scarico_2bar_4V_bocca2(:,3)*0.1+Pa;
load 'Scarico_2bar_45V_bocca2.txt'
Q2_45V=Scarico_2bar_45V_bocca2(:,1);
Pserb_45V=Scarico_2bar_45V_bocca2(:,3)*0.1+Pa;
timeadd=(0:size(Q2_45V,1)-1)*1/fs;

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

%% Svuotamento nel tempo
figure
plot((time),Pserb_0V,'b','MarkerSize',12,'linewidth',4),hold on
plot((time),Pserb_05V,':','MarkerSize',12,'linewidth',3),hold on
plot((time),Pserb_1V,':','MarkerSize',12,'linewidth',3),hold on
plot((time),Pserb_15V,':','MarkerSize',12,'linewidth',3),hold on
plot((time),Pserb_2V,':','MarkerSize',12,'linewidth',3),hold on %*0.712+1.777
plot((time),Pserb_25V,':','MarkerSize',12,'linewidth',3),hold on %0.567+2.1945
plot((time),Pserb_3V,':','MarkerSize',12,'linewidth',3),hold on %0.43+2.588
plot((time),Pserb_35V,':','MarkerSize',12,'linewidth',3),hold on %0.294+2.6
plot((time),Pserb_4V,':','MarkerSize',12,'linewidth',3),hold on %0.166+2.87
plot((timeadd),Pserb_45V,':','MarkerSize',12,'linewidth',3),hold on %0.055+3.04
plot((timeadd),Pserb_5V,'b','MarkerSize',12,'linewidth',2),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('$ 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 $','$ 5 V $','FontSize',32,'interpreter','latex')


%%
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:10
    a(1,1)=1;
    a(i,1)=A(i,1)/A(1,1);
end


%% Sovrapposto con rapporto aree
figure
plot((time)*a(1,1),Pserb_0V,'b','MarkerSize',12,'linewidth',4),hold on
plot((time)*a(2,1)+0.166,Pserb_05V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(3,1)-0.136,Pserb_1V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(4,1)-1.109,Pserb_15V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(5,1)-0.039,Pserb_2V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(6,1)+0.3153,Pserb_25V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(7,1)+0.5897,Pserb_3V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(8,1)+0.8943,Pserb_35V,':','MarkerSize',12,'linewidth',3),hold on %
plot((time)*a(9,1)+0.9716,Pserb_4V,':','MarkerSize',12,'linewidth',3),hold on %
plot((timeadd)*a(10,1)+1.15158,Pserb_45V,':','MarkerSize',12,'linewidth',3),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('$ 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);

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_{2}$] $m^2$','interpreter','latex'), xlabel(' V ')
grid on, zoom on,title('Sezione passaggio uscita 2')
set(gca,'FontName','Times New Roman'),set(gca,'FontSize',32)

%% Svuotamento sovrapposto a mano ALBENDER
figure
plot((time)*a(1,1),Pserb_0V,'b','MarkerSize',12,'linewidth',4),hold on
plot((time)*a(2,1)+0.166,Pserb_05V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*a(3,1)-0.136,Pserb_1V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.77-0.999,Pserb_15V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.67+0.007,Pserb_2V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.533+0.3592,Pserb_25V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.385+0.6448,Pserb_3V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.242+0.9477,Pserb_35V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.12+1.0348,Pserb_4V,':','MarkerSize',12,'linewidth',3),hold on
plot((timeadd)*0.04+1.18,Pserb_45V,':','MarkerSize',12,'linewidth',3),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('$ 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')

A_man=[A(1,1);A(2,1);A(3,1);A(1,1)*0.77;A(1,1)*0.67;A(1,1)*0.533;A(1,1)*0.385;A(1,1)*0.242;A(1,1)*0.12;A(1,1)*0.04;A(11,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_{2}$ [$m^2$]','interpreter','latex'), xlabel(' V [$V$] ','interpreter','latex')
grid on, zoom on  %title,('Sezione passaggio uscita 2 manulae')
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)/Pa<b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_0V(i)*1e6*A_eff*sqrt(1-((Pserb_0V(i)/Pa-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_0V(i)*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)/Pa<b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_05V(i)*1e6*A_eff*sqrt(1-((Pserb_05V(i)/Pa-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_05V(i)*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)/Pa<b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_1V(i)*1e6*A_eff*sqrt(1-((Pserb_1V(i)/Pa-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_1V(i)*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)/Pa<b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_15V(i)*1e6*A_eff*sqrt(1-((Pserb_15V(i)/Pa-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_15V(i)*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)/Pa<b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_2V(i)*1e6*A_eff*sqrt(1-((Pserb_2V(i)/Pa-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_2V(i)*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)/Pa<b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_25V(i)*1e6*A_eff*sqrt(1-((Pserb_25V(i)/Pa-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_25V(i)*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)/Pa<b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_3V(i)*1e6*A_eff*sqrt(1-((Pserb_3V(i)/Pa-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_3V(i)*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)/Pa<b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_35V(i)*1e6*A_eff*sqrt(1-((Pserb_35V(i)/Pa-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_35V(i)*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)/Pa<b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_4V(i)*1e6*A_eff*sqrt(1-((Pserb_4V(i)/Pa-b)/(1-b))^2); %kg/s
else
%Sonico

g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_4V(i)*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)/Pa<b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Pserb_45V(i)*1e6*A_eff*sqrt(1-((Pserb_45V(i)/Pa-b)/(1-b))^2); %kg/s
else
%Sonico

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

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

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