//www.BasicAirData.eu (c)2014 JLJ
//Adiabatic process
//http://en.wikipedia.org/wiki/Adiabatic_process
clear all;
xdel(winsid())

gammaair=1.4; //Table
ri=3e-3;//Internal pipe radius
lt=0.4;//lenght of tubing
D=21e-3;//Internal diameter of cylinder
L=25e-3;//Cylinder full stroke
screwpitch=1e-3;//Selected screw pitch
Vcil=(D)^2*L/4*3.14;//Volume of cylinder
P1=101325;//Initial pressure Pa
V1=ri^2*3.14*lt+Vcil;
V2=V1-Vcil;
ratio=V2/V1;
P2=P1*(V2/V1)^(-gammaair)
deltap=P2-P1
Vcil=(D)^2*screwpitch/4*3.14;//Volume of cylinder after one screw turn
V1=ri^2*3.14*lt+Vcil;
V2=V1-Vcil;
Pturn=P1*(V2/V1)^(-gammaair)
printf("Cylinder internal diameter %.3f mm\n",1000*D)
printf("Cylinder full stroke lenght %.3f mm\n",1000*L)
printf("P1 Initial pressure %.0f Pa\n",P1)
printf("P2 terminal pressure %.0f Pa\n",P2)
printf("Maximum differential pressure %.0f Pa\n",deltap)
printf("Differential pressure after one knob turn %.0f Pa\n",Pturn-P1)
printf("P2/P1 %.1f\n", P2/P1)

//Plot section Sizing chart stroke-diameter, deltap constant
fig1=scf(0)

ncurves=15; //Number of curves to plot
nlenght=50;//Number of different lenghts to plot
il=linspace(5e-3,50e-3,nlenght) //A set of piston stroke lenght values
deltapset=linspace(1000,150000,ncurves) //A set of deltap values
P1=101325*ones(1,ncurves);//Initial pressure Pa

P2=deltapset+P1 //Pressure after compression
V1=(ri^2*3.14*lt)*ones(1,ncurves)
V2=V1.*((P1./P2)^(1/gammaair))

deltav=V1-V2;
//D=(deltav*4/3.14/L)^0.5
//ROW contains different lenghts values
//COLUMN contains different deltaps values
//row=1;
//column=1;

for column=1:ncurves
    row=1
    for len=il
        D(column,row)=(deltav(1,column)*4/3.14/len)^0.5;
        row=row+1
    end
end
for idiff=1:ncurves
    toutp(:,idiff)=1000*D(idiff,:)'
end
plot(1000*il,toutp)
xlabel("Piston stroke lenght [mm]","fontsize",4)
ylabel("Piston internal diameter [mm]","fontsize",4)
title("Sizing chart, constant differential pressure for full stroke","fontsize",4)
eti=string(floor(deltapset))+" Pa"
legend(eti);

fig2=scf(1)
//Plot section Sizing chart stroke-diameter, deltap per turn

ncurves=10; //Number of curves to plot
nlenght=50; //Number of therad pitchs to consider
il2=linspace(0.25e-3,2e-3,nlenght) 
//il=[0.25e-3,0.5e-3,0.75e-3,1e-3,1.25e-3,1.5e-3]//A set of thread lenght

deltapset2=linspace(500,5000,ncurves) //A set of deltap values
P12=101325*ones(1,ncurves);//Initial pressure Pa

P22=deltapset2+P12 //Pressure after compression
V12=(ri^2*3.14*lt)*ones(1,ncurves)
V22=V12.*((P12./P22)^(1/gammaair))

deltav2=V12-V22;
//D=(deltav*4/3.14/L)^0.5
//ROW contains different lenghts values
//COLUMN contains different deltaps values
//row=1;
//column=1;

for column=1:ncurves
    row=1
    for len=il2
        D2(column,row)=(deltav2(1,column)*4/3.14/len)^0.5;
        row=row+1
    end
end
for idiff=1:ncurves
    toutp2(:,idiff)=1000*D2(idiff,:)'
end
plot(1000*il2,toutp2)
xlabel("Screw pitcht [mm]","fontsize",4)
ylabel("Piston internal diameter [mm]","fontsize",4)
title("Sizing chart, constant differential pressure for one knob turn","fontsize",4)
eti=string(deltapset2)+" Pa"
legend(eti);
//Closed form solution
desiredpturn=4000 //Pa
desireddeltap=50000//Pa
desiredD=21e-3;//m desired piston internal diameter
P12=101325 //Initial pressure
//pressure, one knob turn
P22=desiredpturn+P12
V12=(ri^2*3.14*lt)
V22=V12*((P12/P22)^(1/gammaair))
deltav=V12-V22 //Piston volume
screwpitch=4*deltav/(3.14*desiredD^2)
//pressure, full stroke
P22=desireddeltap+P12
V12=(ri^2*3.14*lt)
V22=V12*((P12/P22)^(1/gammaair))
deltav=V12-V22 //Piston volume
lrequired=(4*deltav/(3.14*desiredD^2))