% STRATUS_2017_buoy_comparison.m

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

%% initialize
clear;
close all;

% matlab script path
path_prog = fullfile('/Users/byronb/Documents/MATLAB/2017_STRATUS/STRATUS_comparison');
restoredefaultpath              % reset matlab path
addpath(genpath(path_prog));	% add scripts to path
savepath; rehash;

%% define data paths & constants
path_data = '/Users/byronb/Raw_Data_Files.noindex/STRATUS_2017/RHB/flux/Processed_STRATUS_only';
path_whoi = '/Users/byronb/Raw_Data_Files.noindex/STRATUS_2017/WHOI';
path_plots = '/Users/byronb/Raw_Data_Files.noindex/STRATUS_2017/RHB/flux/buoy_comparison';
mkdir(path_plots);
path_report = '/Users/byronb/Raw_Data_Files.noindex/STRATUS_2017/RHB/Scientific_analysis/report';
mkdir(path_report);

% start and end dates for the buoy intercomparison
%  118:140 for whole STRATUS leg
%  132:137.5 in vicinity of buoy
stjd = 132;
endjd = 138;

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

% sensor heights
% as measured on CalWater 01/13/2015: waterline to base of tower = 7.16 m.
zu = 18.16;     % wind speed measurement height (m): 7.16m + 11m
zt = 16.5;      % air T measurement height (m)
zq = 16.5;      % air q measurement height (m)
zP = 17;        % psd pressure sensor height (m)
zWXT = 17;
zuim = 15;      % ship SCS sensor heights
ztim = 15;
zqim = 15;
zpim = 15;

% depth of ship tsgm, tssm
dtsg = 5;

%% load 10 min and 1 hr met/flux data arrays
%{
    flux10 is 3312x64
    COL   DESCRIPTION
    1     decimal day of year (JD)
    2     uship
    3     u true
    4     true dir
    5     u rel
    6     rel wind dir, 0-360 deg
    7     ship heading
    8     Ts sea snake
    9     T
    10    qs
    11    qa
    12    hs covariance
    13    hs ID
    14    hs bulk
    15    hl covariance from licor
    16    hl ID
    17    hl bulk
    18    stress streamwise covariance
    19    stress  xstream covariance
    20    stress ID
    21    stress bulk
    22    solar flux
    23    IR flux
    24    rain rate
    25    ship plume contam index
    26    tilt angle
    27    maneuver index
    28    ct
    29    cq from Licor
    30    cu
    31    cw
    32    Rain heat flux
    33    Hl webb flux
    34    lat
    35    lon
    36    zu
    37    zt
    38    zq
    39    uship scs
    40    u true scs
    41    true dir scs
    42    ship course scs
    43    Ts scs
    44    T scs
    45    qs scs
    46    qa scs
    47    solar flux  IMET
    48    IR flux IMET
    49    <w'c'> Licor - NOT COMPUTED
    50    Specific Humidity from Licor (g/kg)
    51    Stv of Specific Humidity from Licor (g/kg)
    52    CO2 concentration from Licor (umol/mol)
    53    Stv of CO2 concentration from Licor (umol/mol)
    54    atmospheric pressure (mbar), sea level
    55    Wind speed (m/s) relative to earth
    56    Wind direction (deg) from relative to earth
    57    Longwave downwelling clear sky model (W/m2)
    58    Shortwave solar radiance clear sky model (W/m2)
    59    10m neutral windspeed, COARE 3.5
    60    10m neutral temperature, COARE 3.5
    61    10m neutral specific humidity, COARE 3.5
    62    3m neutral windspeed, COARE 3.5
    63    3m neutral temperature, COARE 3.5
    64    3m RH, COARE 3.5
    65    SCS RH, %
    66    SCS Pmb, sea level
    67    PSD WXT RH, %
    
    flux1hr is 552x62
    COL   DESCRIPTION
    1     decimal day of year at start of hour (JD)
    2     uship
    3     u true
    4     true dir
    5     u rel
    6     rel wind dir, 0-360 deg
    8     SST sea snake
    9     Ta best
    10    qs best
    11    qa best
    12    hs covariance
    13    hs ID
    14    hs bulk
    15    hl blend of cov results
    16    hl ID
    17    hl bulk
    18    stress stream-wise covariance
    19    stress cross-stream covariance
    20    stress ID
    21    stress bulk
    22    rs best
    23    rl best
    24    ORG rain rate
    25    ship plume contam index
    26    tilt of good data
    27    maneuver index
    28    ct
    29    cql
    30    cu
    31    cw
    32    Rain heat flux
    33    Hl webb flux
    34    lat
    35    lon
    36    zu
    37    zt
    38    zq
    39    uship scs
    40    u true scs
    41    true dir scs
    42    ship course scs
    43    Ts scs, raw ship measurement
    44    T scs
    45    qs scs
    46    qa scs
    47    rs scs
    48    rl scs
    49    w'c' licor - NOT COMPUTED
    50    Specific Humidity from Licor (g/kg)
    51    Stv of Specific Humidity from Licor (g/kg)
    52    CO2 concentration from Licor (umol/mol)
    53    Stv of CO2 concentration from Licor (umol/mol)
    54    atmospheric pressure best (mbar)
    55    Wind speed (m/s) relative to earth
    56    Wind direction (deg) from relative to earth
    57    net longwave heat flux
    58    net heat flux, bulk, W/m2 hourly
    59    net heat flux, turb, W/m2 hourly
    60    3m neutral windspeed, COARE 3.5
    61    3m neutral temperature, COARE 3.5
    62    3m RH, COARE 3.5
    63    SCS RH, %
    64    SCS Pmb, sea level
    65    PSD WXT RH, %
%}

