%{
    process_nav.m
    
    reads Investigator seapath 1 Hz nav netCDF files, computes component
    vectors and saves to daily .mat files

    netCDF file name format: yyyymmdd-hhmmss-shipnav-shipnav_01.nav
%}

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

nav_path = 'D:\Data\CAPRICORN_2016\Investigator\flux\SCS\nav';
mat_path = 'D:\Data\CAPRICORN_2016\Investigator\flux\SCS\nav\mat';
mkdir(mat_path);
t0 = datenum(1899,12,31);  % zero time for netCDF timestamps
t00 = datenum(2015,12,31); % zero time for 2016 day-of-year values

cd(nav_path);
files = dir('*shipnav-shipnav_01.nav');

jd_sp = []; surfspeed = []; gndspeed = [];
surfcourse = []; gndcourse = [];

%% read all files, convert speeds to m/s

for ii = 1:length(files)
    fname = files(ii).name;
    disp(['reading netCDF file ',fname])
    temp = ncread(fname,'time') + t0 - t00;
    jd_sp = [jd_sp;temp];
    temp = ncread(fname,'surfspeed')*0.514444;
    surfspeed = [surfspeed;temp];
    temp = ncread(fname,'gndspeed')*0.514444;
    gndspeed = [gndspeed;temp];
    temp = ncread(fname,'surfcourse');
    surfcourse = [surfcourse;temp];
    temp = ncread(fname,'gndcourse');
    gndcourse = [gndcourse;temp];
end;

% remove bad surfspeed/surfcourse points
ii = surfspeed==0;
surfspeed(ii) = NaN; surfcourse(ii) = NaN;

%% compute surface speec components
Nsfc = surfspeed .* cos(surfcourse*pi/180);
Esfc = surfspeed .* sin(surfcourse*pi/180);

%% filter noise from gps ground spd/dir & compute components
[Gspd,Gdir,Ngnd,Egnd] = filt_spd_dir(gndspeed,gndcourse,1);

figure; plot(jd_sp,gndspeed,'r.-',jd_sp,Gspd,'b-'); grid;
title('Ground Speed, m/s');
figure; plot(jd_sp,gndcourse,'r.-',jd_sp,Gdir,'b-'); grid;
title('Ground Course'); ylim([0,360]);
%% compute current (drift) spd/dir

Ndft = Ngnd - Nsfc;
Edft = Egnd - Esfc;
Ndft = filt_lowpass(Ndft,1,0.04); % despike & filter current components
Edft = filt_lowpass(Edft,1,0.04);
Ndft(ii) = NaN; Edft(ii) = NaN;   % remove bad surfspeed points
Cspd = sqrt(Ndft.^2 + Edft.^2);
Cdir = mod(atan2(Edft,Ndft)+2*pi, 2*pi)*180/pi;

figure; plot(jd_sp,Cspd,'b-'); grid;
title('Current Speed, m/s');
figure; plot(jd_sp,Cdir,'b-'); grid;
title('Current Direction'); ylim([0,360]);

%% save 1 Hz components as daily .mat files
clear Gspd Gdir Cspd Cdir surfspeed surfcourse gndspeed gndcourse ii
nav0 = [Ngnd,Egnd,Nsfc,Esfc,Ndft,Edft];
stjd = floor(jd_sp(1));
endjd = floor(jd_sp(end));
for ii = stjd:endjd
    jdref = ii:1/86400:ii+86399/86400; jdref = jdref(:);
    jd = sprintf('%03i',floor(jdref(1)));
    nav = NaN(86400,7);
    nav(:,1) = jdref;
    nav(:,2:7) = interp1(jd_sp,nav0,jdref,'linear');
    fout = fullfile(mat_path,['nav_',jd,'.mat']);
    disp(['saving nav file for jd ',jd])
    save(fout,'nav');
end;

