function [a, ph, res, ind] = btcal(cal, label, idelete)
%  function [a, ph, res, ind] = btcal(cal, label, idelete)
%        input:
%     cal      = the matrix output by refabsbt
%     label    = a plot title element
%     idelete  = array of rows to delete
%        output:
%     a        = the least squares best fit amplitude factor
%     ph       = the least squares best fit angle (degrees)
%     res      = the position residuals of the fixes to the best fit
%     ind      = indices in cal of rows used in the fit
%     E.F.        Wed  09-27-1989
%
% the heart of this program is the line aa = bt \ nav, which is a
% (least squares) fit of the bottom track and the sat fix displacements

% i = sqrt(-1); no longer needed in Matlab 3.5

[nfix, nn] = size(cal);

day = cal(:,1);   % Fix times in decimal days.
ind = 1:nfix;     % This will maintain the original fix indices,
                  % regardless of deletions.

% cal(:,2) and cal(:,3) are the differences in meters between successive
%   fixes in the x and y directions, respectively.  Row 1 is zero, row
%   2 is the displacement from time 1 to time 2, etc.  cal(:,4) and cal(:,5)
%   are similar, but determined from ADCP bottom track.
% nav and bt will be the displacements (in the complex plane, where
%   the imaginary part is north) of the ship since the time of the first fix;
%   nav is this displacement as determined from the fixes, and bt is
%   as determined from the bottom track.
nav = cumsum(cal(:,2) + i * cal(:,3));
bt = cumsum(cal(:,4) + i * cal(:,5));

if (nargin > 2)   % Delete rows listed in input array "idelete".
   if (~isempty(idelete))
      nav(idelete) = [];
      bt(idelete) = [];
      day(idelete) = [];
      ind(idelete) = [];
   end
else              % There was no input "idelete", so create an empty one.
   idelete = [];
end

% Scaling and rotation must be done relative to the "middle" of
%  the cruise track, as defined by the mean of the fixes.
nav = nav - mean(nav);
bt = bt - mean(bt);
% Now find best fit coefficient (aa) for  nav = aa * bt + error.
aa = bt \ nav;
res = nav - aa * bt;

a = abs(aa);
ph = angle(aa) * 180/pi;

% Now everything has been calculated, and we are just displaying the results.
% First, build a string which lists the indices of the deleted fixes.
if (~isempty(idelete))
   nd = length(idelete);
   idstring = int2str(idelete(1));
   if nd>1
      for ii = 2:nd;
         idstring = [idstring  ','  int2str(idelete(ii))];
      end
   end
else
   idstring = 'none';
end

% Write the results out to a file named "btcal.out".
fprintf('btcal.out', ...
  ['\n', label, ' (deleting ', idstring, '); amp= %7.4f, ph= %7.3f\n'], a, ph);
nres = length(res);
fprintf('btcal.out', '           day       dx       dy\n');
for ii = 1:nres;
   fprintf('btcal.out', '  %12.5f  %7.0f  %7.0f\n',...
                     day(ii), real(res(ii)), imag(res(ii)));
end

% Plot the residuals and save the plot.
plot(res)
title([label, ' (deleting ', idstring, ') cal residuals'])
xlabel('x (m)')
ylabel('y (m)')
ax=axis;
text(-0.15, -0.1, pctime, 'units', 'normalized', 'rotation', 90.0);
text(ax(1), ax(4), ['amp= ' num2str(a) ' ph= ' num2str(ph)], ...
  'VerticalAlignment', 'top')
print btcal.ps