load(fullfile(path_data,'STRATUS_2017_flux_10_v4_decorr.mat'));
ii = flux10(:,1)>=stjd & flux10(:,1)<endjd;
flux10 = flux10(ii,:);

flux1hr = load(fullfile(path_data,'STRATUS_2017_flux_hr_v4_decorr.txt'));
ii = flux1hr(:,1)>=stjd & flux1hr(:,1)<endjd;
flux1hr = flux1hr(ii,:);

jd_hr = flux1hr(:,1);

%% load buoy data

% load STRATUS-15 1 min data: converted to N x 14 array
% yday BPR ATMP SST WNDE WNDN WSPD WDIR RH SWR LWR PRATE SAL COND
%  1    2   3    4   5    6    7    8   9   10  11   12  13   14
load(fullfile(path_whoi,'s15met_may2017.mat'));  % structures sys1 and sys2
S15_1 = NaN(19770,14); S15_2 = S15_1;
S15_1(:,1) = sys1.mday - datenum(2016,12,31); % jd
S15_1(:,2) = sys1.bpr + 0.125*3; % sea level
S15_1(:,3) = sys1.atmp;
S15_1(:,4) = sys1.sst;
S15_1(:,5) = sys1.wnde;
S15_1(:,6) = sys1.wndn;
S15_1(:,7) = sys1.wspd;
[~,temp] = uv_to_sd(sys1.wndn,-sys1.wnde);
S15_1(:,8) = temp;  % wdir
S15_1(:,9) = sys1.hrh;
S15_1(:,10) = sys1.swr;
S15_1(:,11) = sys1.lwr;
S15_1(:,12) = sys1.precip;
S15_1(:,13) = sys1.sal;
S15_1(:,14) = sys1.cond;
zz = S15_1(:,1)>135.8;
S15_1(zz,2:end) = NaN;

S15_2(:,1) = sys2.mday - datenum(2016,12,31);
S15_2(:,2) = sys2.bpr + 0.125*3; % sea level
S15_2(:,3) = sys2.atmp;
S15_2(:,4) = sys2.sst;
S15_2(:,5) = sys2.wnde;
S15_2(:,6) = sys2.wndn;
S15_2(:,7) = sys2.wspd;
[~,temp] = uv_to_sd(sys2.wndn,-sys2.wnde);
S15_2(:,8) = temp;  % wdir
S15_2(:,9) = sys2.hrh;
S15_2(:,10) = sys2.swr;
S15_2(:,11) = sys2.lwr;
S15_2(:,12) = sys2.precip;
S15_2(:,13) = sys2.sal;
S15_2(:,14) = sys2.cond;
zz = S15_2(:,1)>135.8;
S15_2(zz,2:end) = NaN;

% average S15 to 10-min & 1 hr timebase
jd_10_bin = stjd:10/1440:endjd; jd_10_bin = jd_10_bin';
S15_1_10 = interval_avg(S15_1(:,1),S15_1(:,2:end),jd_10_bin);
[~,temp] = uv_to_sd(S15_1_10(:,6),-S15_1_10(:,5));
S15_1_10(:,8) = temp;  % wdir
S15_2_10 = interval_avg(S15_2(:,1),S15_2(:,2:end),jd_10_bin);
[~,temp] = uv_to_sd(S15_2_10(:,6),-S15_2_10(:,5));
S15_2_10(:,8) = temp;  % wdir
jd_hr_bin = stjd:1/24:endjd; jd_hr_bin = jd_hr_bin';
S15_1_hr = interval_avg(S15_1(:,1),S15_1(:,2:end),jd_hr_bin);
[~,temp] = uv_to_sd(S15_1_hr(:,6),-S15_1_hr(:,5));
S15_1_hr(:,8) = temp;  % wdir
S15_2_hr = interval_avg(S15_2(:,1),S15_2(:,2:end),jd_hr_bin);
[~,temp] = uv_to_sd(S15_2_hr(:,6),-S15_2_hr(:,5));
S15_2_hr(:,8) = temp;  % wdir

% 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
%    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/hr ppt  S/m

S16_1_hr = load(fullfile(path_whoi,'Stratus16_sys1_nohdr.txt'));
S16_1_hr(:,5) = datenum(S16_1_hr(:,1),S16_1_hr(:,2),S16_1_hr(:,3),S16_1_hr(:,4),0,0) - datenum(2016,12,31);
S16_1_hr = [S16_1_hr(:,5:15),S16_1_hr(:,17:19)];
ii = S16_1_hr(:,1)>=stjd & S16_1_hr(:,1)<endjd;
S16_1_hr = S16_1_hr(ii,:);
S16_1_hr = [NaN(46,14);S16_1_hr];
S16_1_hr(:,1) = S15_1_hr(:,1);
S16_1_hr(:,8) = S16_1_hr(:,8)+180;      % convert wdir to 'from' convention
S16_1_hr(:,2) = S16_1_hr(:,2)+0.125*3;  % sea level pressure
S16_1_hr(:,8) = mod(S16_1_hr(:,8),360);

S16_2_hr = load(fullfile(path_whoi,'Stratus16_sys2_nohdr.txt'));
S16_2_hr(:,5) = datenum(S16_2_hr(:,1),S16_2_hr(:,2),S16_2_hr(:,3),S16_2_hr(:,4),0,0) - datenum(2016,12,31);
S16_2_hr = [S16_2_hr(:,5:15),S16_2_hr(:,17:19)];
ii = S16_2_hr(:,1)>=stjd & S16_2_hr(:,1)<endjd;
S16_2_hr = S16_2_hr(ii,:);
S16_2_hr = [NaN(46,14);S16_2_hr];
S16_2_hr(:,1) = S15_2_hr(:,1);
S16_2_hr(:,8) = S16_2_hr(:,8)+180;      % convert wdir to 'from' convention
S16_2_hr(:,2) = S16_2_hr(:,2)+0.125*3;  % sea level pressure
S16_2_hr(:,8) = mod(S16_2_hr(:,8),360);

