%{
 process_usonic.m  -  reads Investigator raw 1 Hz NMEA meteo files and
    reprocesses Gill 2D sonic data to m/s.
%}

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

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

% move to raw met directory
cd ('D:\Data\CAPRICORN_2016\Investigator\flux\SCS\met_nmea\');
files = dir('*.METEO');

udat = [];

for ii = 1:length(files)
    dfl = files(ii).name;
    disp(['processing ',dfl]);
    
    met = importdata(dfl);  % load as cell array
    
    % scan for sonic data, col 7-8 will be dir & spd kts
    raw = {}; praw = []; temp = [];
    temp = met(cellfun(@(x) ~isempty(strfind(x,'$WIMWV')),met));
    for jj=1:length(temp)
        raw =textscan(temp{jj,:},'%4f-%2f-%2f%*c%2f:%2f:%6.3f %*6c %f %*c %f %*[^\n]','delimiter',',');
        try
            praw(jj,:)=[raw{1},raw{2},raw{3},raw{4},raw{5},raw{6},raw{7}*1.0,raw{8}*0.514444]; % kts to m/s
        catch
            praw(jj,:)=NaN(1,8);
        end;
    end;
    praw(isnan(praw(:,1)),:) = [];
    jd_p = datenum(praw(:,1:6))-datenum(2015,12,31);
    praw = [jd_p,praw(:,7:8)];
        
    udat = [udat;praw];
end;

clear praw sraw temp raw met jd_p

%% convert to spd/dir to UV
[U,V] = sd_to_uv(udat(:,3),udat(:,2));
udat(:,2) = U;
udat(:,3) = V;

%% reference 5 sec timestamp
stjd = floor(udat(1,1));
endjd = ceil(udat(end,1));
jdref = stjd:5/86400:endjd; jdref = jdref(:); % 5-sec bin edges

% compute 5 sec mean of raw usonic data and save to single array
% all speeds are m/s
% cols: jd, U, V, rwspd, rwdir
%        1, 2, 3,   4,     5
uson5 = NaN(length(jdref)-1,5);
uson5(:,1) = jdref(1:end-1);
avg = interval_avg(udat(:,1),udat(:,2:3),jdref);
uson5(:,2:3) = avg(:,2:3);
[rwspd,rwdir] = uv_to_sd(uson5(:,2),uson5(:,3));
uson5(:,4:5) = [rwspd,rwdir];

clear udat U V jdref avg

%% load psd wind for selected day
%{
1 jd 1min
2 U
3 V
4 W
5 Tsonic
6 rel wind dir
7 rel wind spd
8 true wind dir
9 true wind spd
%}
psd = load('D:\Data\CAPRICORN_2016\Investigator\flux\Processed_In_Field\CAPRICORN_2016_proc_son_1min_2016_03_21_081.txt');
jdpsd = psd(:,1);

%% load ship 5 sec data for selected day
%{
1       jd                  decimal day of year
2       lat                 seapath latitude
3       lon                 seapath longitude, deg E
4       sog                 seapath sog, m/s
5       cog                 seapath cog, deg
6       hed                 gyro heading, deg
7       wspd_port           true wspd, port, m/s @ 25m
8       wdir_port           true wdir, port, deg
9       rwspd_port          rel wspd, port, m/s @ 25m
10      rwdir_port          rel wdir, port, 0-360 deg
11      wspd_stbd           true wspd, stbd, m/s @ 25m
12      wdir_stbd           true wdir, stbd, deg
13      rwspd_stbd          rel wspd, stbd, m/s @ 25m
14      rwdir_stbd          rel wdir, stbd, 0-360 deg
15      wspd_sonic          true wspd, sonic, m/s @ 22.5m
16      wdir_sonic          true wdir, sonic, deg
17      rwspd_sonic         rel wspd, sonic, m/s @ 22.5m
18      rwdir_sonic         rel wdir, sonic, 0-360 deg
19      ta_port             air temp, port, C @ 31.4m
20      ta_stbd             air temp, stbd, C 31.4m
21      rh_port             rel humidity, port, % 31.4m
22      rh_stbd             rel humidity, stbd, % 31.4m
23      Pmb                 atm preessure, mb 20m
24      rain                accumulated rain, RM Young, mm
25      rain_rate           rain rate, RM Young, mm/hr
26      rain_optical        accumulated rain, optical disdrometer, mm
27      rain_rate_optical   rain rate, optical disdrometer, mm/hr
28      PSP_port            solar rad flux, port, W/m2
29      PSP_stbd            solar rad flux, stbd, W/m2
30      PIR_port            IR rad flux, port, W/m2
31      PIR_stbd            IR rad flux, stbd, W/m2
32      SST_IR              ISAR infrared sea surface temp, C
33      SST                 TSG water temp, C @ depth of port drop keel
34      Sal                 TSG salinity, psu
35      uship_long          longitudinal speed through water, m/s
36      uship_trans         transverse speed through water, m/s
37      SST_depth           depth of SST on port drop keel, m
38      sog_n               N-S component of ship speed
39      sog_w               E-W component of ship speed
40      hed_n               N-S component of ship heading
41      hed_w               E-W component of ship heading
%}

path_scs_data = 'D:\Data\CAPRICORN_2016\Investigator\flux\SCS\in2016_v02uwy.nc';
scs = read_scs_day_CAPRICORN_2016(path_scs_data,81);
jdscs = scs(:,1);
rwdir5 = scs(:,18);

%% plot rwspd
figure; plot(uson5(:,1),uson5(:,4),'r.',jdpsd,psd(:,7),'g.',jdscs,scs(:,17),'b.');
grid; xlim([jdpsd(1),jdpsd(end)]); xlabel('day-of-year'); ylabel('RWSPD m/s');
title('Relative Wind Speeds'); legend('Gill 2D Raw','NOAA','Gill 2D 5sec','location','best');

% plot rwdir
rwdir(rwdir>180) = rwdir(rwdir>180) - 360;
rwdir5(rwdir5>180) = rwdir5(rwdir5>180) - 360;
figure; plot(uson5(:,1),rwdir,'r.',jdpsd,psd(:,6),'g.',jdscs,rwdir5,'b.');
grid; xlim([jdpsd(1),jdpsd(end)]); xlabel('day-of-year'); ylabel('RWDIR');
title('Relative Wind Direction'); legend('Gill 2D Raw','NOAA','Gill 2D 5sec','location','best');



