% STRATUS_2019_buoy_comparison.m
% NOT REWRITTEN FOR STRATUS YET...
% NEED BUOY DATA FILE FORMAT AND SHIP DATA FORMAT IF AVAILABLE

% loads data & runs comparison plots for buoy vs ship data

%% initialize
%clear;
%close all;

% matlab script path
restoredefaultpath; userpath reset
%path_prog = fullfile('D:\Data\STRATUS_2019\CabodeHornos\flux\buoy_comparison');
path_prog = fullfile('D:\Data\STRATUS_2019\CabodeHornos\Scientific_analysis\programs\buoy_comparison');
addpath(genpath(path_prog));	% add scripts to path
savepath; rehash toolboxcache

%% define data paths & constants
path_data = 'D:\Data\STRATUS_2019\CabodeHornos\flux\processed\';
%%path_samos = '/Users/byronb/Raw_Data_Files.noindex/STRATUS_2019/CaboDeHornos/samos/netcdf';
path_whoi = 'D:\Data\STRATUS_2019\SCS\S18';
path_plots = 'D:\Data\STRATUS_2019\CabodeHornos\flux\buoy_comparison\';
mkdir(path_plots);
path_report = 'D:\Data\STRATUS_2019\CabodeHornos\Scientific_analysis\report\';
mkdir(path_report);

% start and end dates for the buoy intercomparison
%  206:215.125 for whole WHOTS cruise
%  208.0:215.125 at station ALOHA
stjd = 110;
endjd = 111;
Ndays = endjd-stjd+1;
jd_ref_1 = stjd:1/1440:endjd+1-1/1440; jd_ref_1 = jd_ref_1';
jd_ref_hr = stjd:1/24:endjd+1-1/24; jd_ref_hr = jd_ref_hr';
jd_ref_hr_bin = [jd_ref_hr;jd_ref_hr(end)+1/24];

graphformat = '.png';
graphdevice = '-dpng';
plotit = true;         % to run plots
printit = true;        % to save plots

% sensor heights
% WXT height (wind,P,T,RH,rain), top of bow tower
zwxt = 15.4;
% PSD pressure sensor height (m), top of wheelhouse
zp = 16.4;
% SCS pressure sensor height (m),
zpim = 24;
% SCS sonic height (m), top of bridge
zuim = 15.88;
% depth of ship tsgm, tssm
dtsg = 5;

%% load 1 hour met/flux data
%{
    flux1hr is nnn x 33
    COL   ID           DESCRIPTION
    1     jd_1min      Decimal julian day at beginning of line
    2     wspd         wxt true wind speed at zu, m/s
    3     wdir         wxt true wind direction, deg
    4     ts           sst seasnake T, C
    5     nd           NO DATA ship theromsalinograph T, C bow
    6     nd           NO DATA ship theromsalinograph salinity, psu bow
    7     ta           wxt air T at zwxt, C
    8     qs           air specific humidity at sea surface, g/kg
    9     qa           air specific humidity, g/kg
    10    rs           solar flux, w/m^2
    11    rl           IR flux, w/m^2
    12    rain         wxt precip rate, mm/hr
    13    cog          cog from psd gps
    14    sog          sog from psd gps, m/s
    15    hed          heading from psd gps, deg
    16    rwspd        wxt rel wind speed, m/s
    17    rwdir        wxt rel wind direction, deg
    18    lat          decimal latitude, deg
    19    lon          decimal longtude, deg
    20    zts          Depth of SST sensor used in heat flux calc, m
    21    usr          u-star, COARE 3.5
    22    taub         wind stress, COARE 3.5, N/m^2
    23    hsb          sensible heat flux, COARE 3.5, w/m^2
    24    hlb          latent heat flux, COARE 3.5, w/m^2
    25    rf           rain heat flux, COARE 3.5, w/m^2
    26    Urf          PSD wspd at 3m, COARE 3.5
    27    Trf          PSD ta at 3m, COARE 3.5
    28    RHrf         PSD rh at 3m, COARE 3.5
    29    Pmb_sl       wxt sea level pressure
    30    rh           wxt RH at height zq
    31    rslcr        relative wind, U vector, m/s
    32    rlclr        relative wind, V vector, m/s
    33    rwdir_std    N component of ship's speed, m/s
%}

% create empty 1-hr array for entire day range
flux1hr = NaN(Ndays*24,33);
flux1hr(:,1) = jd_ref_hr;

% load hourly data
%STRATUS_2019_PSD_flux_1hr_2019_04_14_104
%STRATUS_2019_PSD_flux_1hr_all
temp = load(fullfile(path_data,'STRATUS_2019_PSD_flux_1hr_all.txt'));

% copy into flux1hr array
for jj = 1:length(temp) 
    diff = jd_ref_hr - temp(jj,1);
    [val,idx] = min(abs(diff));
    if val<double(2.0/86400) && length(idx)==1
        flux1hr(idx,2:33) = temp(jj,2:33);
    end
end
clear temp

%% load buoy data

% load STRATUS-16 1-hr data: converted to N x 14 array:
%  1    2    3    4   5    6    7    8    9    10   11   12    13   14  15
%  yday BPR  ATMP SST WNDE WNDN WSPD WDIR RH   SRAD LRAD PRATE SAL  COND 
%  yday BPR  ATMP SST LWR  SWR  WNDN WNDE WSPD WDIR PRC  HRH   COND SAL SH
%  UTC  mBar C    C   w/m2 w/m2 m/s  m/s  m/s  deg  mm/hr %    S/m  ppm