%% define variables
% PSD
lat10 = flux10(:,34); lat1Hr = flux1hr(:,34);
lon10 = flux10(:,35); lon1Hr = flux1hr(:,35);
rwdir10 = flux10(:,6); rwdir10(rwdir10>180) = rwdir10(rwdir10>180)-360;
rwdir1hr = flux1hr(:,6); rwdir1hr(rwdir1hr>180) = rwdir1hr(rwdir1hr>180)-360;
spd10 = flux10(:,3); spd10_3m = flux10(:,62); dir10 = flux10(:,4);
spd1hr = flux1hr(:,3); spd1hr_3m = flux1hr(:,60); dir1hr = flux1hr(:,4);
ta10 = flux10(:,9); ta1hr = flux1hr(:,9);
ta10_3m = flux10(:,63); ta1hr_3m = flux1hr(:,61);
ts10 = flux10(:,8); ts1hr = flux1hr(:,8);
rh10 = flux10(:,67); rh1hr = flux1hr(:,65);
rh10_3m = flux10(:,64); rh1hr_3m = flux1hr(:,62);
rs10 = flux10(:,22); rs1hr = flux1hr(:,22);
rl10 = flux10(:,23); rl1hr = flux1hr(:,23);
Pmb10 = flux10(:,54); Pmb1hr = flux1hr(:,54);
qa10 = flux10(:,11); qa1hr = flux1hr(:,11);
qa10_3m = qair_p([ta10_3m,rh10_3m,Pmb10]);
qa1hr_3m = qair_p([ta1hr_3m,rh1hr_3m,Pmb1hr]);

% ship
imspd10 = flux10(:,40); imdir10 = flux10(:,41);
imspd1hr = flux1hr(:,40); imdir1hr = flux1hr(:,41);
imta10 = flux10(:,44); imta1hr = flux1hr(:,44);
imts10 = flux10(:,43); imts1hr = flux1hr(:,43);
imrh10 = flux10(:,65); imrh1hr = flux1hr(:,63);
imrs10 = flux10(:,47); imrs1hr = flux1hr(:,47);
imrl10 = flux10(:,48); imrl1hr = flux1hr(:,48);
imPmb10 = flux10(:,66); imPmb1hr = flux1hr(:,64);
imqa10 = flux10(:,46); imqa1hr = flux1hr(:,46);


%S15
s15spd10_1 = S15_1_10(:,7); s15spd1hr_1 = S15_1_hr(:,7);
s15dir10_1 = S15_1_10(:,8); s15dir1hr_1 = S15_1_hr(:,8);
s15ta10_1 = S15_1_10(:,3); s15ta1hr_1 = S15_1_hr(:,3);
s15ts10_1 = S15_1_10(:,4); s15ts1hr_1 = S15_1_hr(:,4);
s15rh10_1 = S15_1_10(:,9); s15rh1hr_1 = S15_1_hr(:,9);
s15rs10_1 = S15_1_10(:,10); s15rs1hr_1 = S15_1_hr(:,10);
s15rl10_1 = S15_1_10(:,11); s15rl1hr_1 = S15_1_hr(:,11);
s15Pmb10_1 = S15_1_10(:,2); s15Pmb1hr_1 = S15_1_hr(:,2);
s15qa10_1 = qair_p([s15ta10_1,s15rh10_1,s15Pmb10_1]);
s15qa1hr_1 = qair_p([s15ta1hr_1,s15rh1hr_1,s15Pmb1hr_1]);

s15spd10_2 = S15_2_10(:,7); s15spd1hr_2 = S15_2_hr(:,7);
s15dir10_2 = S15_2_10(:,8); s15dir1hr_2 = S15_2_hr(:,8);
s15ta10_2 = S15_2_10(:,3); s15ta1hr_2 = S15_2_hr(:,3);
s15ts10_2 = S15_2_10(:,4); s15ts1hr_2 = S15_2_hr(:,4);
s15rh10_2 = S15_2_10(:,9); s15rh1hr_2 = S15_2_hr(:,9);
s15rs10_2 = S15_2_10(:,10); s15rs1hr_2 = S15_2_hr(:,10);
s15rl10_2 = S15_2_10(:,11); s15rl1hr_2 = S15_2_hr(:,11);
s15Pmb10_2 = S15_2_10(:,2); s15Pmb1hr_2 = S15_2_hr(:,2);
s15qa10_2 = qair_p([s15ta10_2,s15rh10_2,s15Pmb10_2]);
s15qa1hr_2 = qair_p([s15ta1hr_2,s15rh1hr_2,s15Pmb1hr_2]);

%S16 - hourly avg only
s16spd1hr_1 = S16_1_hr(:,7);
s16dir1hr_1 = S16_1_hr(:,8);
s16ta1hr_1 = S16_1_hr(:,3);
s16ts1hr_1 = S16_1_hr(:,4);
s16rh1hr_1 = S16_1_hr(:,9);
s16rs1hr_1 = S16_1_hr(:,10);
s16rl1hr_1 = S16_1_hr(:,11);
s16Pmb1hr_1 = S16_1_hr(:,2);
s16qa1hr_1 = qair_p([s16ta1hr_1,s16rh1hr_1,s16Pmb1hr_1]);


