function Q = nozzleFlow(d,dplen,dpn,p,T)
% nozzleFlow  Volumetric flow rate of air through a long radius nozzle.
%
% Synopsis:  Q = nozzleFlow(d,dplen,dpn)
%            Q = nozzleFlow(d,dplen,dpn,p)
%            Q = nozzleFlow(d,dplen,dpn,p,T)
%
% Input: d = throat diameter,  (m)
%        dplen = upstream pipe (or plenum) diameter (m)
%        dpn = pressure drop across nozzle, (Pa)
%        p = (optional) pressure upstream of nozzle;  Default: p = 101325 (Pa)
%        T = (optional) temperature upstream of nozzle;  Default: T = 20 (C)
%
% Output: Q = volumetric flow rate, (m^3/s)
if nargin<4,  p = 101325;  end   %  standard atmosphere, (Pa)
if nargin<5,  T = 20;      end   %  degrees C
if dpn<eps || d<eps,  Q=0;  return;  end  %  allows for no-flow case on fan curves

% --- Evaluate fluid properties and other constants
mu = airViscosity(T);        %  kinematic viscosity
rho = p/(287*(T+273.15));    %  air density from ideal gas law
bbeta = d/dplen;
y = expansionFactor(p,dpn,bbeta,1.4);
area = 0.25*pi*d^2;
qcon = area*y*sqrt(2*dpn/(rho*(1-bbeta^4)));  rcon = rho*d/(area*mu);

% --- Initialize and loop until cd converges
tol = 5e-6;  it = 0;  maxit = 25;   cdold = 0;  cd = 0.9;
while abs(cdold-cd)>tol && it<maxit
  cdold = cd;
  Q = cd*qcon;
  Re = rcon*Q;
  cd = 0.9986 - 7.006/sqrt(Re) + 134.6/Re;
  it = it + 1;  
end
if it>=maxit
  warning(sprintf('cd computations not converged after %d iterations',it));
end