% DEFINE EMPTY ARRAY FOR W13 LOGGER 1 and 2
W13_1_hr = NaN(length(flux1hr),33);
W13_1_hr(:,1) = jd_ref_hr(:,1);
W13_2_hr = NaN(length(flux1hr*24),33);
W13_2_hr(:,1) = jd_ref_hr(:,1);

load('D:\Data\STRATUS_2019\CabodeHornos\SCS\S18\W13_2_hr.mat');
%W13_2_hr = load(fullfile(path_whoi,'WHOTS13sys2_nohdr.txt'));
%W13_2_hr(:,5) = datenum(W13_2_hr(:,1),W13_2_hr(:,2),W13_2_hr(:,3),W13_2_hr(:,4),0,0);
%W13_2_hr = [W13_2_hr(:,5:15),W13_2_hr(:,17:19)];
%ii = mday(1,1)>=stjd+datenum(2017,12,31) & mday(1,1)<endjd+datenum(2017,12,31)+1;
%W13_2_hr = W13_2_hr(ii,:);              % select cruise days
W13_2_hr(:,8) = wnde(:,1)+180;      % convert wdir to 'from' convention
W13_2_hr(:,2) = bpr(:,1)+0.125*3;  % sea level pressure
W13_2_hr(:,8) = mod(wndn(:,1),360); % wdir
W13_2_hr(:,1) = mday(:,1) - datenum(2016,12,31);  % convert datenum to jd
 
%pad end of array to endjd
ii = find(abs((jd_ref_hr)-(W13_2_hr(end,1)))<2/86400);
N = length(jd_ref_hr) - ii;
W13_2_hr = [W13_2_hr;NaN(N,14)];
W13_2_hr(ii+1:end,1) = jd_ref_hr(ii+1:end);
%}
% load STRATUS-18 1-hr data: converted to N x 14 array:
%     1    2     3    4   5    6    7    8   9   10   11  12   13   14
%    yday BPR   ATMP SST WNDE WNDN WSPD WDIR RH SRAD LRAD PRATE SAL COND
%    UTC  mBar   C    C  m/s  m/s  m/s  deg  %  w/m2 w/m2 mm  mm/hr ppt  S/m
load('D:\Data\STRATUS_2019\CabodeHornos\SCS\S18\W13_1_hr.mat');
%load(fullfile(path_whoi,'W13_1_hr.mat'));
W13_1_hr(:,8) = wnde(:,1)+180;      % convert wdir to 'from' convention
W13_1_hr(:,2) = bpr(:,1)+0.125*3;  % sea level pressure
W13_1_hr(:,8) = mod(wndn(:,1),360); % wdir

load('D:\Data\STRATUS_2019\CabodeHornos\SCS\S18\W13_2_hr.mat');
%load(fullfile(path_whoi,'W13_2_hr.mat'));
W13_2_hr(:,8) = wnde(:,1)+180;      % convert wdir to 'from' convention
W13_2_hr(:,2) = bpr(:,1)+0.125*3;  % sea level pressure
W13_2_hr(:,8) = mod(wndn(:,1),360); % wdir

%% define variables
% PSD
lat1Hr = flux1hr(:,18);
lon1Hr = flux1hr(:,19);
sog1Hr = flux1hr(:,14);
rwdir1hr = flux1hr(:,17);
rwdir1hr_std = flux1hr(:,33);
spd1hr = flux1hr(:,2);
spd1hr_3m = flux1hr(:,26);
dir1hr = flux1hr(:,3);
ta1hr = flux1hr(:,7);
ta1hr_3m = flux1hr(:,27);
ts1hr = flux1hr(:,4);
rh1hr = flux1hr(:,30);
rh1hr_3m = flux1hr(:,28);
rs1hr = flux1hr(:,10);
rl1hr = flux1hr(:,11);
Pmb1hr = flux1hr(:,29);
qa1hr = flux1hr(:,9);
qa1hr_3m = qair_p([ta1hr_3m,rh1hr_3m,Pmb1hr]);

% ship

imspd1hr = scsm_hr(:,9);
imdir1hr = scsm_hr(:,10);
imta1hr = scsm_hr(:,8);
imts1hr = scsm_hr(:,14);
imtssg1hr = scsm_hr(:,16);
imrh1hr = scsm_hr(:,7);
imPmb1hr = scsm_hr(:,11);
imqa1hr = qair_p([imta1hr,imrh1hr,imPmb1hr]);

% W13
w13spd1hr_1 = W13_1_hr(:,7);
w13dir1hr_1 = W13_1_hr(:,8);
w13ta1hr_1 = W13_1_hr(:,3);
w13ts1hr_1 = W13_1_hr(:,4);
w13rh1hr_1 = W13_1_hr(:,9);
w13rs1hr_1 = W13_1_hr(:,10);
w13rl1hr_1 = W13_1_hr(:,11);
w13Pmb1hr_1 = W13_1_hr(:,2);
w13qa1hr_1 = qair_p([w13ta1hr_1,w13rh1hr_1,w13Pmb1hr_1]);

