//Pitot Simulation
//www.basicairdata.eu 2013 (c) JLJ
clear ;
//Setup variable
graphics=0; //{1,0}  graphics plots on/off
//Pitot calibration factor c
cpitot=0.9965790;
//Model parameters
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.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))
[phi(i),mag(i)]=phasemag(tfunction(i))
//mag(i)=10^(mag(i)/20)
end


fdt=frep2tf(nu/6.28,tfunction,2);
fdt=factors(fdt)
gain=syslin('c',1/0.9963631,1); //set the dc gain case 1
q=fdt*gain; //Pressure line transfer function


if graphics==1 then
for i=1:c
[phi2(i),mag2(i)]=phasemag(repfreq(fdt,nu(i)/6.28))
end
end
//Launch xcos
xcos('pitotdyn3.zcos')

//Visualitation tasks
if graphics==1 then
xsetech([0,0,1,0.5]);
plot2d(nu/6.28,mag',1)
plot2d(nu/6.28,mag2',2)
xtitle('Frequency response, single pressure tap line  Original tf (black) Vs 2nd order identified tf (blue) www.basicairdata.eu JLJ')
xlabel('Frequency [Hz]')
ylabel('Magnitude OUT/IN [dB]')
xgrid(1)
xsetech([0.,0.5,1,0.5]);
plot2d(nu/6.28,phi',1);
plot2d(nu/6.28,phi2',2);
xlabel('Frequency [Hz]')
ylabel('Phase lag [Degree]')
xgrid(1)
end
