% buoy_comparison_WHOTS_2022.m

% run this after running flux_eval_WHOTS_2022.m
% also, convert hourly WHOTS data from txt to matlab table format
% and save in WHOI data directory

% June 2023 LB

%%  set up parameters, edit as desired and run
clear;
close('all');
fclose('all');

% date range: JD 204 to 208 for complete WHOTS cruise
firstday = 204 + datenum(2021,12,31);  % start datenum
lastday = 208 + datenum(2021,12,31);   % end datenum

plotit = true;  % for doing plots.

graphformat = '.png';  % select graphics format files
graphdevice = '-dpng'; % select graphic device
cruise = 'WHOTS_2022'; % Cruise name
ship = 'RV_Sette';    % Research vessel

% system specific path definitions
sysType = computer;
if strncmp(sysType,'MACI64',7)       % set Mac paths
    data_drive = '/Users/ludo/Documents/Work/';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','flux');
    path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed');
    path_whoi_data = fullfile(data_drive,cruise,ship,'flux','WHOI');
    path_images = fullfile(data_drive,cruise,ship,'flux','Processed_Images','buoy_comparison');
    mkdir(path_images);
elseif strncmp(sysType,'PCWIN',7)    % set PC Windows paths
    data_drive = 'C:\Data';
    path_prog = fullfile(data_drive,cruise,ship,'Scientific_analysis','programs','flux');
    path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed');
    path_whoi_data = fullfile(data_drive,cruise,ship,'flux','WHOI');
    path_images = fullfile(data_drive,cruise,ship,'flux','Processed_Images','buoy_comparison');
    mkdir(path_images);
end

% matlab script path
restoredefaultpath
cd(path_prog);
addpath(genpath(path_prog));
savepath;

%% load 10-min and hourly data
% this assumes you have already run flux_eval_WHOTS_2021.m and you have
% converted the hourly mooring data files to matlab table format.

load(fullfile(path_proc_data,'WHOTS_2022_PSL_flux_1hr_all.mat'),'flxhr');
load(fullfile(path_proc_data,'WHOTS_2022_PSL_flux_10min_all.mat'),'flx');

W_18_1=load('/Users/ludo/Documents/Work/WHOTS_2022/RV_Sette/WHOI/WHOTS-XVIII_MET_sys1QC.txt');
W_18_2=load('/Users/ludo/Documents/Work/WHOTS_2022/RV_Sette/WHOI/WHOTS-XVIII_MET_sys2QC.txt');
W_17_1=load('/Users/ludo/Documents/Work/WHOTS_2022/RV_Sette/WHOI/WHOTS-XVII_MET_sys1QC.txt');
W_17_2=load('/Users/ludo/Documents/Work/WHOTS_2022/RV_Sette/WHOI/WHOTS-XVII_MET_sys2QC.txt');

%load previous year mat for vars names
load('/Users/ludo/Documents/Work/WHOTS_2022/RV_Sette/WHOI/WHOTS_16_sys1.mat','W_16_1');
wvars = W_16_1.Properties.VariableNames;
clear W_16_1

W_18_1 = array2table(W_18_1,'VariableNames',wvars);
W_18_2 = array2table(W_18_2,'VariableNames',wvars);
W_17_1 = array2table(W_17_1,'VariableNames',wvars);
W_17_2 = array2table(W_17_2,'VariableNames',wvars);

%% run some offset corrections to PSL based on lab data
%adjust sea snake to TSG at night
% flxhr.tsea=flxhr.tsea-nanmedian(flxhr.tsea(flxhr.rs<10)-flxhr.tsg2(flxhr.rs<10));


%%  set up selection criteria

% set up timestamps in datenum format
flxhr.dnum = flxhr.jd + datenum(2021,12,31);
W_18_1.dnum = datenum(W_18_1.Year,W_18_1.mo,W_18_1.day,W_18_1.hr,0,0);
W_18_2.dnum = datenum(W_18_2.Year,W_18_2.mo,W_18_2.day,W_18_2.hr,0,0);
W_17_1.dnum = datenum(W_17_1.Year,W_17_1.mo,W_17_1.day,W_17_1.hr,0,0);
W_17_2.dnum = datenum(W_17_2.Year,W_17_2.mo,W_17_2.day,W_17_2.hr,0,0);

% selector for entire time at HOT site: Aug 27 0600 to Aug 31 1900
tt = find(flxhr.dnum>datenum(2022,7,24,8,0,0) & flxhr.dnum<datenum(2022,7,27,5,0,0));
flxhr = flxhr(tt,:);  % select out ship data for comparison period
N = height(flxhr);
dt = datetime(flxhr.dnum,'convertfrom','datenum'); % datetime formate timestamp for plots...

% select mooring data covering same time period as PSL measurements
% and copy over into the same timebase (some hours will be missing).
% getting rid of the data we don't need simplifies the comparison code to follow.

% first, create empty arrays for the mooring data, same length as PSL table
C = width(W_18_1);
w181 = NaN(N,C);
w181(:,C) = flxhr.dnum;
w182 = w181; w171 = w181; w172 = w181;