w13spd1hr_2 = W13_1_hr(:,7);
w13dir1hr_2 = W13_1_hr(:,8);
w13ta1hr_2 = W13_1_hr(:,3);
w13ts1hr_2 = W13_1_hr(:,4);
w13rh1hr_2 = W13_1_hr(:,9);
w13rs1hr_2 = W13_1_hr(:,10);
w13rl1hr_2 = W13_1_hr(:,11);
w13Pmb1hr_2 = W13_1_hr(:,2);
w13qa1hr_2 = qair_p([w13ta1hr_2,w13rh1hr_2,w13Pmb1hr_2]);

% W14
%{
w14spd1hr_1 = W14_1_hr(:,7);
w14dir1hr_1 = W14_1_hr(:,8);
w14ta1hr_1 = W14_1_hr(:,3);
w14ts1hr_1 = W14_1_hr(:,4);
w14rh1hr_1 = W14_1_hr(:,9);
w14rs1hr_1 = W14_1_hr(:,10);
w14rl1hr_1 = W14_1_hr(:,11);
w14Pmb1hr_1 = W14_1_hr(:,2);
w14qa1hr_1 = qair_p([w14ta1hr_1,w14rh1hr_1,w14Pmb1hr_1]);

w14spd1hr_2 = W13_1_hr(:,7);
w14dir1hr_2 = W13_1_hr(:,8);
w14ta1hr_2 = W13_1_hr(:,3);
w14ts1hr_2 = W13_1_hr(:,4);
w14rh1hr_2 = W13_1_hr(:,9);
w14rs1hr_2 = W13_1_hr(:,10);
w14rl1hr_2 = W13_1_hr(:,11);
w14Pmb1hr_2 = W13_1_hr(:,2);
w14qa1hr_2 = qair_p([w14ta1hr_2,w14rh1hr_2,w14Pmb1hr_2]);
%}
%% define selectors

% define mooring locations and position tolerance

lat13 = 24.787; lon13 = 269.093;
lat18 = 22.6699; lon14 = 275.3859;
tol = 0.03; % ~ 1.8 N miles or 3.33 km

% near STRATUS-16 lat/lon selector:
aa = find(lat1Hr>=lat13-tol & lat1Hr<=lat13+tol & lon1Hr>=lon13-tol & lon1Hr<=lon13+tol);

% near STRATUS-18 lat/lon selector:
bb = find(lat1Hr>=lat14-tol & lat1Hr<=lat14+tol & lon1Hr>=lon14-tol & lon1Hr<=lon14+tol);

% relative wind selector for ship measurements
cc = find((rwdir1hr>-60 & rwdir1hr<60) & rwdir1hr_std<30);

% filter for STRATUS-16 - use position and rel wind selectors only
ff = intersect(aa,cc);
% filter for STRATUS-18 - use position and rel wind selectors only
gg = intersect(bb,cc);

% selectors for entire period at station ALOHA:
dd = find(jd_ref_hr>208.0 & jd_ref_hr<215.125); % all data
hh = intersect(dd,cc);  % with rel wind restriction

%% build data arrays of [Ta,RH,SST,U,dir,Rs,Rl,P] for each platform
none = NaN(size(jd_ref_hr));
psd     = [ta1hr,rh1hr,ts1hr,spd1hr,unwrap(dir1hr*pi/180),rs1hr,rl1hr,Pmb1hr];
psd_3m  = [ta1hr_3m,rh1hr_3m,ts1hr,spd1hr_3m,unwrap(dir1hr*pi/180),rs1hr,rl1hr,Pmb1hr];
ship    = [imta1hr,imrh1hr,imts1hr,imspd1hr,unwrap(imdir1hr*pi/180),none,none,imPmb1hr];

w13_1 = [w13ta1hr_1,w13rh1hr_1,w13ts1hr_1,w13spd1hr_1,unwrap(w13dir1hr_1*pi/180),...
       w13rs1hr_1,w13rl1hr_1,w13Pmb1hr_1];
w13_2 = [w13ta1hr_2,w13rh1hr_2,w13ts1hr_2,w13spd1hr_2,unwrap(w13dir1hr_2*pi/180),...
       w13rs1hr_2,w13rl1hr_2,w13Pmb1hr_2];
%w14_1 = [w14ta1hr_1,w14rh1hr_1,w14ts1hr_1,w14spd1hr_1,unwrap(w14dir1hr_1*pi/180),...
 %      w14rs1hr_1,w14rl1hr_1,w14Pmb1hr_1];
%w14_2 = [w14ta1hr_2,w14rh1hr_2,w14ts1hr_2,w14spd1hr_2,unwrap(w14dir1hr_2*pi/180),...
  %     w14rs1hr_2,w14rl1hr_2,w14Pmb1hr_2];

%% psd vs ship stats for entire period at site