s16spd1hr_2 = S16_2_hr(:,7);
s16dir1hr_2 = S16_2_hr(:,8);
s16ta1hr_2 = S16_2_hr(:,3);
s16ts1hr_2 = S16_2_hr(:,4);
s16rh1hr_2 = S16_2_hr(:,9);
s16rs1hr_2 = S16_2_hr(:,10);
s16rl1hr_2 = S16_2_hr(:,11);
s16Pmb1hr_2 = S16_2_hr(:,2);
s16qa1hr_2 = qair_p([s16ta1hr_2,s16rh1hr_2,s16Pmb1hr_2]);

%% define selectors

% define mooring locations
lat15 = -19.6262; lon15 = 275.0530;
lat16 = -19.4302; lon16 = 274.9262;

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

% date selector for STRATUS-15 intercomparison
%  134.1 to 135.9 on station at STRATUS-15
bb = find(jd_hr>=134.1 & jd_hr<=135.9);

% date selector for STRATUS-16 intercomparison
%  136.1 to 136.8 on station at STRATUS-16
cc = find((jd_hr>=136.1 & jd_hr<=136.9));

% selector for entire period at site
dd = find(jd_hr>=132.2 & jd_hr<=136.8);

ff = intersect(aa,bb); % S15 wind sector and proximity selector
gg = intersect(aa,cc); % S16 wind sector and proximity selector

hh = intersect(aa,dd); % general wind sector for entire period

%% build hourly data arrays of [Ta,RH,SST,U,dir,Rs,Rl,P] for each platform
none = NaN(size(jd_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),imrs1hr,imrl1hr,imPmb1hr];

s15_1 = [s15ta1hr_1,s15rh1hr_1,s15ts1hr_1,s15spd1hr_1,unwrap(s15dir1hr_1*pi/180),...
       s15rs1hr_1,s15rl1hr_1,s15Pmb1hr_1];
s15_2 = [s15ta1hr_2,s15rh1hr_2,s15ts1hr_2,s15spd1hr_2,unwrap(s15dir1hr_2*pi/180),...
       s15rs1hr_2,s15rl1hr_2,s15Pmb1hr_2];
s16_1 = [s16ta1hr_1,s16rh1hr_1,s16ts1hr_1,s16spd1hr_1,unwrap(s16dir1hr_1*pi/180),...
       s16rs1hr_1,s16rl1hr_1,s16Pmb1hr_1];