% now, find coincident hours for each system and copy data into arrays
for ii = 1:N
    jj = find(abs(W_18_1.dnum - w181(ii,C))<2/1440);
    if ~ isempty(jj)
        w181(ii,1:19) = W_18_1{jj,1:19};
    end
    jj = find(abs(W_18_2.dnum - w182(ii,C))<1/1440);
    if ~ isempty(jj)
        w182(ii,1:19) = W_18_2{jj,1:19};
    end
    jj = find(abs(W_17_1.dnum - w171(ii,C))<1/1440);
    if ~ isempty(jj)
        w171(ii,1:19) = W_17_1{jj,1:19};
    end
    jj = find(abs(W_17_2.dnum - w172(ii,C))<1/1440);
    if ~ isempty(jj)
        w172(ii,1:19) = W_17_2{jj,1:19};
    end
end

% convert arrays to tables
wvars = W_18_1.Properties.VariableNames;
w181 = array2table(w181,'VariableNames',wvars);
w182 = array2table(w182,'VariableNames',wvars);
w171 = array2table(w171,'VariableNames',wvars);
w172 = array2table(w172,'VariableNames',wvars);

% compute hourly distance from the ship to each mooring anchor locations
% WHOTS-17 anchor  22 46.002 N, 157 53.958 W
% WHOTS-18 anchor  22 40.021 N, 157 57.018 W

d18 = NaN(N,1); d17 = d18;
lat18 = 22 + 40.021/60; lon18 = -(157 + 57.018/60); % lat/lon WHOTS 18
lat17 = 22 + 46.002/60; lon17 = -(157 + 53.958/60); % lat/lon WHOTS 17
for ii = 1:N
    [~,d18(ii)] = legs([flxhr.lat(ii),lat18],[flxhr.lon(ii)-360,lon18]);
    [~,d17(ii)] = legs([flxhr.lat(ii),lat17],[flxhr.lon(ii)-360,lon17]);    
end

% selector for period closest to WHOTS17: Start: 0800 UTC 7/24/22 .   End:   1400 UTC 7/25/22
tt17 = find(flxhr.dnum>datenum(2022,7,24,8,0,0) & flxhr.dnum<datenum(2022,7,25,14,0,0));

% selector for period closest to WHOTS18:   Start: 0800 UTC 7/26/22 .   End:   0500 UTC 7/27/22
tt18 = find(flxhr.dnum>datenum(2022,7,26,8,0,0) & flxhr.dnum<datenum(2022,7,27,5,0,0));

% selector for both intercomparison periods
all = union(tt18,tt17);

%% met comparison plots

% Distance to Moorings, timeseries
figure('position',[300,300,800,600],'name','Distance_to_moorings');
plot(dt,d18,'.',dt,d17,'.','markersize',10); grid on; grid minor;
ylabel('Distance, N. Miles','fontsize',12);
title('WHOTS 2022, Distance, Ship to Mooring Anchor Location','fontsize',16);
legend('WHOTS-18','WHOTS-17','location','northeast','fontsize',12);
if plotit
    ppath = fullfile(path_images,['01_dist_to_mooring',graphformat]);
    print(graphdevice,ppath);
end

% SL Pressure
padj = 3*0.125; % adjustment of mooring P from 3m to sea level
figure('position',[300,300,800,700],'name','Pressure_all');
plot([1016,1022],[1016,1022],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Pmb(tt18),w181.BPR(tt18)+padj,'o',flxhr.Pmb(tt18),w182.BPR(tt18)+padj,'o',...
          flxhr.Pmb(all),w171.BPR(all)+padj,'o',flxhr.Pmb(all),w172.BPR(all)+padj,'o',...
          flxhr.Pmb(all),flxhr.Pmb_scs(all),'o','linewidth',1.5); hold off;