jj = [1,2,4];  % columns for Ta, RH, wspd
mean_psd_ship_all = nanmean1(psd(dd,:) - ship(dd,:));           % for Ts,Rs,Rl,P
mean_psd_ship_all(jj) = nanmean1(psd(hh,jj) - ship(hh,jj));     % for Ta,RH,wspd
mean_psd_ship_all(5) = nanmean1(psd(hh,5) - ship(hh,5))*180/pi; % wdir
if mean_psd_ship_all(5)>180; mean_psd_ship_all(5) = mean_psd_ship_all(5)-360; end
med_psd_ship_all = nanmedian1(psd(dd,:) - ship(dd,:));
med_psd_ship_all(jj) = nanmedian1(psd(hh,jj) - ship(hh,jj));
med_psd_ship_all(5) = nanmedian1(psd(hh,5) - ship(hh,5))*180/pi;
if med_psd_ship_all(5)>180; med_psd_ship_all(5) = med_psd_ship_all(5)-360; end
std_psd_ship_all = nanstd1(psd(dd,:) - ship(dd,:));
std_psd_ship_all(jj) = nanstd1(psd(hh,jj) - ship(hh,jj));
std_psd_ship_all(5) = nanstd1(psd(hh,5) - ship(hh,5))*180/pi;

% psd vs STRATUS-16 stats for entire period at site
mean_psd_w13_1_all = nanmean1(psd(dd,:) - w13_1(dd,:));
mean_psd_w13_1_all(jj) = nanmean1(psd(hh,jj) - w13_1(hh,jj));
mean_psd_w13_1_all(5) = nanmean1(psd(hh,5) - w13_1(hh,5))*180/pi;
med_psd_w13_1_all = nanmedian1(psd(dd,:) - w13_1(dd,:));
med_psd_w13_1_all(jj) = nanmedian1(psd(hh,jj) - w13_1(hh,jj));
med_psd_w13_1_all(5) = nanmedian1(psd(hh,5) - w13_1(hh,5))*180/pi;
std_psd_w13_1_all = nanstd1(psd(dd,:) - w13_1(dd,:));
std_psd_w13_1_all(jj) = nanstd1(psd(hh,jj) - w13_1(hh,jj));
std_psd_w13_1_all(5) = nanstd1(psd(hh,5) - w13_1(hh,5))*180/pi;

mean_psd_w13_2_all = nanmean1(psd(dd,:) - w13_2(dd,:));
mean_psd_w13_2_all(jj) = nanmean1(psd(hh,jj) - w13_2(hh,jj));
mean_psd_w13_2_all(5) = nanmean1(psd(hh,5) - w13_2(hh,5))*180/pi;
med_psd_w13_2_all = nanmedian1(psd(dd,:) - w13_2(dd,:));
med_psd_w13_2_all(jj) = nanmedian1(psd(hh,jj) - w13_2(hh,jj));
med_psd_w13_2_all(5) = nanmedian1(psd(hh,5) - w13_2(hh,5))*180/pi;
std_psd_w13_2_all = nanstd1(psd(dd,:) - w13_2(dd,:));
std_psd_w13_2_all(jj) = nanstd1(psd(hh,jj) - w13_2(hh,jj));
std_psd_w13_2_all(5) = nanstd1(psd(hh,5) - w13_2(hh,5))*180/pi;

% psd vs STRATUS-16 stats for period near bouy
mean_psd_w13_1 = nanmean1(psd(aa,:) - w13_1(aa,:));
mean_psd_w13_1(jj) = nanmean1(psd(ff,jj) - w13_1(ff,jj));
mean_psd_w13_1(5) = nanmean1(psd(ff,5) - w13_1(ff,5))*180/pi;
med_psd_w13_1 = nanmedian1(psd(aa,:) - w13_1(aa,:));
med_psd_w13_1(jj) = nanmedian1(psd(ff,jj) - w13_1(ff,jj));
med_psd_w13_1(5) = nanmedian1(psd(ff,5) - w13_1(ff,5))*180/pi;
std_psd_w13_1 = nanstd1(psd(aa,:) - w13_1(aa,:));
std_psd_w13_1(jj) = nanstd1(psd(ff,jj) - w13_1(ff,jj));
std_psd_w13_1(5) = nanstd1(psd(ff,5) - w13_1(ff,5))*180/pi;

mean_psd_w13_2 = nanmean1(psd(aa,:) - w13_2(aa,:));
mean_psd_w13_2(jj) = nanmean1(psd(ff,jj) - w13_2(ff,jj));
mean_psd_w13_2(5) = nanmean1(psd(ff,5) - w13_2(ff,5))*180/pi;
med_psd_w13_2 = nanmedian1(psd(aa,:) - w13_2(aa,:));
med_psd_w13_2(jj) = nanmedian1(psd(ff,jj) - w13_2(ff,jj));
med_psd_w13_2(5) = nanmedian1(psd(ff,5) - w13_2(ff,5))*180/pi;
std_psd_w13_2 = nanstd1(psd(aa,:) - w13_2(aa,:));
std_psd_w13_2(jj) = nanstd1(psd(ff,jj) - w13_2(ff,jj));
std_psd_w13_2(5) = nanstd1(psd(ff,5) - w13_2(ff,5))*180/pi;

