function mu = airViscosity(T,units)
% airViscosity  Dynamic viscosity of air
%
% Synopsis:   mu = airViscosity
%             mu = airViscosity(T)
%             mu = airViscosity(T,units)
%
% Input:   T = (optional) temperature (deg. C) at which properties are
%              evaluated.  Default: T = 20C. If units='BG' then T is in deg. F
%          units = (optional) string to indicate the system of units
%                  units = 'SI' for System International (kg-m-s) units,
%                  units = 'BG' for British Gravitational (lbm-ft-s) units
%
% Output: mu  = dynamic viscosity, kg/(m*s) = Pa*s or lbm/(ft*s)
%
% Ref:  Incropera and Dewitt, Introduction to Heat Transfer, 3rd ed., Table A.4

if nargin<1         %  return viscosity at 20 C without evaluating curve fit
  mu = 1.81316e-05;
  return
end
if nargin==1,  units='SI';          end
if nargin>1 & isempty(T),  T = 20;  end

if units=='BG'
  Tin = (T-32)*5/9;    %  convert deg F to deg C
else
  Tin = T;
end
if Tin<-23.15 | Tin > 726.85
  error(sprintf('T = %f (deg C) is out of range for viscosity curve fit',Tin));
end
% --- Store polynomial curve fit coefficients
pm = [ 1.183165303908385553e-14 -2.838455962871034222e-11 ...
       4.921892043466875411e-08  1.715849963309301011e-05 ];
% --- Evaluate the polynomials
mu  = polyval(pm,Tin);

if units=='BG',  mu = mu*0.67197;   end