% beamdiff.m
%
% Quick check of shears in each beam.

nsub = 1;
profrange = [5500 6000];
min_snr = 3;
max_std = 3;
filename = 'hale4001';
ylim = [0 10];

[ts, vel, amp, dday] = adp(filename, 'ts', 'vel', 'amp', nsub, profrange);
%[ts, vel, amp, dday] = adp(filename, 'ts', 'vel', 'amp');

   figure
   orient tall


for ibeam =1:4;

   v = squeeze(vel(:,ibeam,:));
   a = squeeze(amp(:,ibeam,:));
   [ncells, nprofs] = size(v);

   a_noise = median(a(end,:))
   snr = a - a_noise;
   i_nosignal = find(snr < min_snr);
   v(i_nosignal) = NaN;

   dv = diff(v);
   cells = (1:(ncells-1)).';  % for plotting dv

   [dvm, dvstd, dvn] = mstdgap(dv.');
   dvm = dvm.'; dvstd = dvstd.'; dvn = dvn.';


if 0
   subplot(2,2,1)
   plot(dvm, cells, dvstd, cells)
   set(gca,'YDir', 'reverse', 'XLim', [-0.2 0.2]);
   set(gca,'YLim', ylim);
   grid
   ylabel('cell')
   xlabel('beam vel diff (m/s)')
   % legend('mean','std')
   title(['beam ' num2str(ibeam) ' snr > ' num2str(min_snr)])

   subplot(2,2,2)
   plot(dvn, cells)
   set(gca,'YDir', 'reverse');
   set(gca,'YLim', ylim);
   grid
   xlabel('N samples')
   ylabel('cell')
   title(['snr > ' num2str(min_snr)])
end

   i_outlier = find(abs(dv) > max_std * dvstd(2));

   dv(i_outlier) = NaN;

   [dvm, dvstd, dvn] = mstdgap(dv.');
   dvm = dvm.'; dvstd = dvstd.'; dvn = dvn.';

% Repeat the 3 std check:

   i_outlier = find(abs(dv) > max_std * dvstd(2));

   dv(i_outlier) = NaN;

   [dvm, dvstd, dvn] = mstdgap(dv.');
   dvm = dvm.'; dvstd = dvstd.'; dvn = dvn.';


   subplot(4,2,1 + (ibeam-1) * 2)
   plot(dvm, cells, dvstd, cells)
   set(gca,'YDir', 'reverse', 'XLim', [-0.2 0.2]);
   set(gca,'YLim', ylim);
   grid
   ylabel('cell')
   if ibeam == 4, xlabel('beam vel diff (m/s)'); end
   %legend('mean','std')
   title(['beam ' num2str(ibeam) '  max std: ' num2str(max_std)])

   subplot(4,2,2+(ibeam-1)*2)
   plot(dvn, cells)
   set(gca,'YDir', 'reverse');
   set(gca,'YLim', ylim, 'XLim', [0 nprofs]);
   grid
   if ibeam == 4, xlabel('N samples'); end
   ylabel('cell')
   title(['max std: ' num2str(max_std) ' applied twice'])

   date_it([filename ', profs ' num2str(profrange) ', nsub ' num2str(nsub)])

end