s16_2 = [s16ta1hr_2,s16rh1hr_2,s16ts1hr_2,s16spd1hr_2,unwrap(s16dir1hr_2*pi/180),...
       s16rs1hr_2,s16rl1hr_2,s16Pmb1hr_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
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;
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-15 stats for entire period at site
mean_psd_s15_1_all = nanmean1(psd(dd,:) - s15_1(dd,:));
mean_psd_s15_1_all(jj) = nanmean1(psd(hh,jj) - s15_1(hh,jj));
mean_psd_s15_1_all(5) = nanmean1(psd(hh,5) - s15_1(hh,5))*180/pi;
med_psd_s15_1_all = nanmedian1(psd(dd,:) - s15_1(dd,:));
med_psd_s15_1_all(jj) = nanmedian1(psd(hh,jj) - s15_1(hh,jj));
med_psd_s15_1_all(5) = nanmedian1(psd(hh,5) - s15_1(hh,5))*180/pi;
std_psd_s15_1_all = nanstd1(psd(dd,:) - s15_1(dd,:));
std_psd_s15_1_all(jj) = nanstd1(psd(hh,jj) - s15_1(hh,jj));
std_psd_s15_1_all(5) = nanstd1(psd(hh,5) - s15_1(hh,5))*180/pi;

mean_psd_s15_2_all = nanmean1(psd(dd,:) - s15_2(dd,:));
mean_psd_s15_2_all(jj) = nanmean1(psd(hh,jj) - s15_2(hh,jj));
mean_psd_s15_2_all(5) = nanmean1(psd(hh,5) - s15_2(hh,5))*180/pi;
med_psd_s15_2_all = nanmedian1(psd(dd,:) - s15_2(dd,:));
med_psd_s15_2_all(jj) = nanmedian1(psd(hh,jj) - s15_2(hh,jj));
med_psd_s15_2_all(5) = nanmedian1(psd(hh,5) - s15_2(hh,5))*180/pi;
std_psd_s15_2_all = nanstd1(psd(dd,:) - s15_2(dd,:));
std_psd_s15_2_all(jj) = nanstd1(psd(hh,jj) - s15_2(hh,jj));
std_psd_s15_2_all(5) = nanstd1(psd(hh,5) - s15_2(hh,5))*180/pi;

% psd vs STRATUS-15 stats for period near bouy
mean_psd_s15_1 = nanmean1(psd(bb,:) - s15_1(bb,:));
mean_psd_s15_1(jj) = nanmean1(psd(ff,jj) - s15_1(ff,jj));
mean_psd_s15_1(5) = nanmean1(psd(ff,5) - s15_1(ff,5))*180/pi;
med_psd_s15_1 = nanmedian1(psd(bb,:) - s15_1(bb,:));
med_psd_s15_1(jj) = nanmedian1(psd(ff,jj) - s15_1(ff,jj));
med_psd_s15_1(5) = nanmedian1(psd(ff,5) - s15_1(ff,5))*180/pi;
std_psd_s15_1 = nanstd1(psd(bb,:) - s15_1(bb,:));
std_psd_s15_1(jj) = nanstd1(psd(ff,jj) - s15_1(ff,jj));
std_psd_s15_1(5) = nanstd1(psd(ff,5) - s15_1(ff,5))*180/pi;

mean_psd_s15_2 = nanmean1(psd(bb,:) - s15_2(bb,:));
mean_psd_s15_2(jj) = nanmean1(psd(ff,jj) - s15_2(ff,jj));
mean_psd_s15_2(5) = nanmean1(psd(ff,5) - s15_2(ff,5))*180/pi;
med_psd_s15_2 = nanmedian1(psd(bb,:) - s15_2(bb,:));
med_psd_s15_2(jj) = nanmedian1(psd(ff,jj) - s15_2(ff,jj));
med_psd_s15_2(5) = nanmedian1(psd(ff,5) - s15_2(ff,5))*180/pi;
std_psd_s15_2 = nanstd1(psd(bb,:) - s15_2(bb,:));
std_psd_s15_2(jj) = nanstd1(psd(ff,jj) - s15_2(ff,jj));
std_psd_s15_2(5) = nanstd1(psd(ff,5) - s15_2(ff,5))*180/pi;

% psd vs STRATUS-16 stats for entire period at site
mean_psd_s16_1_all = nanmean1(psd(dd,:) - s16_1(dd,:));
mean_psd_s16_1_all(jj) = nanmean1(psd(hh,jj) - s16_1(hh,jj));
mean_psd_s16_1_all(5) = nanmean1(psd(hh,5) - s16_1(hh,5))*180/pi;
med_psd_s16_1_all = nanmedian1(psd(dd,:) - s16_1(dd,:));
med_psd_s16_1_all(jj) = nanmedian1(psd(hh,jj) - s16_1(hh,jj));
med_psd_s16_1_all(5) = nanmedian1(psd(hh,5) - s16_1(hh,5))*180/pi;
std_psd_s16_1_all = nanstd1(psd(dd,:) - s16_1(dd,:));
std_psd_s16_1_all(jj) = nanstd1(psd(hh,jj) - s16_1(hh,jj));
std_psd_s16_1_all(5) = nanstd1(psd(hh,5) - s16_1(hh,5))*180/pi;

mean_psd_s16_2_all = nanmean1(psd(dd,:) - s16_2(dd,:));
mean_psd_s16_2_all(jj) = nanmean1(psd(hh,jj) - s16_2(hh,jj));
mean_psd_s16_2_all(5) = nanmean1(psd(hh,5) - s16_2(hh,5))*180/pi;
med_psd_s16_2_all = nanmedian1(psd(dd,:) - s16_2(dd,:));
med_psd_s16_2_all(jj) = nanmedian1(psd(hh,jj) - s16_2(hh,jj));
med_psd_s16_2_all(5) = nanmedian1(psd(hh,5) - s16_2(hh,5))*180/pi;
std_psd_s16_2_all = nanstd1(psd(dd,:) - s16_2(dd,:));
std_psd_s16_2_all(jj) = nanstd1(psd(hh,jj) - s16_2(hh,jj));
std_psd_s16_2_all(5) = nanstd1(psd(hh,5) - s16_2(hh,5))*180/pi;

% psd vs STRATUS-16 stats for period near bouy
mean_psd_s16_1 = nanmean1(psd(cc,:) - s16_1(cc,:));
mean_psd_s16_1(jj) = nanmean1(psd(gg,jj) - s16_1(gg,jj));
mean_psd_s16_1(5) = nanmean1(psd(gg,5) - s16_1(gg,5))*180/pi;
med_psd_s16_1 = nanmedian1(psd(cc,:) - s16_1(cc,:));
med_psd_s16_1(jj) = nanmedian1(psd(gg,jj) - s16_1(gg,jj));
med_psd_s16_1(5) = nanmedian1(psd(gg,5) - s16_1(gg,5))*180/pi;
std_psd_s16_1 = nanstd1(psd(cc,:) - s16_1(cc,:));
std_psd_s16_1(jj) = nanstd1(psd(gg,jj) - s16_1(gg,jj));
std_psd_s16_1(5) = nanstd1(psd(gg,5) - s16_1(gg,5))*180/pi;

mean_psd_s16_2 = nanmean1(psd(cc,:) - s16_2(cc,:));
mean_psd_s16_2(jj) = nanmean1(psd(gg,jj) - s16_2(gg,jj));
mean_psd_s16_2(5) = nanmean1(psd(gg,5) - s16_2(gg,5))*180/pi;
med_psd_s16_2 = nanmedian1(psd(cc,:) - s16_2(cc,:));
med_psd_s16_2(jj) = nanmedian1(psd(gg,jj) - s16_2(gg,jj));
med_psd_s16_2(5) = nanmedian1(psd(gg,5) - s16_2(gg,5))*180/pi;
std_psd_s16_2 = nanstd1(psd(cc,:) - s16_2(cc,:));
std_psd_s16_2(jj) = nanstd1(psd(gg,jj) - s16_2(gg,jj));
std_psd_s16_2(5) = nanstd1(psd(gg,5) - s16_2(gg,5))*180/pi;

%% save stats to table format for import to excel
% 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_s15_1_all;med_psd_s15_1_all;std_psd_s15_1_all]; % psd - S15-1 for all
Stats = [Stats;mean_psd_s15_2_all;med_psd_s15_2_all;std_psd_s15_2_all]; % psd - S15-2 for all
Stats = [Stats;mean_psd_s16_1_all;med_psd_s16_1_all;std_psd_s16_1_all]; % psd - S16-1 for all
Stats = [Stats;mean_psd_s16_2_all;med_psd_s16_2_all;std_psd_s16_2_all]; % psd - S16-2 for all
Stats = [Stats;mean_psd_s15_1;med_psd_s15_1;std_psd_s15_1];             % psd - S15-1 near bouy
Stats = [Stats;mean_psd_s15_2;med_psd_s15_2;std_psd_s15_2];             % psd - S15-2 near bouy
Stats = [Stats;mean_psd_s16_1;med_psd_s16_1;std_psd_s16_1];             % psd - S15-1 near bouy
Stats = [Stats;mean_psd_s16_2;med_psd_s16_2;std_psd_s16_2];             % psd - S16-2 near bouy

