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

load 'Apertura_4bar_55V_bocca4.txt'
Q4_55V=Apertura_4bar_55V_bocca4(:,2);
Pserb_55V=Apertura_4bar_55V_bocca4(:,3)*0.1+Pa;
load 'Apertura_4bar_6V_bocca4.txt'
Q4_6V=Apertura_4bar_6V_bocca4(:,2);
Pserb_6V=Apertura_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 'Apertura_4bar_65V_bocca4.txt'
Q4_65V=Apertura_4bar_65V_bocca4(:,2);
Pserb_65V=Apertura_4bar_65V_bocca4(:,3)*0.1+Pa;
load 'Apertura_4bar_7V_bocca4.txt'
Q4_7V=Apertura_4bar_7V_bocca4(:,2);
Pserb_7V=Apertura_4bar_7V_bocca4(:,3)*0.1+Pa;
load 'Apertura_4bar_75V_bocca4.txt'
Q4_75V=Apertura_4bar_75V_bocca4(:,2);
Pserb_75V=Apertura_4bar_75V_bocca4(:,3)*0.1+Pa;
load 'Apertura_4bar_8V_bocca4.txt'
Q4_8V=Apertura_4bar_8V_bocca4(:,2);
Pserb_8V=Apertura_4bar_8V_bocca4(:,3)*0.1+Pa;
load 'Apertura_4bar_85V_bocca4.txt'
Q4_85V=Apertura_4bar_85V_bocca4(:,2);
Pserb_85V=Apertura_4bar_85V_bocca4(:,3)*0.1+Pa;
load 'Apertura_4bar_9V_bocca4.txt'
Q4_9V=Apertura_4bar_9V_bocca4(:,2);
Pserb_9V=Apertura_4bar_9V_bocca4(:,3)*0.1+Pa;
load 'Apertura_4bar_95V_bocca4.txt'
Q4_95V=Apertura_4bar_95V_bocca4(:,2);
Pserb_95V=Apertura_4bar_95V_bocca4(:,3)*0.1+Pa;
load 'Apertura_4bar_10V_bocca4.txt'
Q4_10V=Apertura_4bar_10V_bocca4(:,2);
Pserb_10V=Apertura_4bar_10V_bocca4(:,3)*0.1+Pa;



massimi=[max(Q4_55V),max(Q4_6V),max(Q4_65V),max(Q4_7V),max(Q4_75V),max(Q4_8V),max(Q4_85V),max(Q4_9V),max(Q4_95V),max(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
R(i,1)=sqrt(A(i,1)/pi); % m
end
A_max=max(A);
for i=1:length(massimi)
Rapp_A(i,1)=A(i,1)/A_max;
end




figure
plot((time)*Rapp_A(10,1),Pserb_10V,'b','MarkerSize',12,'linewidth',4),hold on
plot((timeadd)*Rapp_A(1,1)+3.4051,Pserb_55V,':','MarkerSize',12,'linewidth',4),hold on
plot((time)*Rapp_A(2,1)+3.1147,Pserb_6V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*Rapp_A(3,1)+2.8155,Pserb_65V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*Rapp_A(4,1)+2.397,Pserb_7V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*Rapp_A(5,1)+1.967,Pserb_75V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*Rapp_A(6,1)+1.551,Pserb_8V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*Rapp_A(7,1)+1.113,Pserb_85V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*Rapp_A(8,1)+0.99,Pserb_9V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*Rapp_A(9,1)+0.509,Pserb_95V,':','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('$ 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',20,'interpreter','latex')






%%Sovrapp coeff ALBENDER

figure
plot((time)*Rapp_A(10,1),Pserb_10V,'b','MarkerSize',12,'linewidth',4),hold on
plot((timeadd)*0.04+3.4492,Pserb_55V,':','MarkerSize',12,'linewidth',4),hold on
plot((time)*0.122+3.1696,Pserb_6V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.235+2.8762,Pserb_65V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*0.359+2.451,Pserb_7V,':','MarkerSize',12,'linewidth',3),hold on
plot((time)*Rapp_A(5,1)+1.967,Pserb_75V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*0.622+1.638,Pserb_8V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*Rapp_A(7,1)+1.113,Pserb_85V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*Rapp_A(8,1)+0.99,Pserb_9V,':','MarkerSize',12,'linewidth',3),hold on 
plot((time)*Rapp_A(9,1)+0.509,Pserb_95V,':','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')



V=linspace(5.5,10,10);
A_man=[A_max*0.04;A_max*0.122;A_max*0.235;A_max*0.359;A_max*Rapp_A(5,1);A_max*0.622;A_max*Rapp_A(7,1);A_max*Rapp_A(8,1);A_max*Rapp_A(9,1);A_max*Rapp_A(10,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_{4}$ [$m^2$]','interpreter','latex'), xlabel(' V [$V$] ','interpreter','latex')
grid on, zoom on,  title('Sezione passaggio uscita 4 man')
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)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_6V(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_6V=G_ISOt';
end

%V=6.5
for i=1:length(Pserb_65V)
    A_eff=A(3,1);
%Subsonico
if Pserb_65V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_65V(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_65V=G_ISOt';
end

%V=7
for i=1:length(Pserb_7V)
    A_eff=A(4,1);
%Subsonico
if Pserb_7V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_7V(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_7V=G_ISOt';
end

%V=7.5
for i=1:length(Pserb_75V)
    A_eff=A(5,1);
%Subsonico
if Pserb_75V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_75V(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_75V=G_ISOt';
end

%V=8
for i=1:length(Pserb_8V)
    A_eff=A(6,1);
%Subsonico
if Pserb_8V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_8V(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_8V=G_ISOt';
end

%V=8.5
for i=1:length(Pserb_85V)
    A_eff=A(7,1);
%Subsonico
if Pserb_85V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_85V(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_85V=G_ISOt';
end

%V=9
for i=1:length(Pserb_9V)
    A_eff=A(8,1);
%Subsonico
if Pserb_9V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_9V(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_9V=G_ISOt';
end

%V=9.5
for i=1:length(Pserb_95V)
    A_eff=A(9,1);
%Subsonico
if Pserb_95V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_95V(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_95V=G_ISOt';
end

%V=10
for i=1:length(Pserb_10V)
    A_eff=A(10,1);
%Subsonico
if Pserb_10V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_10V(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_10V=G_ISOt';
end

%V=5.5
for i=1:length(Pserb_55V)
    A_eff=A(1,1);
%Subsonico
if Pserb_55V(i)/Ps>b;
g_ISO=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*1e6*A_eff*sqrt(1-((Pserb_55V(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_55V=G_ISOt';
end


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