% ntubes.m
%
% This matlab script is used to calculate and plot the theoretical input
% impedance of a sequence of cylindrical sections.  The input impedance is
% calculated using transmission matrices.
%
% by Gary P. Scavone, CCRMA
% November 2001

% Physical Constants and evaluation frequencies

f = [1:4:5000];
omega = 2*pi*f;
c = 347.23;          % speed of sound (m/sec)
rho = 1.1769;        % density of air (kg/m^3)
eta = 1.846*10^(-3); % shear viscosity coefficient (kg/m sec)
gamma = 1.4017;      % ratio of specific heats
nu = 0.841;          % sqrt of Prandtl number
kappaOverCp = eta / nu^2;
k = omega/c;

% Viscous boundary layer thickness

dv = sqrt(2*eta ./ (rho*omega));

% --------------------------------------------------- %
% ADD/DELETE CYLINDRICAL ELEMENT VALUES HERE
% The section lengths should add up to about 17 cm
% for the vocal tract.

ra = [0.005, 0.003 0.03].';  % radii
L = [0.05 0.02 0.10].';     % lengths

% --------------------------------------------------- %

% Bore characteristic impedance

Ro = rho*c./(pi*ra.^2);

% Bore Losses (approximations for 300 K)

rv = ra*(rho*omega/eta).^(0.5);
Z_c = (Ro*ones(size(omega))).*(1 + 0.369*rv.^(-1))- ...
    j*(Ro*ones(size(omega))).*(0.369*rv.^(-1) ...
		+ 1.149*rv.^(-2) + 0.303*rv.^(-3));
alpha = (ones(size(ra))*k).*(1.045*rv.^(-1) + 1.08*rv.^(-2) + 0.75*rv.^(-3));
vpinv = (1 + 1.045*rv.^(-1))./c;
Gamma = alpha + j*((ones(size(ra))*omega).*vpinv);

% Lossy Cylindrical Transmission Matrix Elements

sinhL = sinh((L*ones(size(omega))).*Gamma);
coshL = cosh((L*ones(size(omega))).*Gamma);
A = coshL ;
B = Z_c.*sinhL ;
C = (Z_c.^(-1)).*sinhL ;
D = coshL ;

% Levine & Schwinger Radiation Impedance Approximation

z = k*ra(length(ra));
[r,l] = reflpoly(z);
R = -0.98*r.*exp(-2*j*z.*l);
Z_r = Ro(length(ra))*(1+R)./(1-R);

% Matrix Multiplications

MAS = 1;
MBS = 0;
MCS = 0;
MDS = 1;

for i=1:length(ra),
   MA = MAS.*A(i,:)+MBS.*C(i,:);
   MB = MAS.*B(i,:)+MBS.*D(i,:);
   MC = MCS.*A(i,:)+MDS.*C(i,:);
   MD = MCS.*B(i,:)+MDS.*D(i,:);
   MAS = MA;
   MBS = MB;
   MCS = MC;
   MDS = MD;
end

Zin = (MB + MA.*Z_r)./(MD + MC.*Z_r);

% Clear figure window, bring it to front, and plot

figure(1)
clf

% Plot Impedance

plot(f,abs(Zin)./Ro(1));
grid
xlabel('Frequency (Hz)')
ylabel('| Input Impedance | / Z0')