savef = fullfile(path_report,'Stats.txt');
save(savef,'Stats','-ascii');

%% plots

if plotit
    % STRATUS-15 intercomparison map plot
    figure; plot(lon10,lat10,'k-'); hold on;
    plot(lon1Hr(ff),lat1Hr(ff),'ro',lon15,lat15,'kx','MarkerSize',8,'LineWidth',1.5); grid;
    xlabel('Lon'); ylabel('Lat'); axis([275.04,275.16,-19.76,-19.66]);
    title('STRATUS-15 Intercomparison Points');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_S15_Intercomp',graphformat]);
        print(graphdevice,ppath);
    end
    
    % STRATUS-16 intercomparison map plot
    figure; plot(lon10,lat10,'k-'); hold on;
    plot(lon1Hr(cc),lat1Hr(cc),'ro',lon16,lat16,'kx','MarkerSize',8,'LineWidth',1.5); grid;
    xlabel('Lon'); ylabel('Lat'); axis([274.86,274.96,-19.46,-19.38]);
    title('STRATUS-16 Intercomparison Points');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_S16_Intercomp',graphformat]);
        print(graphdevice,ppath);
    end
    
    % wind speed
    figure; plot(jd_hr,spd1hr,'bo',jd_hr,imspd1hr,'bx',jd_hr,s15spd1hr_1,'r+');
    hold on; plot(jd_hr,s15spd1hr_2,'ro',jd_hr,s16spd1hr_1,'g+',jd_hr,s16spd1hr_2,'go');
    grid; xlim([stjd endjd]); xlabel('DOY, 2017'); ylabel('U m/s @ z=3m');
    legend('PSD','SCS','S15-1','S15-2','S16-1','S16-2','location','northwest');
    title('STRATUS-2017 Wind Speed');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_WindSpd',graphformat]);
        print(graphdevice,ppath);
    end
    
    minx = 3; maxx = 9;
    figure; plot(spd1hr_3m(ff),s15spd1hr_1(ff),'r+',spd1hr_3m(ff),s15spd1hr_2(ff),'ro');
    hold on; plot(spd1hr_3m(gg),s16spd1hr_1(gg),'b+',spd1hr_3m(gg),s16spd1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); grid;
    legend('S15-1','S15-2','S16-1','S16-2','location','northwest');
    ylabel('U m/s, Buoy'); xlabel('U m/s @ 3m, PSD');
    title('STRATUS-2017 Wind Speed, Buoy vs PSD U_{ref} @ 3m');
    if printit
        ppath = fullfile(path_plots,['STRATUS-2017_WindSpd_v_Uref',graphformat]);
        print(graphdevice,ppath);
    end

    % wind dir
    figure; plot(jd_hr,dir1hr,'bx',jd_hr,imdir1hr,'kd',jd_hr,s15dir1hr_1,'r+');
    hold on; plot(jd_hr,s15dir1hr_2,'ro',jd_hr,s16dir1hr_1,'g+',jd_hr,s16dir1hr_2,'go');
    grid; axis([stjd endjd 0 360]); ylabel('True Wind Direction'); xlabel('DOY, 2017');
    legend('PSD','SCS','S15-1','S15-2','S16-1','S16-2','location','northeast');
    title('STRATUS-2017 Wind Direction Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_WindDir',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 80; maxx = 160;
    figure; plot(dir1hr(ff),s15dir1hr_1(ff),'r+',dir1hr(ff),s15dir1hr_2(ff),'ro');
    hold on; plot(dir1hr(gg),s16dir1hr_1(gg),'b+',dir1hr(gg),s16dir1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); grid; axis([minx,maxx,minx,maxx]);
    legend('S15-1','S15-2','S16-1','S16-2','location','northwest');
    ylabel('Dir True, Buoy'); xlabel('Dir True, PSD');
    title('STRATUS-2017 Wind Direction, Buoy vs PSD');
    if printit
        ppath = fullfile(path_plots,['STRATUS-2017_WindDir_v_PSD',graphformat]);
        print(graphdevice,ppath);
    end

    % air temperature
    figure; plot(jd_hr,ta1hr,'bx',jd_hr,imta1hr,'bd',jd_hr,s15ta1hr_1,'r+');
    hold on; plot(jd_hr,s15ta1hr_2,'ro',jd_hr,s16ta1hr_1,'g+',jd_hr,s16ta1hr_2,'go');
    grid; axis([stjd,endjd,19,23]); xlabel('DOY, 2017'); ylabel('Air Temp, C @ z=3m');
    legend('PSD','SCS','S15-1','S15-2','S16-1','S16-2','location','south');
    ylabel('T air, C'); xlabel('DOY, 2017');
    title('STRATUS-2017 Air Temperature');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_Tair',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 19; maxx = 23;
    figure; plot(ta1hr_3m(ff),s15ta1hr_1(ff),'r+',ta1hr_3m(ff),s15ta1hr_2(ff),'ro');
    hold on; plot(ta1hr_3m(gg),s16ta1hr_1(gg),'b+',ta1hr_3m(gg),s16ta1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); grid;
    legend('S15-1','S15-2','S16-1','S16-2','location','northwest');
    ylabel('T air C, Buoy'); xlabel('T air C @ 3m, PSD');
    title('STRATUS-2017 Air Temperature, Buoy vs PSD Tref @ 3m');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_Tair_v_Tref',graphformat]);
        print(graphdevice,ppath);
    end

    % sst
    figure; plot(jd_hr,ts1hr,'bx',jd_hr,imts1hr,'bd',jd_hr,s15ts1hr_1,'r+');
    hold on; plot(jd_hr,s15ts1hr_2,'ro',jd_hr,s16ts1hr_1,'g+',jd_hr,s16ts1hr_2,'go');
    grid; axis([stjd endjd 22 24]); xlabel('DOY, 2017'); ylabel('Sea Temp, C');
    legend('tsnk','tsg','S15-1','S15-2','S16-1','S16-2','location','southwest');
    title('STRATUS-2017 Sea Temperature Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_SST',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 22.6; maxx = 23.2;
    figure; plot(ts1hr(bb),s15ts1hr_1(bb),'r+',ts1hr(bb),s15ts1hr_2(bb),'ro');
    hold on; plot(ts1hr(cc),s16ts1hr_1(cc),'b+',ts1hr(cc),s16ts1hr_2(cc),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('S15-1','S15-2','S16-1','S16-2','location','northwest');
    ylabel('SST C, Buoy'); xlabel('SST C, Sea Snake');
    title('STRATUS-2017 SST, Buoy vs PSD');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_SST_Buoy_v_PSD',graphformat]);
        print(graphdevice,ppath);
    end
    
%     minx = 22.6; maxx = 23.75;
%     figure; plot(imts1hr(bb),s15ts1hr_1(bb),'r+',imts1hr(bb),s15ts1hr_2(bb),'ro');
%     hold on; plot(imts1hr(cc),s16ts1hr_1(cc),'b+',imts1hr(cc),s16ts1hr_2(cc),'bo');
%     plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
%     legend('S15-1','S15-2','S16-1','S16-2','location','northwest');
%     ylabel('SST C, Buoy'); xlabel('SST C, SCS TSG');
%     title('STRATUS-2017 SST, Buoy vs SCS');
%     if printit
%         ppath = fullfile(path_plots,['STRATUS_2017_SST_Buoy_v_SCS',graphformat]);
%         print(graphdevice,ppath);
%     end

    % specific humidity
    figure; plot(jd_hr,qa1hr,'bx');
    hold on; plot(jd_hr,s15qa1hr_1,'r+',jd_hr,s15qa1hr_2,'ro');
    plot(jd_hr,s16qa1hr_1,'g+',jd_hr,s16qa1hr_2,'go');
    grid; xlim([stjd endjd]); xlabel('DOY, 2017'); ylabel('Sp Hum g/kg');
    legend('PSD','S15-1','S15-2','S16-1','S16-2','location','south');
    title('STRATUS-2017 Humidity Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_Qair',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 10; maxx = 12;
    figure; plot(qa1hr_3m(ff),s15qa1hr_1(ff),'r+',qa1hr_3m(ff),s15qa1hr_2(ff),'ro');
    hold on; plot(qa1hr_3m(gg),s16qa1hr_1(gg),'b+',qa1hr_3m(gg),s16qa1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('S15-1','S15-2','S16-1','S16-2','location','northwest');
    ylabel('q g/kg, Buoy'); xlabel('q g/kg @ 3m, PSD');
    title('STRATUS-2017 Specific Humidity, Buoy vs PSD Qref @ 3m');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_Qair_v_Qref',graphformat]);
        print(graphdevice,ppath);
    end

    % RH
    figure; plot(jd_hr,rh1hr,'bx',jd_hr,imrh1hr,'kx');
    hold on; plot(jd_hr,s15rh1hr_1,'r+',jd_hr,s15rh1hr_2,'ro');
    plot(jd_hr,s16rh1hr_1,'go',jd_hr,s16rh1hr_2,'g+');
    grid; xlim([stjd endjd]); xlabel('DOY, 2017'); ylabel('RH %');
    legend('PSD','SCS','S15-1','S15-2','S16-1','S16-2','location','southeast');
    title('STRATUS-2017 RH Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_RH',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 63; maxx = 83;
    figure; plot(rh1hr_3m(ff),s15rh1hr_1(ff),'r+',rh1hr_3m(ff),s15rh1hr_2(ff),'ro');
    hold on; plot(rh1hr_3m(gg),s16rh1hr_1(gg),'b+',rh1hr_3m(gg),s16rh1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('S15-1','S15-2','S16-1','S16-2','location','northwest');
    ylabel('RH, Buoy'); xlabel('RH @ 3m, PSD');
    title('STRATUS-2017 RH, Buoy vs PSD RHref @ 3m');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_RH_v_RHref',graphformat]);
        print(graphdevice,ppath);
    end

    % pressure
    minx = min(Pmb1hr)-1; maxx = max(Pmb1hr)+1;
    figure; plot(jd_hr,Pmb1hr,'bx',jd_hr,imPmb1hr,'kx'); hold on;
    plot(jd_hr,s15Pmb1hr_1,'r+',jd_hr,s15Pmb1hr_2,'ro');
    plot(jd_hr,s16Pmb1hr_1,'g+',jd_hr,s16Pmb1hr_2,'go');
    grid; axis([stjd,endjd,minx,maxx]); xlabel('DOY, 2017'); ylabel('P mb @ sea level');
    legend('PSD','scs','S15-1','S15-2','S16-1','S16-2','location','south');
    title('STRATUS-2017 Pressure Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_Pressure',graphformat]);
        print(graphdevice,ppath);
    end
    
    minx = 1016; maxx = 1021;
    figure; plot(Pmb1hr(ff),s15Pmb1hr_1(ff),'r+',Pmb1hr(ff),s15Pmb1hr_2(ff),'ro'); hold on;
    plot(Pmb1hr(gg),s16Pmb1hr_1(gg),'g+',Pmb1hr(gg),s16Pmb1hr_2(gg),'go');
    hold on; plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('S15-1','S15-2','S16-1','S16-2','location','northwest');
    ylabel('P, mb, Buoy'); xlabel('P, mb, PSD');
    title('STRATUS-2017 Sea Level Pressure, Buoy vs PSD');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_Pressure_v_Pref',graphformat]);
        print(graphdevice,ppath);
    end

    % solar radiation
    figure; plot(jd_hr,rs1hr,'bx',jd_hr,imrs1hr,'bo');
    hold on; plot(jd_hr,s15rs1hr_1,'r+',jd_hr,s15rs1hr_2,'ro');
    plot(jd_hr,s16rs1hr_1,'g+',jd_hr,s16rs1hr_2,'go');
    grid; axis([stjd endjd -20 900]); xlabel('DOY, 2017'); ylabel('Solar Downwelling, W/m^2');
    legend('PSD','SCS','S15-1','S15-2','S16-1','S16-2','location','northwest');
    title('STRATUS-2017 Solar Radiation Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_PSP',graphformat]);
        print(graphdevice,ppath);
    end
    
    xx = find(mod(jd_hr,1)<=0.95 & mod(jd_hr,1)>=0.5); % daylight hours only
    yy = intersect(xx,bb);
    zz = intersect(xx,cc);
    minx = min(rs1hr(yy)); maxx = max(rs1hr(yy));
    figure; plot(rs1hr(xx),imrs1hr(xx),'kd',rs1hr(yy),s15rs1hr_1(yy),'r+',rs1hr(yy),s15rs1hr_2(yy),'ro');
    hold on; plot(rs1hr(zz),s16rs1hr_1(zz),'b+',rs1hr(zz),s16rs1hr_2(zz),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('SCS','S15-1','S15-2','S16-1','S16-2','location','northwest');
    ylabel('PSP W/m^2, Buoy'); xlabel('PSP W/m^2, PSD');
    title('STRATUS-2017 Solar Radiation, Buoy vs PSD, daylight hours');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_PSP_Buoy_v_PSD',graphformat]);
        print(graphdevice,ppath);
    end

    % infrared radiation
    figure; plot(jd_hr,rl1hr,'bx',jd_hr,imrl1hr,'bo');
    hold on; plot(jd_hr,s15rl1hr_1,'r+',jd_hr,s15rl1hr_2,'ro');
    plot(jd_hr,s16rl1hr_1,'g+',jd_hr,s16rl1hr_2,'go');
    grid; axis([stjd endjd 320 420]); xlabel('DOY, 2017'); ylabel('IR Downwelling, W/m^2');
    legend('PSD','SCS','S15-1','S15-2','S16-1','S16-2','location','southwest');
    title('STRATUS-2017 Infrared Radiation Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_PIR',graphformat]);
        print(graphdevice,ppath);
    end
    
    minx = min(rl1hr)-10; maxx = max(rl1hr)+10;
    figure; plot(rl1hr,imrl1hr,'kd',rl1hr(bb),s15rl1hr_1(bb),'r+',rl1hr(bb),s15rl1hr_2(bb),'ro');
    hold on; plot(rl1hr(cc),s16rl1hr_1(cc),'b+',rl1hr(cc),s16rl1hr_2(cc),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('SCS','S15-1','S15-2','S16-1','S16-2','location','northwest');
    ylabel('PIR W/m^2, Buoy'); xlabel('PIR W/m^2, PSD');
    title('STRATUS-2017 IR Radiation, Buoy vs PSD');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_PIR_Buoy_v_PSD',graphformat]);
        print(graphdevice,ppath);
    end
    
    % analyze PIR bias
    tod = mod(jd_hr,1)*24;
    [tod_bin,imir_diff_bin,~,~,~,~,~] = binave(tod,imrl1hr-rl1hr,2,0,24);
    [~,S15_1_diff_bin,~,~,~,~,~] = binave(tod(bb),s15rl1hr_1(bb)-rl1hr(bb),2,0,24);
    [~,S15_2_diff_bin,~,~,~,~,~] = binave(tod(bb),s15rl1hr_2(bb)-rl1hr(bb),2,0,24);
    [~,S16_1_diff_bin,~,~,~,~,~] = binave(tod(cc),s16rl1hr_1(cc)-rl1hr(cc),2,0,24);
    [~,S16_2_diff_bin,~,~,~,~,~] = binave(tod(cc),s16rl1hr_2(cc)-rl1hr(cc),2,0,24);
    tod_bin = tod_bin+1;
    
    figure; plot(tod_bin,imir_diff_bin,tod_bin,S15_1_diff_bin,tod_bin,S15_2_diff_bin); hold on;
    plot(tod_bin,S16_1_diff_bin,tod_bin,S16_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('SCS','S15-1','S15-2','S16-1','S16-2','location','south');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_PIR_Bias',graphformat]);
        print(graphdevice,ppath);
    end
    
    figure; plot(rs1hr(xx),imrl1hr(xx)-rl1hr(xx),'ro'); hold on;
    plot(rs1hr(yy),s15rl1hr_1(yy)-rl1hr(yy),'bx',rs1hr(yy),s15rl1hr_2(yy)-rl1hr(yy),'gx');
    plot(rs1hr(zz),s16rl1hr_1(zz)-rl1hr(zz),'rx',rs1hr(zz),s16rl1hr_2(zz)-rl1hr(zz),'kx');
    plot([0,900],[0,0],'k-'); grid;
    xlabel('solar flux, W/m^2, daytime only'); ylabel('Rl_{scs} - Rl_{psd}, W/m^2');
    legend('SCS','S15-1','S15-2','S16-1','S16-2','location','best'); title('Daytime Rl Bias');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2017_PIR_Bias_v_solar',graphformat]);
        print(graphdevice,ppath);
    end
    
end
