//Calculation of dynamic response of a pressure tap point
//www.basicairdata.eu 2013 (c) JLJ
//Refer to notations, Pitot page on wwww.basicairdata.eu
clear ;
R=1e-3; //[m] Internal radius of pressure tap tubing
L=0.6;//[m] Line lenght
Vu=50e-9;//[m^3] Sensor internal volume
Vt=3.14*R^2*L;//[m^3] Line internal diameter
sigma=0;//
rhos=1.225; //[kg/m^3] Density of air at standard conditions
mu=1.78938e-5//[kg/m/s] or [Pas]Aboslute fluid viscosity
//mu=0;
a0=340;//[m/s] speed of sound
gammaratio=1.4;//Cp/Cv Specific heat ratio  
Pr=0.712;//Prandt Number
nu=[1:0.5:1200]';//Frequency value range radiants
[c r]=size(nu)
for i=1:c
alfa(i)=R*((rhos*nu(i)/mu)^0.5)*(complex(-2^0.5/2,2^0.5/2))
//nj(i)=1/((1+(gammaratio-1)/gammaratio*(besselj(2,alfa(i)*Pr^0.5)/besselj(2,alfa(i)*Pr^0.5))))
nj(i)=1.4;
phaselag(i)=nu(i)/a0*(besselj(0,alfa(i))/besselj(2,alfa(i)))^0.5*((gammaratio/nj(i))^0.5)
tfunction(i)=1/(cosh(phaselag(i)*L)+Vu/Vt*(sigma + 1/gammaratio)*nj(i)*phaselag(i)*L*sinh(phaselag(i)*L))
//tfunction(i)=(cos(nu(i)*L/a0)-gammaratio*Vu/Vt*(sigma+1/gammaratio)*nu(i)*L/a0*sin(nu(i)*L/a0))^-1
[phi(i),mag(i)]=phasemag(tfunction(i))
mag(i)=10^(mag(i)/20)
end
xsetech([0,0,1,0.5]);
plot2d(nu/6.28,mag')
xtitle('Frequency response, single pressure tap line www.basicairdata.eu JLJ')
xlabel('Frequency [Hz]')
ylabel('Magnitude P1/P0')
xgrid(1)
xsetech([0.,0.5,1,0.5]);
plot2d(nu/6.28,phi')
//xlabel('Frequency [Hz]')
ylabel('Phase lag [Degree]')
xgrid(1)