function [a, ph, good_mask] = btcaluv(cal_, d_, label, step, trange_arg)
%  function [a, ph, good_mask] = btcaluv(cal, d, label, step, trange_arg)
%   input:
%     cal      = the matrix output by refabsbt; or output file name
%     d        = the matrix output by lst_btrk; or output file name
%     label    = a plot title element
%   optional input:
%     step     = 1 to use individual ensembles; 2 to average
%                adjacent ones; etc. default is 1. (see note below)
%     trange_arg = vector of start, end times desired  (decimal day)
%
%   output:
%     a        = the individual amplitude factor estimates
%     ph       = the individual angle estimates (degrees)
%     good_mask= editing mask for these
%
%  Conclusions are appended to the file "btcaluv.out".
%  Plot is written to "btcaluv.ps".

%     NOTE: recommended step size:
%	for PCODE with 5 minute ensembles, use step = 1
%	for GPS with SA, try step =3  
%     As the step size increases the smoothing causes a decrease
%     in variance, but also decreases the ability to detect outliers
%     and results in fewer independent samples.
%
%                   Eric Firing
%
%     95/11/24 Original version.
%		96/05/06 - EF - major modifications for convenience,
%              improved editing and plotting.
%     97/05/16 - EF - added optional time range input argument.
%     97/05/20 - JF - additional documentation added.

min_speed = 2;
max_sig = 1;
max_gap = 0.1; % max gap in minutes. (For column 8 in refabsbt output.)
tol_dt = 0.02; % fractional tolerance for sample time interval

min_depth = 25;
max_depth = 400;

if nargin < 5,
   trange_arg = [];
end
if nargin < 4, 
        step = 1, 
end
if nargin > 5 | nargin < 3,
        help btcaluv
        error('3, 4, or 5 inputs are required')
end

if isstr(cal_),
   cal = read_asc(cal_, 8);
else
   cal = cal_;
end

if isstr(d_),
   d = read_asc(d_, 4);
else
   d = d_;
end

%-----------------------------------------------------------
% Get what we need from lst_btrk output; partly for plotting,
% partly for selecting the desired depth range.
[nd, nn] = size(d);

i_d = 1:nd; 
if ~isempty(trange_arg)
   i_d = find(d(:,1) >= trange_arg(1) & d(:,1) <= trange_arg(2));
   nd = length(i_d);
end


day_bt = d(i_d,1);
u_bt   = d(i_d,2);
v_bt   = d(i_d,3);
depth_bt = d(i_d,4);
clear d i_d

uv_bt = u_bt + i * v_bt;
heading_bt = 90 - angle(uv_bt)*180/pi;
ii = find(heading_bt > 180);
if ~isempty(ii), heading_bt(ii) = heading_bt(ii) - 360; end
speed_bt = abs(uv_bt);
% End of lst_btrk output extraction.

%---------------------------------------------------------------
[ndcal, nn] = size(cal);
i_d = 1:ndcal; 
if ~isempty(trange_arg)
   i_d = find(cal(:,1) >= trange_arg(1) & cal(:,1) <= trange_arg(2));
   cal = cal(i_d,:);
end


% At the outset, eliminate any points with significant gaps in the BT:
i_gap = find(cal(:,8) > max_gap);
cal(i_gap,:) = [];

% Eliminate any points outside the desired depth range:
% Estimate the depths at the fix times from the lst_btrk depths.
% Normally, the fix times will nearly, but not perfectly, coincide
% with the BT times.  Hence the interpolation step is required, but
% will not greatly change the depth estimates.
i_no_depth = find(cal(:,1) < day_bt(1) | cal(:,1) > day_bt(nd));
cal(i_no_depth, :) = [];

dtab = [day_bt depth_bt];
cal_depth = table1(dtab, cal(:,1));
i_wrong_depth = find(cal_depth < min_depth | cal_depth > max_depth);
cal(i_wrong_depth,:) = [];

% With this reduced cal array, select the indices for the BT versus
% nav comparison.  Velocities will be averaged from i0 to i1.
[nfix, nn] = size(cal);
i0 = 1:(nfix - step);
i1 = i0 + step;


day = cal(:,1);   % Fix times in decimal days.
dayuv = 0.5 * (day(i0) + day(i1));
trange = day([1 nfix]);

% 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.


xynav = cumsum(cal(:,2) + i * cal(:,3));
xybt  = cumsum(cal(:,4) + i * cal(:,5));

dt = (day(i1) - day(i0)) * 86400;
uvnav = (xynav(i1) - xynav(i0)) ./ dt; 
uvbt  = (xybt(i1)  - xybt(i0))  ./ dt; 
med_dt = median(dt);

n_raw = length(dt);

speed_nav = abs(uvnav);

aa = uvnav./uvbt;

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

good_mask = speed_nav > min_speed & ...
               abs(dt - med_dt)/med_dt < tol_dt;