% psd vs STRATUS-18 stats for entire period at site
%{
mean_psd_w14_1_all = nanmean1(psd(dd,:) - w14_1(dd,:));
mean_psd_w14_1_all(jj) = nanmean1(psd(hh,jj) - w14_1(hh,jj));
mean_psd_w14_1_all(5) = nanmean1(psd(hh,5) - w14_1(hh,5))*180/pi;
med_psd_w14_1_all = nanmedian1(psd(dd,:) - w14_1(dd,:));
med_psd_w14_1_all(jj) = nanmedian1(psd(hh,jj) - w14_1(hh,jj));
med_psd_w14_1_all(5) = nanmedian1(psd(hh,5) - w14_1(hh,5))*180/pi;
std_psd_w14_1_all = nanstd1(psd(dd,:) - w14_1(dd,:));
std_psd_w14_1_all(jj) = nanstd1(psd(hh,jj) - w14_1(hh,jj));
std_psd_w14_1_all(5) = nanstd1(psd(hh,5) - w14_1(hh,5))*180/pi;

mean_psd_w14_2_all = nanmean1(psd(dd,:) - w14_2(dd,:));
mean_psd_w14_2_all(jj) = nanmean1(psd(hh,jj) - w14_2(hh,jj));
mean_psd_w14_2_all(5) = nanmean1(psd(hh,5) - w14_2(hh,5))*180/pi;
med_psd_w14_2_all = nanmedian1(psd(dd,:) - w14_2(dd,:));
med_psd_w14_2_all(jj) = nanmedian1(psd(hh,jj) - w14_2(hh,jj));
med_psd_w14_2_all(5) = nanmedian1(psd(hh,5) - w14_2(hh,5))*180/pi;
std_psd_w14_2_all = nanstd1(psd(dd,:) - w14_2(dd,:));
std_psd_w14_2_all(jj) = nanstd1(psd(hh,jj) - w14_2(hh,jj));
std_psd_w14_2_all(5) = nanstd1(psd(hh,5) - w14_2(hh,5))*180/pi;

% psd vs STRATUS-18 stats for period near bouy
mean_psd_w14_1 = nanmean1(psd(aa,:) - w14_1(aa,:));
mean_psd_w14_1(jj) = nanmean1(psd(gg,jj) - w14_1(gg,jj));
mean_psd_w14_1(5) = nanmean1(psd(gg,5) - w14_1(gg,5))*180/pi;
med_psd_w14_1 = nanmedian1(psd(aa,:) - w14_1(aa,:));
med_psd_w14_1(jj) = nanmedian1(psd(gg,jj) - w14_1(gg,jj));
med_psd_w14_1(5) = nanmedian1(psd(gg,5) - w14_1(gg,5))*180/pi;
std_psd_w14_1 = nanstd1(psd(aa,:) - w14_1(aa,:));
std_psd_w14_1(jj) = nanstd1(psd(gg,jj) - w14_1(gg,jj));
std_psd_w14_1(5) = nanstd1(psd(gg,5) - w14_1(gg,5))*180/pi;

mean_psd_w14_2 = nanmean1(psd(aa,:) - w14_2(aa,:));
mean_psd_w14_2(jj) = nanmean1(psd(gg,jj) - w14_2(gg,jj));
mean_psd_w14_2(5) = nanmean1(psd(gg,5) - w14_2(gg,5))*180/pi;
med_psd_w14_2 = nanmedian1(psd(aa,:) - w14_2(aa,:));
med_psd_w14_2(jj) = nanmedian1(psd(gg,jj) - w14_2(gg,jj));
med_psd_w14_2(5) = nanmedian1(psd(gg,5) - w14_2(gg,5))*180/pi;
std_psd_w14_2 = nanstd1(psd(aa,:) - w14_2(aa,:));
std_psd_w14_2(jj) = nanstd1(psd(gg,jj) - w14_2(gg,jj));
std_psd_w14_2(5) = nanstd1(psd(gg,5) - w14_2(gg,5))*180/pi;
%}
%% save stats to table format for export to table
% 8 columns: Ta  RH  SST  U  dir  Rs  Rl  P

Stats = [];
Stats = [mean_psd_ship_all;med_psd_ship_all;std_psd_ship_all];          % psd - ship for all
Stats = [Stats;mean_psd_w13_1_all;med_psd_w13_1_all;std_psd_w13_1_all]; % psd - w13-1 for all
Stats = [Stats;mean_psd_w13_2_all;med_psd_w13_2_all;std_psd_w13_2_all]; % psd - w13-2 for all
%Stats = [Stats;mean_psd_w14_1_all;med_psd_w14_1_all;std_psd_w14_1_all]; % psd - w14-1 for all
%Stats = [Stats;mean_psd_w14_2_all;med_psd_w14_2_all;std_psd_w14_2_all]; % psd - w14-2 for all
Stats = [Stats;mean_psd_w13_1;med_psd_w13_1;std_psd_w13_1];             % psd - w13-1 near bouy
Stats = [Stats;mean_psd_w13_2;med_psd_w13_2;std_psd_w13_2];             % psd - w13-2 near bouy
%Stats = [Stats;mean_psd_w14_1;med_psd_w14_1;std_psd_w14_1];             % psd - w14-1 near bouy
%Stats = [Stats;mean_psd_w14_2;med_psd_w14_2;std_psd_w14_2];             % psd - w14-2 near bouy

