%To be changed: test deltap 450Pa
%jlj@BasicAirData.eu
%Script to plot time behaviour 
%Newsletter http:\www.basicairdata.eu\newsletter.htlm
nvol=30; %%Sets the square calculation grid
span=10;
vlimit=28e-6*span;
alimit=1e-10*span;
tf=60;
V=(3e-3)^2*3.14*1;
al=1e-10;
R=287.058;
T=288.15;
k=V/al/R/T;
p=101325;
p1=p+3400;
p2=p+400;
M1=p1*V/R/T;
M2=M1*p2/p1;
deltap=(M1-M2)/V*R*T;
%%deltam=(M1-M2)=deltap*V/(R*T);
tend=V/R/T/al*((p1-p)^0.5-(p2-p)^0.5);
p2attf=p+((p1-p)^0.5-tf/k)^2;
diffperc=(((p1-p)^0.5-tf/k)^2)/3600*100;
%%Calculates % variation of p2 respect 3600 Pa
Vl=linspace(vlimit,vlimit*10,nvol);
al=linspace(alimit,alimit*10,nvol);
for i = 1:nvol
    for j=1:nvol
    ki=Vl(1,i)/al(1,j)/R/T;
    if ((p1-p)^0.5-tf/ki)>0
    pointdiff(i,j)=(((p1-p)^0.5-tf/ki)^2)/3600*100;
    end
    end
end
figure(1);
%axis([0 vlimit 0 alimit 0 100]);
surf(Vl,al,pointdiff);
xlabel('$V[m^3]$','Interpreter','latex')
ylabel('$a_l[\frac{Kg}{s\sqrt{Pa}}$]','Interpreter','latex')
title('Pressure Loss Percentage'); 

figure(2);
contourf(Vl,al,pointdiff,15);
%%axis([0 vlimit*span 0 alimit*span]);
xlabel('$V[m^3]$','Interpreter','latex')
ylabel('$a_l[\frac{Kg}{s\sqrt{Pa}}$]','Interpreter','latex')
title('Pressure Loss Percentage'); 

%%Pass not pass, plot al, v plane
figure(3)
psw=linspace(101325,110000,nvol);
i=1;
j=1;
hold on;
%%for i=1:nvol     
for j=1:nvol
    %Vout(1,j)=al(1,j)*(tf*R*T)/((psw(1,i)-p)^0.5);
    Vout(1,j)=al(1,j)*(tf*R*T)/((p1-p)^0.5);
end
%plot(Vout,al)
stem(al,Vout)
xlabel('$V[m^3]$','Interpreter','latex')
ylabel('$a_l[\frac{Kg}{s\sqrt{Pa}}$]','Interpreter','latex')
title('Existence Surface'); 

%%end
hold off;

%export_fig -r600 with export_fig matlab pippone