%
% PURPOSE: Enter, process, and plot the calibration data.
% AUTHOR:  Rob McEwen
% DATE:    98/12/17
%
% 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
% Positive direction (Tension on the same side as the motor):
% Beware that Gary's box flips the sign of the conditioning amp.
%         weight increase    decrease  
  cal1   = [ 0, hex2dec('7fe'), hex2dec('7fe')
            20, hex2dec('839'), hex2dec('83b')
            40, hex2dec('889'), hex2dec('88a')
            60, hex2dec('8d8'), hex2dec('8e0')
            80, hex2dec('928'), hex2dec('92f')
           100, hex2dec('97b'), hex2dec('988')
           120, hex2dec('9c8'), hex2dec('9d9')
           140, hex2dec('a1c'), hex2dec('a2b')
           160, hex2dec('a6c'), hex2dec('A6d')
           180, hex2dec('aae'), hex2dec('abc')
           200, hex2dec('b03'), hex2dec('b01')
           220, hex2dec('b51'), hex2dec('b50')
           240, hex2dec('b9d'), hex2dec('b9a')
           260, hex2dec('be4'), hex2dec('be4')
           280, hex2dec('c1e'), hex2dec('c36')
           300, hex2dec('c81'), hex2dec('c81')];
%
% NOTE: The 1,2 and 16,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)];
  param = (A'*A)\A'*y;
  slope1  = param(1);     
  yInter1 = param(2);     
%
  !rm cal.txt
  diary('cal.txt')
  disp(' ')
  disp(' Positive Direction: ')
  disp([' The slope is ',num2str(slope1),' Counts per pound.'])
  disp([' The y intercept is ',num2str(yInter1),' Counts.'])
%
% Minus direction: (tension on opposite side of motor)
%         weight increase    decrease  
  cal2   = [ 0 hex2dec('800') hex2dec('7fe')
            20 hex2dec('7b2') hex2dec('7ac')
            40 hex2dec('767') hex2dec('760')
            60 hex2dec('716') hex2dec('712')
            80 hex2dec('6ca') hex2dec('6c3')
           100 hex2dec('676') hex2dec('671')
           120 hex2dec('619') hex2dec('621')
           140 hex2dec('5ca') hex2dec('5cc')
           160 hex2dec('582') hex2dec('57e')
           180 hex2dec('533') hex2dec('530')
           200 hex2dec('4e9') hex2dec('4e8')
           220 hex2dec('49c') hex2dec('498')
           240 hex2dec('449') hex2dec('444')
           260 hex2dec('3fc') hex2dec('3f8')
           280 hex2dec('3ac') hex2dec('3a9')
           300 hex2dec('358') hex2dec('358')];
%
% 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 the no-load bias, with the motor:
%
  bias = (yInter1 + yInter2)/2;
  disp( ' ')
  disp(' Second case: minus direction: ')
  disp([' The slope is ',num2str(slope2),' Counts per pound.'])
  disp([' The y intercept is ',num2str(yInter2),' Counts.'])
  disp([' The no-load bias is ',num2str(bias),' Counts'])
%
% Now put the two (strain guage up, down) together:
%
  cal = [ flipud([-cal2(:,1) cal2(:,2:3)])
                   cal1                  ];
  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
%
  calbrNeg =  flipud([-cal2(:,1)  cal2(:,2:3)-ones(sz2(1),2)*yInter2]);
  calbrPos =  [ cal1(:,1)  cal1(:,2:3)-ones(sz1(1),2)*yInter1];
  calbr = [calbrNeg;calbrPos];
%
% 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),' Counts per pound.'])
  disp([' The y intercept is ',num2str(yInterbr),' Counts.'])
%
  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) ), ' Counts.'])
  disp([' The standard deviation is ', num2str( std(error) ), ' Counts'])
  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('Counts')
  xlabel('Pounds')
  title('Raw Strain Data')
  grid on

  figure(2)
  clf

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

