%
% PURPOSE: Enter, process, and plot the calibration data.
% AUTHOR:  Rob McEwen
% DATE:    98/10/8
%
% Add lables below so that they match the order of the columns.
%
  lbls = { 'weight'
           'increase'
           'decrease'
         };
  sz = size(lbls);
  i = 1;
  while( i <= sz(1) ),
    eval([ char(lbls(i))  ' = '  num2str(i) ';']);
    i = i + 1;
  end

% Strain gauge up.
%          weight  increase decrease  log line #
  cal1   = [ 0      -.500  -.500
            8.5     -.580  -.610      %with plate and level
            31.0    -.819  -.855      %1
            55.5    -1.070 -1.101     %2
            79.6    -1.292 -1.334     %3
            104.8   -1.513 -1.555     %4
            127.7   -1.711 -1.752     %5
            151.2   -1.905 -1.982     %6
            171.7   -2.080 -2.154     %7
            195.7   -2.278 -2.357     %8
            218.45  -2.468 -2.551     %9
            240.55  -2.660 -2.745     %10
            261.35  -2.840 -2.940     %11
            285.35  -3.055 -3.173];
%
% NOTE: The 1,3 and 14,3 numbers were not measured, but were added just
% to fill out the matrix.
%
% Compute the slope and y intercept from least-squares.  Stack the 
% "increasing" and "decreasing" measurements to form a 2nx2 A matrix
% properly arranged for least squares.
%
%
  sz1    = size(cal1);
  A     = [cal1(:,weight) ones(sz1(1),1)];       
  A     = [A;A];                              
  y     = [cal1(:,increase);cal1(:,decrease)];
%
% Remove fake data point cal1(1,2):
% A     = A([1:14 16:27],:);
% y     = y([1:14 16:27]);
% sz1   = sz1 - 2;
%
  param = (A'*A)\A'*y;
  slope1  = param(1);     
  yInter1 = param(2);     
%
% Compute offset due to the weight of the structure and the duct, w0.
% The average reading with the stand horizontal is -.18.
%
  bias = -.201;    
  w1   = (yInter1 - bias)/abs(slope1);
  !rm cal.txt
  diary('cal.txt')
  disp(' ')
  disp([' The sensor bias on level ground is ', num2str(bias), ' volts.'])
  disp( ' ')
  disp(' Strain guage up: ')
  disp([' The slope is ',num2str(slope1),' Volts per pound.'])
  disp([' The y intercept is ',num2str(yInter1),' Volts.'])
  disp([' The weight with the duct and the level, but with the bias '])
  disp([' removed, is ',num2str(w1),' Pounds.'])


%
% Now repeat this for the other direction, that is, strain guage down.
%
%          weight  increase decrease  log line #
  cal2   = [ 0       .123   .07       %w/o plate, with level
            8.5      .181   .312      %with plate and level
            31.0     .401   .524      %1
            55.5     .652   .809      %2
            79.6     .895  1.029      %3
            104.8   1.149  1.250      %4
            127.7   1.378  1.450      %5
            151.2   1.615  1.666      %6
            171.7   1.825  1.824      %7
            195.7   2.063  2.062      %8
            218.45  2.290  2.289      %9
            240.55  2.510  2.512      %10
            261.35  2.714  2.718      %11
            285.35  2.951  2.958];
%
% Compute the slope and y intercept from least-squares.  Stack the 
% "increasing" and "decreasing" measurements to form a 2nx2 A matrix
% properly arranged for least squares.
%
  sz2    = size(cal2);
  A     = [cal2(:,weight) ones(sz2(1),1)];       
  A     = [A;A];                              
  y     = [cal2(:,increase);cal2(:,decrease)];
%
  param = (A'*A)\A'*y;
  slope2  = param(1);     
  yInter2 = param(2);     
%
% Compute offset due to the weight of the structure, the level, and duct, w0.  
%
  w2 = (yInter2 - bias)/abs(slope2);
  disp( ' ')
  disp(' Second case: strain guage down: ')
  disp([' The slope is ',num2str(slope2),' Volts per pound.'])
  disp([' The y intercept is ',num2str(yInter2),' Volts.'])
  disp([' The weight with the duct and the level, but with the bias '])
  disp([' removed, is ',num2str(w2),' Pounds.'])
%
% Now put the two (strain guage up, down) together:
%
  cal = [ flipud([-cal1(:,1) cal1(:,2:3)])
                   cal2                  ];
  yls1 = -slope1*cal(:,weight) + yInter1;  %Eqn of the line overlaying the data
  yls2 =  slope2*cal(:,weight) + yInter2;

%
% Subtract out the bias and the structural weight
%
  calbr1 =  flipud([-cal1(:,1)  cal1(:,2:3)-ones(sz1(1),2)*yInter1]);
  calbr2 =  [ cal2(:,1)  cal2(:,2:3)-ones(sz2(1),2)*yInter2];
  calbr = [calbr1;calbr2];
%
% Compute least squares linear fit of the unbiased measurements:
%
  szbr    = size(calbr);
  A     = [calbr(:,weight) ones(szbr(1),1)];       
  A     = [A;A];                              
  y     = [calbr(:,increase);calbr(:,decrease)];
  param = (A'*A)\A'*y;
  slopebr  = param(1);     
  yInterbr = param(2);     
  yls = slopebr*calbr(:,weight) + yInterbr;  
%
  disp(' ')
  disp(' De-Biased Combination of Both Cases ')
  disp([' The slope is ',num2str(slopebr),' Volts per pound.'])
  disp([' The y intercept is ',num2str(yInterbr),' Volts.'])
%
  error = [calbr(:,2)-yls;
           calbr(:,3)-yls];
  disp(' ')
  disp(' The average error between the least-squares fit and the')
  disp([' measured voltages is ', num2str( mean(error) ), ' volts.'])
  disp([' The standard deviation is ', num2str( std(error) ), ' volts'])
  disp([' or ', num2str(std(error)/slopebr), ' pounds.'])
  disp(' ')
  disp([' The sensor bias on level ground is ', num2str(bias/slopebr), ' lbs.'])

  diary off
%

  figure(1)
  clf

  plot(cal(:,weight),[cal(:,increase:decrease) yls1 yls2 ])
  ylabel('Volts')
  xlabel('Pounds')
  title('Raw Strain Data')
  grid on

  figure(2)
  clf

  plot(calbr(:,weight),[calbr(:,increase:decrease) yls ])
  ylabel('Volts')
  xlabel('Pounds')
  title('De-Biased Strain Data')
  grid on


%  figure(1)
%  clf
%
%  plot(cal1(:,weight),cal1(:,increase:decrease))
%  ylabel('Volts')
%  xlabel('Pounds')
%  title('Strain Gauge Up')
%
%  figure(2)
%  clf
%
%  plot(cal2(2:sz2,weight),cal2(2:sz2,increase:decrease))
%  ylabel('Volts')
%  xlabel('Pounds')
%  title('Strain Gauge Down')

