function [r,l] = reflpoly(z)
% REFLPOLY	Reflection coefficient magnitude and length correction of an
% unflanged cylindrical bore.
%
% [R.L] = REFLPOLY(Z) provides the reflection coefficient magnitude and
% length correction of an unflanged cylindrical bore, as given by Levine &
% Schwinger, for the values Z (= ka).
%
% This Matlab function first uses polynomial approximations to the
% reflection coefficient magnitude and phase (l/a) as determined with the
% script levinepoly.m.  The end correction polynomial is valid only for the
% range z <= 3.8317.  For values of z greater than 3.8317, the end
% correction is fitted with a smooth tail out to z = 20. This enables
% digital filter designs for frequencies out to 22.05 kHz and bore radii to
% 5 cm.
%
% By Gary P. Scavone, CCRMA
% Last Modified: February 1999

rPoly = [-4.498152195980034e-12  5.833298262887810e-10 ...
  -3.355168279123298e-08  1.127761635012252e-06  -2.456794111455048e-05  3.635273412899122e-04  -3.720768280670044e-03     2.627278681664893e-02    -1.246951247288215e-01  3.711841002387674e-01  -5.617225571053632e-01  -1.299579410009597e-02  1.000000000000000e+00];

loaPoly = [ -1.958328696574216e-03  2.114846071848918e-02 ...
  -9.051270191766737e-02  1.949901775176107e-01  -2.219012922642712e-01  1.498735073790306e-02   6.082823259482204e-01];

r = polyval(rPoly,z);

if max(z) <= 3.8317,
  l = polyval(loaPoly,z);
else
  i=1;
  while z(i) <= 3.8317,
    i = i + 1;
  end
  z1 = z(1:i-1);
  z2 = z(i:length(z));
  l = polyval(loaPoly,z1);
  scale = polyval(loaPoly,z(i-1))*(z(i-1)-3.0)^2;
  l = [l, scale./(z2 - 3).^2];
end