savef = fullfile(path_report,'Stats.txt');
save(savef,'Stats','-ascii');
%% plots
endjd = endjd+1;
if plotit
    % STRATUS-16 intercomparison map plot
    figure; plot(lon1Hr,lat1Hr,'k-'); hold on;
    plot(lon1Hr(aa),lat1Hr(aa),'ro',lon13,lat13,'bx','MarkerSize',8,'LineWidth',1.5); grid;
    xlabel('Lon'); ylabel('Lat'); axis([lon13-1.5*tol,lon13+1.5*tol,lat13-1.5*tol,lat13+1.5*tol]);
    title('STRATUS-16 Intercomparison Points');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_W13_Intercomp',graphformat]);
        print(graphdevice,ppath);
    end

    % STRATUS-18 intercomparison map plot
    figure; plot(lon1Hr,lat1Hr,'k-'); hold on;
    plot(lon1Hr(bb),lat1Hr(bb),'ro',lon14,lat14,'bx','MarkerSize',8,'LineWidth',1.5); grid;
    xlabel('Lon'); ylabel('Lat'); axis([lon14-1.5*tol,lon14+1.5*tol,lat14-1.5*tol,lat14+1.5*tol]);
    title('STRATUS-18 Intercomparison Points');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_SW14_Intercomp',graphformat]);
        print(graphdevice,ppath);
    end

    % wind speed
    figure; plot(jd_ref_hr,spd1hr,'bo',jd_ref_hr,imspd1hr,'bx',jd_ref_hr,w13spd1hr_1,'r+');
    hold on; plot(jd_ref_hr,w13spd1hr_2,'ro',jd_ref_hr,w14spd1hr_1,'g+',jd_ref_hr,w14spd1hr_2,'go');
    grid; xlim([stjd endjd]); xlabel('DOY, 2017'); ylabel('U m/s @ z=3m');
    legend('PSD','SCS','W13-1','W13-2','W14-1','W14-2','location','northwest');
    title('WHOTS-2017 Wind Speed');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_WindSpd',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 2; maxx = 11;
    figure; plot(spd1hr_3m(ff),w13spd1hr_1(ff),'r+',spd1hr_3m(ff),w13spd1hr_2(ff),'ro');
    hold on; plot(spd1hr_3m(gg),w14spd1hr_1(gg),'b+',spd1hr_3m(gg),w14spd1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); grid;
    legend('W13-1','W13-2','W14-1','W14-2','location','northwest');
    ylabel('U m/s, Buoy'); xlabel('U m/s @ 3m, PSD');
    title('WHOTS-2017 Wind Speed, Buoy vs PSD U_{ref} @ 3m');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_WindSpd_v_Uref',graphformat]);
        print(graphdevice,ppath);
    end

    % wind dir
    figure; plot(jd_ref_hr,dir1hr,'bx',jd_ref_hr,imdir1hr,'kd',jd_ref_hr,w13dir1hr_1,'r+');
    hold on; plot(jd_ref_hr,w13dir1hr_2,'ro',jd_ref_hr,w14dir1hr_1,'g+',jd_ref_hr,w14dir1hr_2,'go');
    grid; axis([stjd endjd 0 360]); ylabel('True Wind Direction'); xlabel('DOY, 2017');
    legend('PSD','SCS','W13-1','W13-2','W14-1','W14-2','location','northeast');
    title('WHOTS-2017 Wind Direction Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_WindDir',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 45; maxx = 120;
    figure; plot(dir1hr(ff),w13dir1hr_1(ff),'r+',dir1hr(ff),w13dir1hr_2(ff),'ro');
    hold on; plot(dir1hr(gg),w14dir1hr_1(gg),'b+',dir1hr(gg),w14dir1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); grid; axis([minx,maxx,minx,maxx]);
    legend('W13-1','W13-2','W14-1','W14-2','location','northwest');
    ylabel('Dir True, Buoy'); xlabel('Dir True, PSD');
    title('WHOTS-2017 Wind Direction, Buoy vs PSD');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_WindDir_v_PSD',graphformat]);
        print(graphdevice,ppath);
    end

    % air temperature
    minx = 24; maxx = 29;
    figure; plot(jd_ref_hr,ta1hr,'bx',jd_ref_hr,imta1hr,'bd',jd_ref_hr,w13ta1hr_1,'r+');
    hold on; plot(jd_ref_hr,w13ta1hr_2,'ro',jd_ref_hr,w14ta1hr_1,'g+',jd_ref_hr,w14ta1hr_2,'go');
    grid; axis([stjd,endjd,minx,maxx]); xlabel('DOY, 2017'); ylabel('Air Temp, C @ z=3m');
    legend('PSD','SCS','W13-1','W13-2','W14-1','W14-2','location','southwest');
    ylabel('T air, C'); xlabel('DOY, 2017');
    title('WHOTS-2017 Air Temperature');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_Tair',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 23.5; maxx = 27;
    figure; plot(ta1hr_3m(ff),w13ta1hr_1(ff),'r+',ta1hr_3m(ff),w13ta1hr_2(ff),'ro');
    hold on; plot(ta1hr_3m(gg),w14ta1hr_1(gg),'b+',ta1hr_3m(gg),w14ta1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); grid;
    legend('W13-1','W13-2','W14-1','W14-2','location','northwest');
    ylabel('T air C, Buoy'); xlabel('T air C @ 3m, PSD');
    title('WHOTS-2017 Air Temperature, Buoy vs PSD Tref @ 3m');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_Tair_v_Tref',graphformat]);
        print(graphdevice,ppath);
    end

    % sst
    minx = 26.4; maxx = 27.4;
    figure; plot(jd_ref_hr,ts1hr,'bx',jd_ref_hr,imts1hr,'ko',jd_ref_hr,imtssg1hr,'k+',jd_ref_hr,w13ts1hr_1,'r+');
    hold on; plot(jd_ref_hr,w13ts1hr_2,'ro',jd_ref_hr,w14ts1hr_1,'g+',jd_ref_hr,w14ts1hr_2,'go');
    grid; axis([stjd endjd minx maxx]); xlabel('DOY, 2017'); ylabel('Sea Temp, C');
    legend('tsnk','ts38','ts45','W13-1','W13-2','W14-1','W14-2','location','southeast');
    title('WHOTS-2017 Sea Temperature Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_SST',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 26.4; maxx = 27.4;
    figure; plot(ts1hr(dd),imts1hr(dd),'ko',ts1hr(aa),w13ts1hr_1(aa),'r+',ts1hr(aa),w13ts1hr_2(aa),'ro');
    hold on; plot(ts1hr(bb),w14ts1hr_1(bb),'b+',ts1hr(bb),w14ts1hr_2(bb),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('ts38','W13-1','W13-2','W14-1','W14-2','location','northwest');
    ylabel('SST C, Ship/Buoy'); xlabel('SST C, Sea Snake');
    title('WHOTS-2017 SST, Buoy vs PSD');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_SST_Buoy_v_PSD',graphformat]);
        print(graphdevice,ppath);
    end

    % specific humidity
    figure; plot(jd_ref_hr,qa1hr,'bx');
    hold on; plot(jd_ref_hr,w13qa1hr_1,'r+',jd_ref_hr,w13qa1hr_2,'ro');
    plot(jd_ref_hr,w14qa1hr_1,'g+',jd_ref_hr,w14qa1hr_2,'go');
    grid; xlim([stjd endjd]); xlabel('DOY, 2017'); ylabel('Sp Hum g/kg');
    legend('PSD','W13-1','W13-2','W14-1','W14-2','location','south');
    title('WHOTS-2017 Humidity Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_Qair',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 13; maxx = 18;
    figure; plot(qa1hr_3m(ff),w13qa1hr_1(ff),'r+',qa1hr_3m(ff),w13qa1hr_2(ff),'ro');
    hold on; plot(qa1hr_3m(gg),w14qa1hr_1(gg),'b+',qa1hr_3m(gg),w14qa1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('W13-1','W13-2','W14-1','W14-2','location','northwest');
    ylabel('q g/kg, Buoy'); xlabel('q g/kg @ 3m, PSD');
    title('WHOTS-2017 Specific Humidity, Buoy vs PSD Qref @ 3m');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_Qair_v_Qref',graphformat]);
        print(graphdevice,ppath);
    end

    % RH
    figure; plot(jd_ref_hr,rh1hr,'bx',jd_ref_hr,imrh1hr,'kx');
    hold on; plot(jd_ref_hr,w13rh1hr_1,'r+',jd_ref_hr,w13rh1hr_2,'ro');
    plot(jd_ref_hr,w14rh1hr_1,'go',jd_ref_hr,w14rh1hr_2,'g+');
    grid; xlim([stjd endjd]); xlabel('DOY, 2017'); ylabel('RH %');
    legend('PSD','SCS','W13-1','W13-2','W14-1','W14-2','location','northeast');
    title('WHOTS-2017 RH Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_RH',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 60; maxx = 95;
    figure; plot(rh1hr_3m(ff),w13rh1hr_1(ff),'r+',rh1hr_3m(ff),w13rh1hr_2(ff),'ro');
    hold on; plot(rh1hr_3m(gg),w14rh1hr_1(gg),'b+',rh1hr_3m(gg),w14rh1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('W13-1','W13-2','W14-1','W14-2','location','northwest');
    ylabel('RH, Buoy'); xlabel('RH @ 3m, PSD');
    title('WHOTS-2017 RH, Buoy vs PSD RHref @ 3m');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_RH_v_RHref',graphformat]);
        print(graphdevice,ppath);
    end

    % pressure
    minx = min(Pmb1hr)-1; maxx = max(Pmb1hr)+1;
    figure; plot(jd_ref_hr,Pmb1hr,'bx',jd_ref_hr,imPmb1hr,'kx'); hold on;
    plot(jd_ref_hr,w13Pmb1hr_1,'r+',jd_ref_hr,w13Pmb1hr_2,'ro');
    plot(jd_ref_hr,w14Pmb1hr_1,'g+',jd_ref_hr,w14Pmb1hr_2,'go');
    grid; axis([stjd,endjd,minx,maxx]); xlabel('DOY, 2017'); ylabel('P mb @ sea level');
    legend('PSD','SCS','W13-1','W13-2','W14-1','W14-2','location','north');
    title('WHOTS-2017 Pressure Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_Pressure',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 1014; maxx = 1021;
    figure; plot(Pmb1hr(ff),w13Pmb1hr_1(ff),'r+',Pmb1hr(ff),w13Pmb1hr_2(ff),'ro'); hold on;
    plot(Pmb1hr(gg),w14Pmb1hr_1(gg),'g+',Pmb1hr(gg),w14Pmb1hr_2(gg),'go');
    hold on; plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('W13-1','W13-2','W14-1','W14-2','location','northwest');
    ylabel('P, mb, Buoy'); xlabel('P, mb, PSD');
    title('WHOTS-2017 Sea Level Pressure, Buoy vs PSD');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_Pressure_v_Pref',graphformat]);
        print(graphdevice,ppath);
    end

    % solar radiation
    minx = -20; maxx = 1200;
    figure; plot(jd_ref_hr,rs1hr,'bx');
    hold on; plot(jd_ref_hr,w13rs1hr_1,'r+',jd_ref_hr,w13rs1hr_2,'ro');
    plot(jd_ref_hr,w14rs1hr_1,'g+',jd_ref_hr,w14rs1hr_2,'go');
    grid; axis([stjd endjd minx maxx]); xlabel('DOY, 2017'); ylabel('Solar Downwelling, W/m^2');
    legend('PSD','W13-1','W13-2','W14-1','W14-2','location','northwest');
    title('WHOTS-2017 Solar Radiation Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_PSP',graphformat]);
        print(graphdevice,ppath);
    end

    xx = find(mod(jd_ref_hr,1)<=0.95 & mod(jd_ref_hr,1)>=0.5); % daylight hours only
    yy = intersect(xx,aa); % W13
    zz = intersect(xx,bb); % W14
    minx = min(rs1hr(yy)); maxx = max(rs1hr(yy));
    figure; plot(rs1hr(yy),w13rs1hr_1(yy),'r+',rs1hr(yy),w13rs1hr_2(yy),'ro');
    hold on; plot(rs1hr(zz),w14rs1hr_1(zz),'b+',rs1hr(zz),w14rs1hr_2(zz),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('W13-1','W13-2','W14-1','W14-2','location','northwest');
    ylabel('PSP W/m^2, Buoy'); xlabel('PSP W/m^2, PSD');
    title('WHOTS-2017 Solar Radiation, Buoy vs PSD, daylight hours');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_PSP_Buoy_v_PSD',graphformat]);
        print(graphdevice,ppath);
    end

    % infrared radiation
    minx = 360; maxx = 460;
    figure; plot(jd_ref_hr,rl1hr,'bx');
    hold on; plot(jd_ref_hr,w13rl1hr_1,'r+',jd_ref_hr,w13rl1hr_2,'ro');
    plot(jd_ref_hr,w14rl1hr_1,'g+',jd_ref_hr,w14rl1hr_2,'go');
    grid; axis([stjd endjd minx maxx]); xlabel('DOY, 2017'); ylabel('IR Downwelling, W/m^2');
    legend('PSD','W13-1','W13-2','W14-1','W14-2','location','southwest');
    title('WHOTS-2017 Infrared Radiation Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_PIR',graphformat]);
        print(graphdevice,ppath);
    end

    minx = min(rl1hr)-10; maxx = max(rl1hr)+10;
    figure; plot(rl1hr(aa),w13rl1hr_1(aa),'r+',rl1hr(aa),w13rl1hr_2(aa),'ro');
    hold on; plot(rl1hr(bb),w14rl1hr_1(bb),'b+',rl1hr(bb),w14rl1hr_2(bb),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('W13-1','W13-2','W14-1','W14-2','location','northwest');
    ylabel('PIR W/m^2, Buoy'); xlabel('PIR W/m^2, PSD');
    title('WHOTS-2017 IR Radiation, Buoy vs PSD');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_PIR_Buoy_v_PSD',graphformat]);
        print(graphdevice,ppath);
    end

    % daytime PIR bias?
    tod = mod(jd_ref_hr,1)*24;
    [tod_bin,W13_1_diff_bin,~,~,~,~,~] = binave(tod(aa),w13rl1hr_1(aa)-rl1hr(aa),2,0,24);
    [~,W13_2_diff_bin,~,~,~,~,~] = binave(tod(aa),w13rl1hr_2(aa)-rl1hr(aa),2,0,24);
%     [~,W14_1_diff_bin,~,~,~,~,~] = binave(tod(bb),w14rl1hr_1(bb)-rl1hr(bb),2,0,24);
%     [~,W14_2_diff_bin,~,~,~,~,~] = binave(tod(bb),w14rl1hr_2(bb)-rl1hr(bb),2,0,24);
    tod_bin = tod_bin+1;

    figure; plot(tod_bin,W13_1_diff_bin,tod_bin,W13_2_diff_bin); hold on;
    plot(tod_bin,W14_1_diff_bin,tod_bin,W14_2_diff_bin,[0,24],[0,0],'k-'); grid;
    xlim([0,24]); xlabel('Hour of Day'); ylabel('Rl difference to mean PSD Rl');
    title('Daytime Rl Bias'); legend('W13-1','W13-2','W14-1','W14-2','location','south');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_PIR_Bias',graphformat]);
        print(graphdevice,ppath);
    end

    minx = -20; maxx = 1200;
    figure;
    plot(rs1hr(yy),w13rl1hr_1(yy)-rl1hr(yy),'bx',rs1hr(yy),w13rl1hr_2(yy)-rl1hr(yy),'gx'); hold on;
    plot(rs1hr(zz),w14rl1hr_1(zz)-rl1hr(zz),'rx',rs1hr(zz),w14rl1hr_2(zz)-rl1hr(zz),'kx');
    plot([minx,maxx],[0,0],'k-'); xlim([minx,maxx]); grid;
    xlabel('solar flux, W/m^2, daytime only'); ylabel('Rl_{scs} - Rl_{psd}, W/m^2');
    legend('W13-1','W13-2','W14-1','W14-2','location','best'); title('Daytime Rl Bias');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_PIR_Bias_v_solar',graphformat]);
        print(graphdevice,ppath);
    end
end