subplot(2,2,1),hist(a(good_mask),20)
title('Amplitude prior to filtering outliers')
subplot(2,2,2),hist(ph(good_mask),20)
title('Phase prior to filtering outliers')

for ii = 1:3
   a_mode = mode(a(good_mask),20);
   a_std = std(a(good_mask));
   a_res = a - a_mode;

   ph_mode = mode(ph(good_mask),20);
   ph_std = std(ph(good_mask));
   ph_res = ph - ph_mode;


   good_mask = good_mask & ...
                 abs(a_res) < max_sig*a_std & ...
                 abs(ph_res) < max_sig*ph_std;
end

subplot(2,2,3),hist(a(good_mask),20)
title('Amplitude after to filtering outliers')
subplot(2,2,4),hist(ph(good_mask),20)
title('Phase after to filtering outliers')
print -dps btcaluv_hist

a_mn_e = mean(a(good_mask));
a_med_e = median(a(good_mask));
a_std_e = std(a(good_mask));
a_res_e = a(good_mask) - a_mn_e;

ph_mn_e = mean(ph(good_mask));
ph_med_e = median(ph(good_mask));
ph_std_e = std(ph(good_mask));
ph_res_e = ph(good_mask) - ph_mn_e;

n_ed = length(find(good_mask));

%---------------------------------------------------------
% Write the results out to a file named "btcaluv.out".

fid_out = fopen('btcaluv.out', 'a');
fprintf(fid_out, [label '\n   Time range %6.2f to %6.2f'], ...
                     trange(1), trange(2));
fprintf(fid_out, ['\n   Calculation done at ' pctime ]);
fprintf(fid_out, '\n    step: %d', step);
fprintf(fid_out, '\n    min_depth: %3.0f   max_depth: %3.0f', ...
                   min_depth, max_depth);
fprintf(fid_out, '\n    min_speed: %4.1f m/s   max_sig: %4.1f std devs', ...
                   min_speed, max_sig);
fprintf(fid_out, '\n    max_gap: %5.2f minutes   tol_dt: %4.2f (fraction)', ...
                   max_gap, tol_dt);
fprintf(fid_out, '\nunedited: %d points', n_raw);
%fprintf(fid_out, '\n            median     mean      std');
%fprintf(fid_out, '\namplitude  %7.4f  %7.4f  %7.4f',      ...
%                     a_med, a_mn, a_std);
%fprintf(fid_out, '\nphase      %7.4f  %7.4f  %7.4f',      ...
%                     ph_med, ph_mn, ph_std);
%
fprintf(fid_out, '\nedited:   %d points, %3.1f min speed, %3.1f max dev',...
                                                                                         n_ed, min_speed, max_sig);
fprintf(fid_out, '\n            median     mean      std');
fprintf(fid_out, '\namplitude  %7.4f  %7.4f  %7.4f',      ...
                     a_med_e, a_mn_e, a_std_e);
fprintf(fid_out, '\nphase      %7.4f  %7.4f  %7.4f',      ...
                     ph_med_e, ph_mn_e, ph_std_e);
fprintf(fid_out, '\n\n');
fclose(fid_out);


%---------------------plotting section-------------------------%
t0 = min(day);
t1 = max(day);

subplot(6,1,1)
plot(dayuv(good_mask), a(good_mask), '*');
%axis([t0 t1 (a_mn-0.01) (a_mn+0.01)]);
axis([t0 t1 -Inf Inf]);
grid
ylabel('scale');
title(label)

subplot(6,1,2)
plot(dayuv(good_mask), ph(good_mask), '*');
%axis([t0 t1 (ph_mn-0.5) (ph_mn+0.5)]);
axis([t0 t1 -Inf Inf]);
grid
ylabel('Angle')

subplot(6,1,3)
plot(day_bt, speed_bt, '+');
axis([t0 t1 -8 8]);
grid
ylabel('Speed')

subplot(6,1,4)
plot(day_bt, heading_bt, '+');
axis([t0 t1 -200 200]);
set(gca, 'YTick', [-180 -90 0 90 180]);
grid
ylabel('Heading')

subplot(6,1,5)
plot(day_bt, depth_bt, '+');
%set(gca, 'YScale', 'log')
axis([t0 t1 0 400]);
grid
set(gca, 'YDir', 'reverse');
ylabel('Depth')
xlabel('Decimal Day')

date_it(sprintf('btcaluv.m step: %d, mean ph: %7.4f, amp: %7.4f', ...
                                 step, ph_mn_e, a_mn_e));

subplot(6,1,6)
plot(day_bt(1:length(day_bt)-1)+0.5*diff(day_bt), ...
     24*3600*diff(depth_bt)./diff(day_bt), '+');
set(gca,'xlim',[t0 t1]);
grid
ylabel('d(Depth)/dt')
xlabel('Decimal Day')
date_it(sprintf('btcaluv.m step: %d, mean ph: %7.4f, amp: %7.4f', ...
                                 step, ph_mn_e, a_mn_e));

orient('tall')
print btcaluv.ps