xlabel('Pressure, mb (PSL PTB)','fontsize',14);
ylabel('Pressure, mb (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Sea Level Pressure Comparison - both periods','fontsize',16);
% subtitle('For both intercomparison periods: Aug 27 - Aug 31, 2022','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2','SCS', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['02_pressure_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','Pressure_nearby');
plot([1016,1022],[1016,1022],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Pmb(tt18),w181.BPR(tt18)+padj,'o',flxhr.Pmb(tt18),w182.BPR(tt18)+padj,'o',...
          flxhr.Pmb(tt17),w171.BPR(tt17)+padj,'o',flxhr.Pmb(tt17),w172.BPR(tt17)+padj,'o',...
          'linewidth',1.5); hold off;
xlabel('Pressure, mb (PSL PTB)','fontsize',14);
ylabel('Pressure, mb (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Sea Level Pressure Comparison - near mooring','fontsize',16);
% subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['03_pressure_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Air Temperature (3m)
figure('position',[300,300,800,700],'name','Tair_all');
plot([24,26],[24,26],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Trf(tt18),w181.ATMP(tt18),'o',flxhr.Trf(tt18),w182.ATMP(tt18),'o',...
          flxhr.Trf(all),w171.ATMP(all),'o',flxhr.Trf(all),w172.ATMP(all),'o',...
          flxhr.Trf_scs(all),flxhr.ta_scs(all),'o','linewidth',1.5); hold off;
axis([24,26,24,26]);
xlabel('Air Temp, degC (COARE Ref)','fontsize',14);
ylabel('Air Temp, degC (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Air Temperature Comparison - both periods','fontsize',16);
% subtitle('PSL Ta adjusted to 3m for moorings and 21m for SCS','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2','SCS','location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['04_tair_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','Tair_nearby');
plot([24,26],[24,26],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Trf(tt18),w181.ATMP(tt18),'o',flxhr.Trf(tt18),w182.ATMP(tt18),'o',...
          flxhr.Trf(tt17),w171.ATMP(tt17),'o',flxhr.Trf(tt17),w172.ATMP(tt17),'o',...
          'linewidth',1.5); hold off;
axis([24,26,24,26]);
xlabel('Air Temp, degC (COARE 3m Ref)','fontsize',14);
ylabel('Air Temp, degC (moorings)','fontsize',14);
title('WHOTS 2022, Hourly 3m Air Temperature Comparison - near mooring','fontsize',16);
% subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['05_tair_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Relative Humidity (3m)
figure('position',[300,300,800,700],'name','RH_all');
plot([66,84],[66,84],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.RHrf(tt18),w181.RH(tt18),'o',flxhr.RHrf(tt18),w182.RH(tt18),'o',...
          flxhr.RHrf(all),w171.RH(all),'o',flxhr.RHrf(all),w172.RH(all),'o',...
          flxhr.RHrf_scs(all),flxhr.rh_scs(all),'o','linewidth',1.5); hold off;
      axis([66 84 66 84])
xlabel('RH %, (COARE Ref)','fontsize',14);
ylabel('RH %, (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Relative Humidity Comparison - both periods','fontsize',16);
% subtitle('PSL RH adjusted to 3m for moorings and 21m for SCS','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2','SCS','location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['06_rh_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','RH_nearby');
plot([68,86],[68,86],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.RHrf(tt18),w181.RH(tt18),'o',flxhr.RHrf(tt18),w182.RH(tt18),'o',...
          flxhr.RHrf(tt17),w171.RH(tt17),'o',flxhr.RHrf(tt17),w172.RH(tt17),'o',...
          'linewidth',1.5); hold off;
      axis([66 84 66 84])
xlabel('RH %, (COARE 3m Ref)','fontsize',14);
ylabel('RH %, (moorings)','fontsize',14);
title('WHOTS 2022, Hourly 3m Relative Humidity Comparison - near mooring','fontsize',16);
% subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['07_rh_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Shortwave Flux
figure('position',[300,300,800,700],'name','RS_all');
plot([-200,1200],[-200,1200],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.rs(tt18),w181.SRAD(tt18),'o',flxhr.rs(tt18),w182.SRAD(tt18),'o',...
          flxhr.rs(all),w171.SRAD(all),'o',flxhr.rs(all),w172.SRAD(all),'o',...
          'linewidth',1.5); hold off;
xlabel('Shortwave Radiation, W/m^2 (PSL K&Z)','fontsize',14);
ylabel('Shortwave Radiation, W/m^2 (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Shortwave Radiation Comparison - both periods','fontsize',16);
% subtitle('For both intercomparison periods: Aug 27 - Aug 31, 2022','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['08_rs_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','RS_nearby');
plot([-200,1200],[-200,1200],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.rs(tt18),w181.SRAD(tt18),'o',flxhr.rs(tt18),w182.SRAD(tt18),'o',...
          flxhr.rs(tt17),w171.SRAD(tt17),'o',flxhr.rs(tt17),w172.SRAD(tt17),'o',...
          'linewidth',1.5); hold off;
xlabel('Shortwave Radiation, W/m^2 (PSL K&Z)','fontsize',14);
ylabel('Shortwave Radiation, W/m^2 (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Shortwave Radiation Comparison - near mooring','fontsize',16);
% subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['09_rs_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Longwave Flux
figure('position',[300,300,800,700],'name','RL_all');
plot([360,440],[360,440],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.rl(tt18),w181.LRAD(tt18),'o',flxhr.rl(tt18),w182.LRAD(tt18),'o',...
          flxhr.rl(all),w171.LRAD(all),'o',flxhr.rl(all),w172.LRAD(all),'o',...
          'linewidth',1.5); hold off;
xlabel('Longwave Radiation, W/m^2 (PSL PIR)','fontsize',14);
ylabel('Longwave Radiation, W/m^2 (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Longwave Radiation Comparison - both periods','fontsize',16);
% subtitle('For both intercomparison periods: Aug 27 - Aug 31, 2022','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['10_rl_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','RL_nearby');
plot([360,440],[360,440],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.rl(tt18),w181.LRAD(tt18),'o',flxhr.rl(tt18),w182.LRAD(tt18),'o',...
          flxhr.rl(tt17),w171.LRAD(tt17),'o',flxhr.rl(tt17),w172.LRAD(tt17),'o',...
          'linewidth',1.5); hold off;
xlabel('Longwave Radiation, W/m^2 (PSL PIR)','fontsize',14);
ylabel('Longwave Radiation, W/m^2 (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Longwave Radiation Comparison - near mooring','fontsize',16);
% subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['11_rl_nearby',graphformat]);
    print(graphdevice,ppath);
end

% SST
figure('position',[300,300,800,700],'name','SST_all');
plot([25.5,26.2],[25.5,26.2],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.tsea(all),w181.STMP(all),'o',flxhr.tsea(all),w182.STMP(all),'o',...
          flxhr.tsea(all),w171.STMP(all),'o',flxhr.tsea(all),w172.STMP(all),'o',...
          flxhr.tsea(all),flxhr.tsg2(all),'o','linewidth',1.5); hold off;
xlabel('Sea Temp, degC, (PSL Seasnake)','fontsize',14);
ylabel('Sea Temp, degC (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Sea Surface Comparison - both periods','fontsize',16);
% subtitle('For both intercomparison periods: Aug 27 - Aug 31, 2022','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2','SCS','location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['12_sst_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','SST_nearby');
plot([25.5,26.2],[25.5,26.2],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.tsea(tt18),w181.STMP(tt18),'o',flxhr.tsea(tt18),w182.STMP(tt18),'o',...
          flxhr.tsea(tt17),w171.STMP(tt17),'o',flxhr.tsea(tt17),w172.STMP(tt17),'o',...
          'linewidth',1.5); hold off;
xlabel('Sea Temp, degC, (PSL Seasnake)','fontsize',14);
ylabel('Sea Temp, degC (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Sea Surface Temperature Comparison - near mooring','fontsize',16);
% subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['13_sst_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Wind Speed (3m)
figure('position',[300,300,800,700],'name', 'WSPD_all');
plot([0,10],[0,10],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Urf(tt18),w181.WSPD(tt18),'o',flxhr.Urf(tt18),w182.WSPD(tt18),'o',...
          flxhr.Urf(all),w171.WSPD(all),'o',flxhr.Urf(all),w172.WSPD(all),'o',...
          flxhr.Urf_scs(all),flxhr.wspd_scs(all),'o','linewidth',1.5); hold off;
xlabel('Wind Speed, m/s, (COARE Ref)','fontsize',14);
ylabel('Wind Speed, m/s, (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Wind Speed Comparison - both periods','fontsize',16);
% subtitle('PSL Wspd adjusted to 3m for moorings and 21m for SCS','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2','SCS','location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['14_wspd_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','WSPD_nearby');
plot([0,10],[0,10],'k-'); grid on; grid minor; hold on;
h1 = plot(flxhr.Urf(tt18),w181.WSPD(tt18),'o',flxhr.Urf(tt18),w182.WSPD(tt18),'o',...
          flxhr.Urf(tt17),w171.WSPD(tt17),'o',flxhr.Urf(tt17),w172.WSPD(tt17),'o',...
          'linewidth',1.5); hold off;
xlabel('Wind Speed, m/s, (COARE 3m Ref)','fontsize',14);
ylabel('Wind Speed, m/s, (moorings)','fontsize',14);
title('WHOTS 2022, Hourly 3m Wind Speed Comparison - near mooring','fontsize',16);
% subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['15_wspd_nearby',graphformat]);
    print(graphdevice,ppath);
end

% Wind Direction
wdir = flxhr.wdir + 180;   % to oceanographic convention 'direction to'
wdir(wdir>360) = wdir(wdir>360) - 360;
wdir(wdir<0) = wdir(wdir<0) + 360;
wdir_scs = flxhr.wdir_scs + 180;   % to oceanographic convention 'direction to'
wdir_scs(wdir_scs>360) = wdir_scs(wdir_scs>360) - 360;
wdir_scs(wdir_scs<0) = wdir_scs(wdir_scs<0) + 360;
figure('position',[300,300,800,700],'name','WDIR_all');
plot([0,360],[0,360],'k-'); grid on; grid minor; hold on;
h1 = plot(wdir(tt18),w181.WDIR(tt18),'o',wdir(tt18),w182.WDIR(tt18),'o',...
          wdir(all),w171.WDIR(all),'o',wdir(all),w172.WDIR(all),'o',...
          wdir(all),wdir_scs(all),'o','linewidth',1.5); hold off;
axis([200,340,200,340]);
xlabel('Wind Direction ''to'', deg (PSL WXT)','fontsize',14);
ylabel('Wind Direction ''to'', deg (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Wind Direction (Ocn Convention) - both periods','fontsize',16);
% subtitle('For both intercomparison periods: Aug 27 - Aug 31, 2022','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2','SCS','location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['16_wdir_all',graphformat]);
    print(graphdevice,ppath);
end

figure('position',[300,300,800,700],'name','WDIR_nearby');
plot([0,360],[0,360],'k-'); grid on; grid minor; hold on;
h1 = plot(wdir(tt18),w181.WDIR(tt18),'o',wdir(tt18),w182.WDIR(tt18),'o',...
          wdir(tt17),w171.WDIR(tt17),'o',wdir(tt17),w172.WDIR(tt17),'o',...
          'linewidth',1.5); hold off;
axis([200,340,200,340]);
xlabel('Wind Direction ''to'', deg (PSL WXT)','fontsize',14);
ylabel('Wind Direction ''to'', deg (moorings)','fontsize',14);
title('WHOTS 2022, Hourly Wind Direction (Ocn Convention) - near mooring','fontsize',16);
% subtitle('Only direct comparison hours near mooring','fontsize',14);
legend(h1,'w18-1','w18-2','W17-1','W17-2', 'location','northwest','fontsize',14);
if plotit
    ppath = fullfile(path_images,['17_wdir_nearby',graphformat]);
    print(graphdevice,ppath);
end

%% prepare vectors of stats for the ship and buoy comparison:
% comparison is PSL value – ship or mooring value

% for entire period at HOT site:
psl_scs_all_means  = [nanmean1(flxhr.Trf_scs(all) - flxhr.ta_scs(all)),...
                      nanmean1(flxhr.RHrf_scs(all) - flxhr.rh_scs(all)),...
                      nanmean1(flxhr.tsea(all) - flxhr.tsg2(all)),...
                      nanmean1(flxhr.Urf_scs(all) - flxhr.wspd_scs(all)),...
                      nanmean1(flxhr.wdir(all) - flxhr.wdir_scs(all)),...
                      NaN,... % no rs
                      NaN,... % no rl
                      nanmean1(flxhr.Pmb(all) - flxhr.Pmb_scs(all))];

psl_scs_all_meds  =  [nanmedian1(flxhr.Trf_scs(all) - flxhr.ta_scs(all)),...
                      nanmedian1(flxhr.RHrf_scs(all) - flxhr.rh_scs(all)),...
                      nanmedian1(flxhr.tsea(all) - flxhr.tsg2(all)),...
                      nanmedian1(flxhr.Urf_scs(all) - flxhr.wspd_scs(all)),...
                      nanmedian1(flxhr.wdir(all) - flxhr.wdir_scs(all)),...
                      NaN,...
                      NaN,...
                      nanmedian1(flxhr.Pmb(all) - flxhr.Pmb_scs(all))];

psl_scs_all_std  =   [nanstd1(flxhr.Trf_scs(all) - flxhr.ta_scs(all)),...
                      nanstd1(flxhr.RHrf_scs(all) - flxhr.rh_scs(all)),...
                      nanstd1(flxhr.tsea(all) - flxhr.tsg2(all)),...
                      nanstd1(flxhr.Urf_scs(all) - flxhr.wspd_scs(all)),...
                      nanstd1(flxhr.wdir(all) - flxhr.wdir_scs(all)),...
                      NaN,...
                      NaN,...
                      nanstd1(flxhr.Pmb(all) - flxhr.Pmb_scs(all))];

psl_w171_all_means = [nanmean1(flxhr.Trf(all) - w171.ATMP(all)),...
                      nanmean1(flxhr.RHrf(all) - w171.RH(all))...
                      nanmean1(flxhr.tsea(all) - w171.STMP(all))...
                      nanmean1(flxhr.Urf(all) - w171.WSPD(all))...
                      nanmean1(wdir(all) - w171.WDIR(all))...
                      nanmean1(flxhr.rs(all) - w171.SRAD(all))...
                      nanmean1(flxhr.rl(all) - w171.LRAD(all))...
                      nanmean1(flxhr.Pmb(all) - w171.BPR(all))];

psl_w171_all_meds =  [nanmedian1(flxhr.Trf(all) - w171.ATMP(all)),...
                      nanmedian1(flxhr.RHrf(all) - w171.RH(all))...
                      nanmedian1(flxhr.tsea(all) - w171.STMP(all))...
                      nanmedian1(flxhr.Urf(all) - w171.WSPD(all))...
                      nanmedian1(wdir(all) - w171.WDIR(all))...
                      nanmedian1(flxhr.rs(all) - w171.SRAD(all))...
                      nanmedian1(flxhr.rl(all) - w171.LRAD(all))...
                      nanmedian1(flxhr.Pmb(all) - w171.BPR(all))];

psl_w171_all_std =   [nanstd1(flxhr.Trf(all) - w171.ATMP(all)),...
                      nanstd1(flxhr.RHrf(all) - w171.RH(all))...
                      nanstd1(flxhr.tsea(all) - w171.STMP(all))...
                      nanstd1(flxhr.Urf(all) - w171.WSPD(all))...
                      nanstd1(wdir(all) - w171.WDIR(all))...
                      nanstd1(flxhr.rs(all) - w171.SRAD(all))...
                      nanstd1(flxhr.rl(all) - w171.LRAD(all))...
                      nanstd1(flxhr.Pmb(all) - w171.BPR(all))];

psl_w172_all_means = [nanmean1(flxhr.Trf(all) - w172.ATMP(all)),...
                      nanmean1(flxhr.RHrf(all) - w172.RH(all))...
                      nanmean1(flxhr.tsea(all) - w172.STMP(all))...
                      nanmean1(flxhr.Urf(all) - w172.WSPD(all))...
                      nanmean1(wdir(all) - w172.WDIR(all))...
                      nanmean1(flxhr.rs(all) - w172.SRAD(all))...
                      nanmean1(flxhr.rl(all) - w172.LRAD(all))...
                      nanmean1(flxhr.Pmb(all) - w172.BPR(all))];

psl_w172_all_meds =  [nanmedian1(flxhr.Trf(all) - w172.ATMP(all)),...
                      nanmedian1(flxhr.RHrf(all) - w172.RH(all))...
                      nanmedian1(flxhr.tsea(all) - w172.STMP(all))...
                      nanmedian1(flxhr.Urf(all) - w172.WSPD(all))...
                      nanmedian1(wdir(all) - w172.WDIR(all))...
                      nanmedian1(flxhr.rs(all) - w172.SRAD(all))...
                      nanmedian1(flxhr.rl(all) - w172.LRAD(all))...
                      nanmedian1(flxhr.Pmb(all) - w172.BPR(all))];

psl_w172_all_std =   [nanstd1(flxhr.Trf(all) - w182.ATMP(all)),...
                      nanstd1(flxhr.RHrf(all) - w182.RH(all))...
                      nanstd1(flxhr.tsea(all) - w182.STMP(all))...
                      nanstd1(flxhr.Urf(all) - w182.WSPD(all))...
                      nanstd1(wdir(all) - w182.WDIR(all))...
                      nanstd1(flxhr.rs(all) - w182.SRAD(all))...
                      nanstd1(flxhr.rl(all) - w182.LRAD(all))...
                      nanstd1(flxhr.Pmb(all) - w182.BPR(all))];

psl_w181_all_means = [nanmean1(flxhr.Trf(all) - w181.ATMP(all)),...
                      nanmean1(flxhr.RHrf(all) - w181.RH(all))...
                      nanmean1(flxhr.tsea(all) - w181.STMP(all))...
                      nanmean1(flxhr.Urf(all) - w181.WSPD(all))...
                      nanmean1(wdir(all) - w181.WDIR(all))...
                      nanmean1(flxhr.rs(all) - w181.SRAD(all))...
                      nanmean1(flxhr.rl(all) - w181.LRAD(all))...
                      nanmean1(flxhr.Pmb(all) - w181.BPR(all))];

psl_w181_all_meds =  [nanmedian1(flxhr.Trf(all) - w181.ATMP(all)),...
                      nanmedian1(flxhr.RHrf(all) - w181.RH(all))...
                      nanmedian1(flxhr.tsea(all) - w181.STMP(all))...
                      nanmedian1(flxhr.Urf(all) - w181.WSPD(all))...
                      nanmedian1(wdir(all) - w181.WDIR(all))...
                      nanmedian1(flxhr.rs(all) - w181.SRAD(all))...
                      nanmedian1(flxhr.rl(all) - w181.LRAD(all))...
                      nanmedian1(flxhr.Pmb(all) - w181.BPR(all))];

psl_w181_all_std =   [nanstd1(flxhr.Trf(all) - w181.ATMP(all)),...
                      nanstd1(flxhr.RHrf(all) - w181.RH(all))...
                      nanstd1(flxhr.tsea(all) - w181.STMP(all))...
                      nanstd1(flxhr.Urf(all) - w181.WSPD(all))...
                      nanstd1(wdir(all) - w181.WDIR(all))...
                      nanstd1(flxhr.rs(all) - w181.SRAD(all))...
                      nanstd1(flxhr.rl(all) - w181.LRAD(all))...
                      nanstd1(flxhr.Pmb(all) - w181.BPR(all))];

psl_w182_all_means = [nanmean1(flxhr.Trf(all) - w182.ATMP(all)),...
                      nanmean1(flxhr.RHrf(all) - w182.RH(all))...
                      nanmean1(flxhr.tsea(all) - w182.STMP(all))...
                      nanmean1(flxhr.Urf(all) - w182.WSPD(all))...
                      nanmean1(wdir(all) - w182.WDIR(all))...
                      nanmean1(flxhr.rs(all) - w182.SRAD(all))...
                      nanmean1(flxhr.rl(all) - w182.LRAD(all))...
                      nanmean1(flxhr.Pmb(all) - w182.BPR(all))];

psl_w182_all_meds =  [nanmedian1(flxhr.Trf(all) - w182.ATMP(all)),...
                      nanmedian1(flxhr.RHrf(all) - w182.RH(all))...
                      nanmedian1(flxhr.tsea(all) - w182.STMP(all))...
                      nanmedian1(flxhr.Urf(all) - w182.WSPD(all))...
                      nanmedian1(wdir(all) - w182.WDIR(all))...
                      nanmedian1(flxhr.rs(all) - w182.SRAD(all))...
                      nanmedian1(flxhr.rl(all) - w182.LRAD(all))...
                      nanmedian1(flxhr.Pmb(all) - w182.BPR(all))];

psl_w182_all_std =   [nanstd1(flxhr.Trf(all) - w182.ATMP(all)),...
                      nanstd1(flxhr.RHrf(all) - w182.RH(all))...
                      nanstd1(flxhr.tsea(all) - w182.STMP(all))...
                      nanstd1(flxhr.Urf(all) - w182.WSPD(all))...
                      nanstd1(wdir(all) - w182.WDIR(all))...
                      nanstd1(flxhr.rs(all) - w182.SRAD(all))...
                      nanstd1(flxhr.rl(all) - w182.LRAD(all))...
                      nanstd1(flxhr.Pmb(all) - w182.BPR(all))];
                  
                  
                  
% Now compute same for comparison intervals near each mooring
psl_w181_nearby_means = [nanmean1(flxhr.Trf(tt18) - w181.ATMP(tt18)),...
                      nanmean1(flxhr.RHrf(tt18) - w181.RH(tt18))...
                      nanmean1(flxhr.tsea(tt18) - w181.STMP(tt18))...
                      nanmean1(flxhr.Urf(tt18) - w181.WSPD(tt18))...
                      nanmean1(wdir(tt18) - w181.WDIR(tt18))...
                      nanmean1(flxhr.rs(tt18) - w181.SRAD(tt18))...
                      nanmean1(flxhr.rl(tt18) - w181.LRAD(tt18))...
                      nanmean1(flxhr.Pmb(tt18) - w181.BPR(tt18))];

psl_w181_nearby_meds =  [nanmedian1(flxhr.Trf(tt18) - w181.ATMP(tt18)),...
                      nanmedian1(flxhr.RHrf(tt18) - w181.RH(tt18))...
                      nanmedian1(flxhr.tsea(tt18) - w181.STMP(tt18))...
                      nanmedian1(flxhr.Urf(tt18) - w181.WSPD(tt18))...
                      nanmedian1(wdir(tt18) - w181.WDIR(tt18))...
                      nanmedian1(flxhr.rs(tt18) - w181.SRAD(tt18))...
                      nanmedian1(flxhr.rl(tt18) - w181.LRAD(tt18))...
                      nanmedian1(flxhr.Pmb(tt18) - w181.BPR(tt18))];

psl_w181_nearby_std =   [nanstd1(flxhr.Trf(tt18) - w181.ATMP(tt18)),...
                      nanstd1(flxhr.RHrf(tt18) - w181.RH(tt18))...
                      nanstd1(flxhr.tsea(tt18) - w181.STMP(tt18))...
                      nanstd1(flxhr.Urf(tt18) - w181.WSPD(tt18))...
                      nanstd1(wdir(tt18) - w181.WDIR(tt18))...
                      nanstd1(flxhr.rs(tt18) - w181.SRAD(tt18))...
                      nanstd1(flxhr.rl(tt18) - w181.LRAD(tt18))...
                      nanstd1(flxhr.Pmb(tt18) - w181.BPR(tt18))];

psl_w182_nearby_means = [nanmean1(flxhr.Trf(tt18) - w182.ATMP(tt18)),...
                      nanmean1(flxhr.RHrf(tt18) - w182.RH(tt18))...
                      nanmean1(flxhr.tsea(tt18) - w182.STMP(tt18))...
                      nanmean1(flxhr.Urf(tt18) - w182.WSPD(tt18))...
                      nanmean1(wdir(tt18) - w182.WDIR(tt18))...
                      nanmean1(flxhr.rs(tt18) - w182.SRAD(tt18))...
                      nanmean1(flxhr.rl(tt18) - w182.LRAD(tt18))...
                      nanmean1(flxhr.Pmb(tt18) - w182.BPR(tt18))];

psl_w182_nearby_meds =  [nanmedian1(flxhr.Trf(tt18) - w182.ATMP(tt18)),...
                      nanmedian1(flxhr.RHrf(tt18) - w182.RH(tt18))...
                      nanmedian1(flxhr.tsea(tt18) - w182.STMP(tt18))...
                      nanmedian1(flxhr.Urf(tt18) - w182.WSPD(tt18))...
                      nanmedian1(wdir(tt18) - w182.WDIR(tt18))...
                      nanmedian1(flxhr.rs(tt18) - w182.SRAD(tt18))...
                      nanmedian1(flxhr.rl(tt18) - w182.LRAD(tt18))...
                      nanmedian1(flxhr.Pmb(tt18) - w182.BPR(tt18))];

psl_w182_nearby_std =   [nanstd1(flxhr.Trf(tt18) - w182.ATMP(tt18)),...
                      nanstd1(flxhr.RHrf(tt18) - w182.RH(tt18))...
                      nanstd1(flxhr.tsea(tt18) - w182.STMP(tt18))...
                      nanstd1(flxhr.Urf(tt18) - w182.WSPD(tt18))...
                      nanstd1(wdir(tt18) - w182.WDIR(tt18))...
                      nanstd1(flxhr.rs(tt18) - w182.SRAD(tt18))...
                      nanstd1(flxhr.rl(tt18) - w182.LRAD(tt18))...
                      nanstd1(flxhr.Pmb(tt18) - w182.BPR(tt18))];

psl_w171_nearby_means = [nanmean1(flxhr.Trf(tt17) - w171.ATMP(tt17)),...
                      nanmean1(flxhr.RHrf(tt17) - w171.RH(tt17))...
                      nanmean1(flxhr.tsea(tt17) - w171.STMP(tt17))...
                      nanmean1(flxhr.Urf(tt17) - w171.WSPD(tt17))...
                      nanmean1(wdir(tt17) - w171.WDIR(tt17))...
                      nanmean1(flxhr.rs(tt17) - w171.SRAD(tt17))...
                      nanmean1(flxhr.rl(tt17) - w171.LRAD(tt17))...
                      nanmean1(flxhr.Pmb(tt17) - w171.BPR(tt17))];

psl_w171_nearby_meds =  [nanmedian1(flxhr.Trf(tt17) - w171.ATMP(tt17)),...
                      nanmedian1(flxhr.RHrf(tt17) - w171.RH(tt17))...
                      nanmedian1(flxhr.tsea(tt17) - w171.STMP(tt17))...
                      nanmedian1(flxhr.Urf(tt17) - w171.WSPD(tt17))...
                      nanmedian1(wdir(tt17) - w171.WDIR(tt17))...
                      nanmedian1(flxhr.rs(tt17) - w171.SRAD(tt17))...
                      nanmedian1(flxhr.rl(tt17) - w171.LRAD(tt17))...
                      nanmedian1(flxhr.Pmb(tt17) - w171.BPR(tt17))];

psl_w171_nearby_std =   [nanstd1(flxhr.Trf(tt17) - w171.ATMP(tt17)),...
                      nanstd1(flxhr.RHrf(tt17) - w171.RH(tt17))...
                      nanstd1(flxhr.tsea(tt17) - w171.STMP(tt17))...
                      nanstd1(flxhr.Urf(tt17) - w171.WSPD(tt17))...
                      nanstd1(wdir(tt17) - w171.WDIR(tt17))...
                      nanstd1(flxhr.rs(tt17) - w171.SRAD(tt17))...
                      nanstd1(flxhr.rl(tt17) - w171.LRAD(tt17))...
                      nanstd1(flxhr.Pmb(tt17) - w171.BPR(tt17))];

psl_w172_nearby_means = [nanmean1(flxhr.Trf(tt17) - w172.ATMP(tt17)),...
                      nanmean1(flxhr.RHrf(tt17) - w172.RH(tt17))...
                      nanmean1(flxhr.tsea(tt17) - w172.STMP(tt17))...
                      nanmean1(flxhr.Urf(tt17) - w172.WSPD(tt17))...
                      nanmean1(wdir(tt17) - w172.WDIR(tt17))...
                      nanmean1(flxhr.rs(tt17) - w172.SRAD(tt17))...
                      nanmean1(flxhr.rl(tt17) - w172.LRAD(tt17))...
                      nanmean1(flxhr.Pmb(tt17) - w172.BPR(tt17))];

psl_w172_nearby_meds =  [nanmedian1(flxhr.Trf(tt17) - w172.ATMP(tt17)),...
                      nanmedian1(flxhr.RHrf(tt17) - w172.RH(tt17))...
                      nanmedian1(flxhr.tsea(tt17) - w172.STMP(tt17))...
                      nanmedian1(flxhr.Urf(tt17) - w172.WSPD(tt17))...
                      nanmedian1(wdir(tt17) - w172.WDIR(tt17))...
                      nanmedian1(flxhr.rs(tt17) - w172.SRAD(tt17))...
                      nanmedian1(flxhr.rl(tt17) - w172.LRAD(tt17))...
                      nanmedian1(flxhr.Pmb(tt17) - w172.BPR(tt17))];

psl_w172_nearby_std =   [nanstd1(flxhr.Trf(tt17) - w172.ATMP(tt17)),...
                      nanstd1(flxhr.RHrf(tt17) - w172.RH(tt17))...
                      nanstd1(flxhr.tsea(tt17) - w172.STMP(tt17))...
                      nanstd1(flxhr.Urf(tt17) - w172.WSPD(tt17))...
                      nanstd1(wdir(tt17) - w172.WDIR(tt17))...
                      nanstd1(flxhr.rs(tt17) - w172.SRAD(tt17))...
                      nanstd1(flxhr.rl(tt17) - w172.LRAD(tt17))...
                      nanstd1(flxhr.Pmb(tt17) - w172.BPR(tt17))];

% save stats vector array as a txt file for import into spreadsheet
out = [psl_scs_all_means;
       psl_scs_all_meds;
       psl_scs_all_std;
       psl_w171_all_means;
       psl_w171_all_meds;
       psl_w171_all_std;
       psl_w172_all_means;
       psl_w172_all_meds;
       psl_w172_all_std;
       psl_w181_all_means;
       psl_w181_all_meds;
       psl_w181_all_std;
       psl_w182_all_means;
       psl_w182_all_meds;
       psl_w182_all_std;
       psl_w181_nearby_means;
       psl_w181_nearby_meds;
       psl_w181_nearby_std;
       psl_w182_nearby_means;
       psl_w182_nearby_meds;
       psl_w182_nearby_std;
       psl_w171_nearby_means;
       psl_w171_nearby_meds;
       psl_w171_nearby_std;
       psl_w172_nearby_means;
       psl_w172_nearby_meds;
       psl_w172_nearby_std];

sfl = fullfile(data_drive,cruise,ship,'Scientific_analysis','report','comp_table.txt');
save(sfl,'out','-ascii','-double','-tabs');

