function myFanCurves(nfit1,nfit2,verbose)
% myFanCurves  Simple m-file function to put two fan curves on a single plot
%
% Synopsis:  myFanCurves
%            myFanCurves(nfit1)
%            myFanCurves(nfit1,nfit2)
%            myFanCurves(nfit1,nfit2,verbose)
%
% Input:  nfit1 = degree of curve fit for data set 1.  Default: nfit1=4
%         nfit2 = degree of curve fit for data set 2.  Default: nfit2=4
%         verbose = flag to control printing of raw and reduced data.
%                   Default: verbose = true,  print the data
%
% Output:  Print out of curve fit coefficients, plot of fan curve, and
%          optional print out of raw and reduced data

if nargin<1 || isempty(nfit1), nfit1=4;      end
if nargin<2 || isempty(nfit2), nfit2=nfit1;  end
if nargin<3,  verbose = true;  end

% --- Specify the file names and degree of polynomial curve fits
fname1 = 'Panaflow20V.dat';  patm1 = (752/760)*101325;  Rz1 = 0.5*(12002+11980);
fname2 = 'Panaflow24V.dat';  patm2 = (749/760)*101325;  Rz2 = 0.5*(11875+11750);
dplen = 22*2.54e-2;

% --- Read data and let flowBench reduce it
[dn,vdpp,vdpn,emf] = loadFlowBenchData(fname1,verbose);
[Q1,dp1] = flowBench(dn,vdpp,vdpn,emf,Rz1,patm1,dplen,verbose);

[dn,vdpp,vdpn,emf] = loadFlowBenchData(fname2,verbose);
[Q2,dp2] = flowBench(dn,vdpp,vdpn,emf,Rz2,patm2,dplen,verbose);

% --- Plot fan curve data from measurements
plot(Q1,dp1,'o',Q2,dp2,'*');
xlabel('Q  (m^3/s)');    ylabel('\Delta p  (Pa)');   legend(fname1,fname2);

% --- Obtain curve fits and print coefficients to the command window
c1 = polyfit(Q1,dp1,nfit1);
fprintf('\nCurve fit coefficients for %s\n',fname1);
for i=1:length(c1)
  fprintf('\t%20.11e  Q^(%d)\n',c1(i),nfit1-i+1);
end
c2 = polyfit(Q2,dp2,nfit2);
fprintf('\nCurve fit coefficients for %s\n',fname2);
for i=1:length(c1)
  fprintf('\t%20.11e  Q^(%d)\n',c2(i),nfit2-i+1);
end

% --- Evaluate curve fit over the range of data (include Q=0 as lower limit)
Qfit1 = linspace(0,max(Q1));   dpfit1 = polyval(c1,Qfit1);
Qfit2 = linspace(0,max(Q2));   dpfit2 = polyval(c2,Qfit2);

% --- Add curve fits without erasing previous plot
hold('on');  plot(Qfit1,dpfit1,'r--',Qfit2,dpfit2,'b--');  hold('off');