%
% PURPOSE: Enter, process, and plot the calibration data.
% AUTHOR:  Rob McEwen
% DATE:    98/9/29
%
% 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

%          weight  increase decrease  log line #
  cal1   = [ 0       -.311  -.311
            8.0     -.445  -.444
            32.5    -.624  -.668      %1
            57.0    -.871  -.891      %2
            80.0    -1.093 -1.107     %3
            104.5   -1.342 -1.339     %4
            128.5   -1.593 -1.565     %5
            151.5   -1.797 -1.740     %6
            177.5   -2.044 -1.965     %7
            198.5   -2.225 -2.144     %8
            221.75  -2.410 -2.330     %9
            246.75  -2.623 -2.842     %10
            271.25  -2.874 -3.049     %11
            294.0   -3.060 -3.060];
%
% 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 points cal1(1,3) and cal1(14,3):
  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, w0.  The average (of two)
% readings with the stand horizontal is -.1645.
%
  bias = -.1645;    
  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(' First case: strain guage up: ')
  disp([' The slope is ',num2str(slope1),' Volts per pound.'])
  disp([' The y intercept is ',num2str(yInter1),' Volts.'])
  disp([' The tare weight, bias removed, is ',num2str(w1),' Pounds.'])


%
% Now repeat this for the other direction, that is, strain guage down.
%

  cal2 = [0       -0.69      -0.69           %w/o plate
          8.0       .023      -.014          %w/ plate, empty
          32.5      .216       .175          %1
          56.5      .419       .400          %2
          81.       .615       .632          %3
          104       .875       .907          %4
          128.5    1.102       1.170         %5
          151.5    1.355       1.411         %6
          172.5    1.595       1.644         %7
          197.5    1.870       1.873         %8
          223.5    2.110       2.127         %9
          246.75   2.360       2.380         %10
          271.25   2.627       2.644         %11
          294.0    2.867       2.867];       %12  
%
% 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.
%
%
  sz2    = size(cal2);
  A     = [cal2(:,weight) ones(sz2(1),1)];       
  A     = [A;A];                              
  y     = [cal2(:,increase);cal2(:,decrease)];
%
% Remove data points cal2(1,3) cal2(1,2), and cal2(14,3).  Assume that the
% first row of cal2 is an outlyer.
%
  A     = A([2:14 16:27],:);
  y     = y([2:14 16:27]);
% sz2   = sz2 -3;
%
  param = (A'*A)\A'*y;
  slope2  = param(1);     
  yInter2 = param(2);     
%
% Compute offset due to the weight of the structure, w0.  The average (of two)
% reading with the stand horizontal is -.1645.
%
  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 tare weight, bias removed,  is ',num2str(w2),' Pounds.'])
%
% Now put the two (strain guage up, down) together:
%
  cal = [ flipud([-cal1(:,1) cal1(:,2:3)])
           cal2(2:sz2,:)                  ];
  yls1 = -slope1*cal(:,weight) + yInter1;
  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(2:sz2(1),1)  cal2(2:sz2(1),2:3)-ones(sz2(1)-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')

