% compare sfc currents from uppermost ADCP bin with currents
% recomputed from the doppler speed log.

%% initialize
clear
close all

% matlab script path
path_prog = fullfile('c:\Users\Sergio\Documents\MATLAB\2016_CAPRICORN\check_ship_data');
cd(path_prog);                  % cd to matlab directory
restoredefaultpath              % reset matlab path
addpath(genpath(path_prog));	% add scripts to path
savepath; rehash;

fig_path = 'D:\Data\CAPRICORN_2016\Investigator\flux\SCS\nav_analysis';
mat_path = 'D:\Data\CAPRICORN_2016\Investigator\flux\SCS\nav\mat\CAPRICORN_2016_nav10.mat';

stjd = 89+10/24; endjd = 90;  % day range for plots
jd1 = sprintf('%03i',stjd);
jd2 = sprintf('%03i',endjd);

%% load recomputed doppler spd log data
load(mat_path); % load nav10 array of 10 min spd/course data

% nav10 13 columns - Dspd/Ddir is estimate for surface current
% jd Ngnd  Egnd  Nsfc  Esfc  Ndft  Edft  Gspd  Gdir  Sspd  Sdir  Dspd  Ddir
%  1   2     3     4     5     6     7     8     9    10    11    12    13
jd_dop = nav10(:,1);
Dspd = nav10(:,12);
Ddir = nav10(:,13);

%% load adcp data and compute surface current
adcp_path = 'D:\Data\CAPRICORN_2016\Investigator\flux\SCS\adcp\os150nb.nc';
jd_adcp = ncread(adcp_path,'time');
u = ncread(adcp_path,'u');  % +E current component, 50 depth bins
v = ncread(adcp_path,'v');  % +N current component, 50 depth bins
Nc = v(3,:)';               % select +N ~ 25m depth
Ec = u(3,:)';               % select +E ~ 25m depth
Nc(Nc>5) = NaN; Ec(Ec>5) = NaN;
Cspd = sqrt(Nc.^2 + Ec.^2); Cspd = Cspd(:);    % ADCP current speed, m/s
Cdir = mod(atan2(Ec,Nc)+2*pi, 2*pi)*180/pi;    % ADCP current direction

% average to doppler timebase
jdref = [jd_dop; jd_dop(end)+mean(diff(jd_dop))];
VU10 = interval_avg(jd_adcp,[Nc,Ec],jdref);
Cspd2 = sqrt(VU10(:,2).^2 + VU10(:,3).^2);
Cdir2 = mod(atan2(VU10(:,3),VU10(:,2))+2*pi, 2*pi)*180/pi;

%% plot currents for selected days
jj = find(jd_dop>=stjd & jd_dop<=endjd+1);

figure; plot(jd_dop(jj),Cspd2(jj),'r.',jd_dop(jj),Dspd(jj),'bo'); grid;
xlabel('Day of Year'); ylabel('sfc current speed, m/s');
legend('OS75, 25m depth','Doppler Recomputed');
title('ADCP & Doppler Spd Log Surface Current Speeds');
saveas(gcf,fullfile(fig_path,['ADCP_DOP_spd_',jd1,'_',jd2,'.fig']));

figure; plot(jd_dop(jj),Cdir2(jj),'r.',jd_dop(jj),Ddir(jj),'bo'); grid;
xlabel('Day of Year'); ylabel('sfc current direction'); ylim([0,360]);
legend('OS75, 25m depth','Doppler Recomputed');
title('ADCP & Doppler Spd Log Surface Current Direction');
saveas(gcf,fullfile(fig_path,['ADCP_DOP_dir_',jd1,'_',jd2,'.fig']));

kk = find(jd_dop>=stjd & jd_dop<=endjd+1 & Cspd2>0.5);

figure; plot([0,3],[0,3],'k-',Cspd2,Dspd,'ro'); grid;
axis([0,1.5,0,1.5]);
xlabel('OS75 Current Speed, 25m depth, m/s');
ylabel('Doppler Current Speed, 7m depth, m/s');
title('Speed comparison');
saveas(gcf,fullfile(fig_path,['ADCP_v_DOP_spd_',jd1,'_',jd2,'.fig']));

figure; plot([0,360],[0,360],'k-',Cdir2(kk),Ddir(kk),'ro'); grid;
axis([0,360,0,360]);
xlabel('OS75 Current Dir, 25m depth');
ylabel('Doppler Current Dir, 7m depth');
title('Direction Comparison (Spd > 0.5 m/s)');
saveas(gcf,fullfile(fig_path,['ADCP_v_DOP_dir_',jd1,'_',jd2,'.fig']));
