function dataReduction
% dataReduction  Convert viscometer data from MYO, Lab # 1-90
%
%  This m-file shows how to (1) store data in arrays, (2) perform array
%  calculations with .* and ./ operators, (3) plot data and label axes,
%  (4) sort data stored in different arrays, and (5) print data in a 
%  nicely formatted table.

% --- Store data
V = [ 9.2  9.7  9.2  9.1  9.2  9.4  9.1];   %  volume (mL)
t = [19.8 15.8 16.8 21.3 13.1 10.1  8.9];   %  time (s)
T = [15.6 26.3 21.3 12.3 34.3 50.4 58.1];   %  temperature (C)

% --- Begin data reduction
% Compute volumetric flow rate:  V/1.0e6 is volume in m^3
% Divide each volume measurement by time to fill volume:  Q = volume/time
% The ./ operator does element-by-element division of the V and t arrays.
Q = (V/1.0e6) ./ t;   %  volumetric flow rate (m^3/s)

% Store reference data for viscosity of water: MYO, Table B.1, p. 831
% Tnu is temperature (C), nuw is kinematic viscosity (m^2/s)
Tnu = [0      5      10     20     30      40      50      60];
nuw = [1.787  1.519  1.307  1.004  0.8009  0.6580  0.5534  0.4745] * 1e-6;

% --- Use one operating condition to determine viscometer calibration K
%     and compute nu at remaining points
% Linear interpolation to find viscosity of water at 15.6 (C)
nuref = nuw(3) + (T(1) - Tnu(3)) * (nuw(4) - nuw(3))/(Tnu(4) - Tnu(3)) 
% Note: nuref = interp1(Tnu,nuw,T(1))  works for any T(1) in range of Tnu data
K = nuref*Q(1)   % Use reference viscosity to determine K
nu = K ./ Q;     % Convert data.  The ./ does element-by-element division

% --- Plot results
plot(T,nu,'o',T(1),nuref,'+',Tnu,nuw,'r--');
xlabel('T ({}^\circ C)');
ylabel('\nu  (m^2/s)');
axis([0 60 1e-7 20e-7]);
grid on
legend('Measured','Calibration','Reference')

% --- Print results
%  Prepare by sorting data in order of increasing temperature
[T,is] = sort(T);   %  "is" is the sort order
nu = nu(is);        %  nu data is now in same order as T data
%  Interpolate to find reference viscosity at *all* measured temperatures
nur = interp1(Tnu,nuw,T);
fprintf('  T (C)    nu (m^2/s)  nu_ref (m^2/s) Diff (m^2/s)  %% Diff\n');
for i=1:length(T)
  Dabs = nu(i) - nur(i);
  Drel = 100*Dabs/nur(i);
  fprintf(' %6.1f  %12.3e  %12.3e  %12.3e  %6.2f\n',T(i),nu(i),nur(i),Dabs,Drel);
end
