function  lossCoeffStats(fname,patm,dpatm,verbose)
% lossCoeffStats  Obtain system curve from flow bench measurements with statistical sampling
%
% Synopsis:  lossCoeffStats
%            lossCoeffStats(fname)
%            lossCoeffStats(fname,patm)
%            lossCoeffStats(fname,patm,dpatm)
%            lossCoeffStats(fname,patm,dpatm,verbose)
%
% Input:  fname = (optional) name of MAT file containing averaged flow bench data.
%                 Default: fname = 'sampleAve.mat'.
%         patm = (optional) atmospheric pressure, mm Hg, Default: patm = 760
%         dpatm = (optional) uncertainty in atmospheric pressure, mm Hg, Default: dpatm = 0.3
%         verbose = (optional) flag to control printing of reduced results.
%                   Default:  verbose = true;

if nargin<1, fname = 'sampleAve.mat';  end
if nargin<2, patm = 760;     end
if nargin<3, dpatm = 0.3;    end
if nargin<4, verbose = true; end

% --- Assign values not measured by the DAQ system and convert units as needed
patm = (patm/760)*101325;     %  atmospheric pressure, Pa
dpatm = (dpatm/760)*101325;   %  uncertainty in atmospheric pressure
ddn = 0.003*2.54e-2;          %  uncertainty in nozzle diameter

% --- Read data from MAT file
d = loadFBStatsData(fname,verbose);

% --- Convert raw data and print results
[Q,dQ,dp,ddp] = flowBenchStats(patm,dpatm,d.dn,ddn,d.emf,2*d.emfstd,...
                   d.rt,2*d.rtstd,d.vpp,2*d.vppstd,d.vpn,2*d.vpnstd);

if verbose
  fprintf('\nReduced Data\n     dpPlenum (Pa)         Q (m^3/s)\n');
  for i=1:length(dp)
    fprintf(' %8.2f +/- %-6.3f   %11.6f +/- %-9.6f\n',dp(i),ddp(i),Q(i),dQ(i));  
  end
end

% --- Least squares fit to find loss coefficient
c = sum(dp.*Q.^2)/sum(Q.^4);
Qfit = linspace(0,max(Q));
dpfit = c*Qfit.^2;                   %  evaluate fit for plotting
fprintf('\nLoss coefficient for system curve: dp = c*Q^2\n')
fprintf('     c = %e\n',c)
% --- Alternative curve fit
b = [Q(:).^2  Q(:)]\dp(:);
dpfit2 = b(1)*Qfit.^2 + b(2)*Qfit;
fprintf('\nAlternative loss coefficient model: dp = b(1)*Q^2 + b(2)*Q\n')
fprintf('     b = %e  %e\n',b)

% --- Plot dp versus Q with error bars on both dp and Q
figure;
plot(Qfit,dpfit,'r--',Qfit,dpfit2,'b:');           %  first plot data
hold on;  errorbarXY(Q,dp,dQ,ddp,'o');  hold off;  %  then add error bars
xlabel('Q  m^3/s');    ylabel('\Delta p   (Pa)');
legend('fit','alt fit','Data','Location','Northwest')