function  fanCurveStats(fname,nfit,patm,dpatm,verbose)
% fanCurveStats  Obtain fan curve from flow bench measurements with statistical sampling
%
% Synopsis:  fanCurveStats
%            fanCurveStats(fname)
%            fanCurveStats(fname,nfit)
%            fanCurveStats(fname,nfit,patm)
%            fanCurveStats(fname,nfit,patm,dpatm,verbose)
%
% Input:  fname = (optional) name of MAT file containing averaged flow bench data.
%                 Default: fname = 'sampleAve.mat'.
%         nfit = (optional) degree of polynomial in lease squares fit to fan curve
%                Default:  nfit = 4
%         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, nfit = 4;       end
if nargin<3, patm = 760;     end
if nargin<4, dpatm = 0.3;    end
if nargin<5, 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,d.emfstd,...
                   d.rt,d.rtstd,d.vpp,d.vppstd,d.vpn,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 = polyfit(Q,dp,nfit);
Qfit = linspace(min(Q),max(Q));
dpfit = polyval(c,Qfit);                   %  evaluate fit for plotting
fprintf('\nCoefficients of fan curve polynomial: dp = sum of c(i)*Q^n\n')
fprintf('       c(i)          Q^n\n')
for i=1:length(c)
  fprintf('  %15.7e    Q^(%d)\n',c(i),nfit+1-i);
end

% --- Plot dp versus Q with error bars on both dp and Q
figure;
plot(Qfit,dpfit,'r--');                            %  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(sprintf('Degree %d fit',nfit),'Data',1)