function [Q,hL,hs] = pipeFlowSolve(Q0,L,A,Dh,e,KL,Am,nu,dz,pc)
% pipeFlowSolve  Balance point for a single pipe with or without a pump
%
% Synopsis:  [Q,hL,hs] = pipeFlowSolve(Q0,L,A,Dh,e,KL,Am,nu,dz,pc)
%
% Input:  Q0 = initial guess at pipe flow rate (m^3/s)
%              If Q0 is a scalar, fzero will try to find a bracket.
%              If Q0 is a two element vector, Q0(1) must be lowest possible
%              flow rate, and Q0(2) must be highest possible flow rate
%         L = vector of pipe lengths
%         A = vector of cross-sectional areas of ducts.  A(1) is area of
%             pipe with length L(1) and hydraulic diameter Dh(1)
%         Dh = vector of pipe diameters
%         e = vector of pipe roughnesses
%         KL = vector of minor loss coefficients
%         Am = vector of areas associated with minor loss coefficients.
%              For a flow rate, Q, through minor loss element 1, the area,
%              Am(1) gives the appropriate velocity from V = Q/Am(1).
%              Thus, different minor loss elements will have different
%              velocities for the same flow rate.  For a sudden expansion
%              the characteristic velocity is the upstream velocity, so
%              Am for that element is the area of the upstream duct.
%         nu = kinematic viscosity of the fluid 
%         dz = change in elevation, z_out - z_in;  dz>0 if outlet is above
%              the inlet;  dz = 0 for closed loop
%         pc = vector of polynomial coefficients defining the pump curve.
%              Note that the polynomial is defined in decreasing power of Q:
%               hs(Q) = pc(1)Q^n + pc(2)*Q^(n-1) + ... + pc(n)*Q + pc(n+1)
%              If there is no pump in the system, use pc = [] (null matrix)
%
% Output: Q  = flow rate through the pipeline
%         hL = head loss in the system
%         hs = shaft work done on the fluid

  % -- Embedded function for fzero
  %    pipeHeadBal returns the error in the head balance when a flow rate
  %    is guessed.  When this function returns zero (or a number less than
  %    the convergence tolerance, the correct flow rate has been guessed
  %    Other inputs to the pipeLoss() function are shared with the main
  %    function (pipeFlowSolve) as if they are global values.  hs needs
  %    to be computed each time because the head created by the pump
  %    depends on the flow rate.
  function h = pipeHeadBal(Q)
    if isempty(pc)
      hs = 0;               %  no pump curve data -> no pump
    else
      hs = polyval(pc,Q);   %  evaluate polynomial model of pump curve at Q
    end
    h = hs - dz - pipeLoss(Q,L,A,Dh,e,KL,Am,nu);
  end

% --- Find balance point
Q = fzero(@pipeHeadBal,Q0,optimset('Display','none'));

if isempty(pc)
  hs = 0;
else
  hs = polyval(pc,Q);     %  compute pump head at the balance point
end
fprintf('\nSystem balance point at\n');
fprintf('\t Q = %8.4f m^3/s    pump head = %6.1f (m)\n',Q,hs);

% --- Print results
[hv,hm,f] = pipeLoss(Q,L,A,Dh,e,KL,Am,nu);  %  re-evaluate losses for printout
hL = sum(hv) + sum(hm);                     %  total loss

fprintf('\nViscous losses in straight pipe sections\n');
fprintf(' section        Re           f         hL\n');
for k=1:length(L)
  Re = (Q/A(k))*Dh(k)/nu;
  fprintf('  %3d    %12.2e    %8.4f   %6.1f\n',k,Re,f(k),hv(k));
end
fprintf('\nMinor losses in straight pipe sections\n');
fprintf(' fitting     KL        hm\n');
for k=1:length(KL)
  fprintf('  %3d   %8.2f    %6.1f\n',k,KL(k),hm(k));
end
fprintf('\nTotal loss in straight sections  =     %6.1f  m\n',sum(hv));
fprintf('Total minor losses               =     %6.1f  m\n',sum(hm));
fprintf('Total viscous loss               =     %6.1f  m\n',hL);
fprintf('Total losses plus elevation gain =');
if abs(hL+dz)<0.005
  fprintf('%11.2e  m  (should equal pump head)\n',hL+dz);
else
  fprintf('     %6.1f  m  (should equal pump head)\n',hL+dz);
end

if abs(hL+dz-hs)>0.00005*hL     % --- error checking
  fprintf('\n\n\tPotential problem in pipeSysFlow\n');
  fprintf('\tpump head    = %g not equal to\n',hs);
  fprintf('\tsystem head = %g at alleged balance point\n',hL+dz);
  h = pipeHeadBal(Q,L,A,Dh,e,KL,Am,nu,dz,pc);
  fprintf('\n\tHead imbalance in energy equation is %11.2e (m)',h);
  fprintf(' at balance point\n');
end

end
