% STRATUS_2019_buoy_comparison.m
% Sergio 
% 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_stratus = 'D:\Data\STRATUS_2019\CabodeHornos\SCS\S18\S18data\';
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 = 113;
endjd = 115;
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
%{
    flux_hr 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
flux_hr = NaN(Ndays*24,33);
flux_hr(:,1) = jd_ref_hr;

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

%copy into flux_hr array
%flux_hr = temp;
% copy into flux_hr 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
         flux_hr(idx,2:33) = temp(jj,2:33);
     end
 end
clear temp;


%% load Cabo de Hornos data
%% no pasa
%{
    scsm array: nnn x 19
    1  jd       decimal day of year
    2  latm     latitude
    3  lonm     longitude
    4  lrgm     heading
    5  cogm     COG
    6  sogm     SOG
    7  imrh     RH, %
    8  imta     T air, C
    9  imum     true wind speed, m/s
    10 imdm     true wind direction
    11 impress  pressure, mb @ 11.5m
    12 imrwspd  rel wind speed, m/s
    13 imrwdir  rel wind direction
    14 tsgm     SST, thermosalinograph
    15 tssm     salinity, thermosalinograph
    16 urel     relative wind, U vector, m/s
    17 vrel     relative wind, V vector, m/s
    18 Ngps     N component of ship's speed
    19 Egps     E component of ship's speed
    20 utrue    true wind, U vector, m/s
    21 vtrue    true wind, V vector, m/s
%}

% create empty 1-min array for entire day range
%scsm = NaN(Ndays*1440,22);
%scsm(:,1) = jd_ref_1;

% load daily netcdf files and copy data into scsm array
%for ddd=stjd:endjd
 %   [MM,DD] = yd2md(2019,ddd);
  %  fnames = dir(fullfile(path_samos,['WTEY_2019',sprintf('%02i',MM),...
%        sprintf('%02i',DD),'*.nc']));  % e.g.WTEY_20170726v10001.nc
%    if isempty(fnames)
%        continue
%    else
%        temp = read_samos_nc_STRATUS_2019(ddd,path_samos,fnames);
%        for jj = 1:length(temp)
%            diff = jd_ref_1 - temp(jj,1);
%            [val,idx] = min(abs(diff));
%            if val<double(2.0/86400) && length(idx)==1
%                scsm(idx,2:16) = temp(jj,2:16);
%            end
%        end
%    end
%end

% add U/V relative wind component vectors to array
% [scsm(:,17),scsm(:,18)] = sd_to_uv(scsm(:,12),scsm(:,13));
% 
% % add Ngps,Egps ship speed vectors to array
% scsm(:,19) = scsm(:,6).*cos(scsm(:,5)./180.*pi); % Ngps
% scsm(:,20) = scsm(:,6).*sin(scsm(:,5)./180.*pi); % Egps
% 
% % add U/V true wind component vectors to array
% [scsm(:,21),scsm(:,22)] = sd_to_uv(scsm(:,9),scsm(:,10));
% 
% % average to 1-hr timebase
% scsm_hr = interval_avg(scsm(:,1),scsm(:,2:end),jd_ref_hr_bin);
% 
% %kill initial jd 206 data
% ii = jd_ref_hr<115.5;
% scsm_hr(ii,2:end) = NaN;

% %recompute COG from 1-hr Egps,Ngps
% scsm_hr(:,5) = atan2(scsm_hr(:,20), scsm_hr(:,19))*180/pi;
% scsm_hr(:,5) = mod(scsm_hr(:,5) + 360, 360);
% 
% %recompute 1-hr wind directions from UV components
%  [~,scsm_hr(:,13)] = uv_to_sd(scsm_hr(:,17),scsm_hr(:,18));    % rwdir
%  [~,scsm_hr(:,10)] = uv_to_sd(scsm_hr(:,21),scsm_hr(:,22));    % rdir
% % 
%  %rel wind back to +/- 180 format
%  ii = scsm(:,13)>180; scsm(ii,13) = scsm(ii,13)-360;
%  ii = scsm_hr(:,13)>180; scsm_hr(ii,13) = scsm_hr(ii,13)-360;

%% load buoy data

% load STRATUS-17 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

% DEFINE EMPTY ARRAY FOR W17 LOGGER 1
  W17_1_hr_txt = [];
 %W17_1_hr = NaN(126,33);
 W17_1_hr = NaN(length(flux_hr),33);
 W17_1_hr_txt = NaN(length(flux_hr),33);
 W17_1_hr(:,1) = jd_ref_hr;
 
 W17_2_hr = load(fullfile(path_stratus,'Stratus18_MET_sys2_nohead.txt'));
 %W17_2_hr = NaN(126,33);
 W17_2_hr(:,5) = datenum(W17_2_hr(:,1),W17_2_hr(:,2),W17_2_hr(:,3),W17_2_hr(:,4),0,0);
 ii = W17_2_hr(:,5)<=stjd+datenum(2018,12,31) & W17_2_hr(:,5)<endjd+datenum(2018,12,31)+1;
 W17_2_hr = [W17_2_hr(:,5:15),W17_2_hr(:,17:19)];

 W17_2_hr = W17_2_hr(ii,:);              % select cruise days
 W17_2_hr(:,8) = W17_2_hr(:,8)+180;      % convert wdir to 'from' convention
 W17_2_hr(:,2) = W17_2_hr(:,2)+0.125*3;  % sea level pressure
 W17_2_hr(:,8) = mod(W17_2_hr(:,8),360); % wdir
 %W17_1_hr(:,5) = W17_1_hr_txt(:,5) - datenum(2018,12,31);  % convert datenum to jd
 %W17_2_hr(:,1) = W17_2_hr_txt(:,5)
 % pad end of array to endjd
 ii = find(abs(jd_ref_hr-W17_1_hr(end,1))<2/86400);
 N = length(jd_ref_hr) - ii;
 N = length(jd_ref_hr);
 W17_2_hr = [W17_2_hr;NaN(N,14)];
 %W17_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
W18_1_hr_txt = [];
W18_1_hr_txt= load(fullfile(path_stratus,'Stratus18_MET_sys2_nohead.txt'));
W18_1_hr = NaN(126,33);
W18_1_hr(:,5) = datenum(W18_1_hr_txt(:,1),W18_1_hr_txt(:,2),W18_1_hr_txt(:,3),W18_1_hr_txt(:,4),0,0);
ii = W18_1_hr(:,5)<=stjd+datenum(2018,12,31) & W18_1_hr(:,5)<endjd+datenum(2018,12,31)+1;
W18_1_hr = [W18_1_hr_txt(:,5:15),W18_1_hr_txt(:,17:19)];
W18_1_hr(:,8) = W18_1_hr(:,8)+180;      % convert wdir to 'from' convention
W18_1_hr(:,2) = W18_1_hr(:,2)+0.125*3;  % sea level pressure
W18_1_hr(:,8) = mod(W18_1_hr(:,8),360); % wdir

W18_2_hr_txt = [];
W18_2_hr_txt= load(fullfile(path_stratus,'Stratus18_MET_sys2_nohead.txt'));
%W18_1_hr = NaN(126,33);
W18_2_hr(:,5) = datenum(W18_2_hr_txt(:,1),W18_2_hr_txt(:,2),W18_2_hr_txt(:,3),W18_2_hr_txt(:,4),0,0);
 ii = W18_2_hr(:,5)<=stjd+datenum(2018,12,31) & W18_2_hr(:,5)<endjd+datenum(2018,12,31)+1;
 W18_2_hr = [W18_2_hr_txt(:,5:15),W18_2_hr_txt(:,17:19)];
W18_2_hr(:,8) = W18_2_hr(:,8)+180;      % convert wdir to 'from' convention
W18_2_hr(:,2) = W18_2_hr(:,2)+0.125*3;  % sea level pressure
W18_2_hr(:,8) = mod(W18_2_hr(:,8),360); % wdir

%% define variables
% PSD
lat1Hr = flux_hr(:,18);
lon1Hr = flux_hr(:,19);
sog1Hr = flux_hr(:,14);
rwdir1hr = flux_hr(:,17);
rwdir1hr_std = flux_hr(:,33);
spd1hr = flux_hr(:,2);
spd1hr_3m = flux_hr(:,26);
dir1hr = flux_hr(:,3);
ta1hr = flux_hr(:,7);
ta1hr_3m = flux_hr(:,27);
ts1hr = flux_hr(:,4);
rh1hr = flux_hr(:,30);
rh1hr_3m = flux_hr(:,28);
rs1hr = flux_hr(:,10);
rl1hr = flux_hr(:,11);
Pmb1hr = flux_hr(:,29);
qa1hr = flux_hr(:,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]);

% W17
W17spd1hr_1 = W17_2_hr(:,7);
W17dir1hr_1 = W17_2_hr(:,8);
W17ta1hr_1 = W17_2_hr(:,3);
W17ts1hr_1 = W17_2_hr(:,4);
W17rh1hr_1 = W17_2_hr(:,9);
W17rs1hr_1 = W17_2_hr(:,10);
W17rl1hr_1 = W17_2_hr(:,11);
W17Pmb1hr_1 = W17_2_hr(:,2);
W17qa1hr_1 = qair_p([W17ta1hr_1,W17rh1hr_1,W17Pmb1hr_1]);

W17spd1hr_2 = W17_2_hr(:,7);
W17dir1hr_2 = W17_2_hr(:,8);
W17ta1hr_2 = W17_2_hr(:,3);
W17ts1hr_2 = W17_2_hr(:,4);
W17rh1hr_2 = W17_2_hr(:,9);
W17rs1hr_2 = W17_2_hr(:,10);
W17rl1hr_2 = W17_2_hr(:,11);
W17Pmb1hr_2 = W17_2_hr(:,2);
W17qa1hr_2 = qair_p([W17ta1hr_2,W17rh1hr_2,W17Pmb1hr_2]);

% W18
W18spd1hr_1 = W18_1_hr(:,7);
W18dir1hr_1 = W18_1_hr(:,8);
W18ta1hr_1 = W18_1_hr(:,3);
W18ts1hr_1 = W18_1_hr(:,4);
W18rh1hr_1 = W18_1_hr(:,9);
W18rs1hr_1 = W18_1_hr(:,10);
W18rl1hr_1 = W18_1_hr(:,11);
W18Pmb1hr_1 = W18_1_hr(:,2);
W18qa1hr_1 = qair_p([W18ta1hr_1,W18rh1hr_1,W18Pmb1hr_1]);

W18spd1hr_2 = W18_2_hr(:,7);
W18dir1hr_2 = W18_2_hr(:,8);
W18ta1hr_2 = W18_2_hr(:,3);
W18ts1hr_2 = W18_2_hr(:,4);
W18rh1hr_2 = W18_2_hr(:,9);
W18rs1hr_2 = W18_2_hr(:,10);
W18rl1hr_2 = W18_2_hr(:,11);
W18Pmb1hr_2 = W18_2_hr(:,2);
W18qa1hr_2 = qair_p([W18ta1hr_2,W18rh1hr_2,W18Pmb1hr_2]);

%% define selectors

% define mooring locations and position tolerance
lat18 = -17.638; lon18 = 275.081;
lat17 = -22.461; lon17 = 274.356;
tol = 0.03; % ~ 1.8 N miles or 3.33 km

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

% near STRATUS-18 lat/lon selector:
bb = find(lat1Hr>=lat18-tol & lat1Hr<=lat18+tol & lon1Hr>=lon18-tol & lon1Hr<=lon18+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 STRATUS:
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];

W17_1 = [W17ta1hr_1,W17rh1hr_1,W17ts1hr_1,W17spd1hr_1,unwrap(W17dir1hr_1*pi/180),...
       W17rs1hr_1,W17rl1hr_1,W17Pmb1hr_1];
W17_2 = [W17ta1hr_2,W17rh1hr_2,W17ts1hr_2,W17spd1hr_2,unwrap(W17dir1hr_2*pi/180),...
       W17rs1hr_2,W17rl1hr_2,W17Pmb1hr_2];
%W18_1 = [W18ta1hr_1,W18rh1hr_1,W18ts1hr_1,W18spd1hr_1,unwrap(W18dir1hr_1*pi/180),W18rs1hr_1,W18rl1hr_1,W18Pmb1hr_1];
%W18_2 = [W18ta1hr_2,W18rh1hr_2,W18ts1hr_2,W18spd1hr_2,unwrap(W18dir1hr_2*pi/180),W18rs1hr_2,W18rl1hr_2,W18Pmb1hr_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_W17_1_all = nanmean1(psd(dd,:) - W17_1(dd,:));
% mean_psd_W17_1_all(jj) = nanmean1(psd(hh,jj) - W17_1(hh,jj));
% mean_psd_W17_1_all(5) = nanmean1(psd(hh,5) - W17_1(hh,5))*180/pi;
% med_psd_W17_1_all = nanmedian1(psd(dd,:) - W17_1(dd,:));
% med_psd_W17_1_all(jj) = nanmedian1(psd(hh,jj) - W17_1(hh,jj));
% med_psd_W17_1_all(5) = nanmedian1(psd(hh,5) - W17_1(hh,5))*180/pi;
% std_psd_W17_1_all = nanstd1(psd(dd,:) - W17_1(dd,:));
% std_psd_W17_1_all(jj) = nanstd1(psd(hh,jj) - W17_1(hh,jj));
% std_psd_W17_1_all(5) = nanstd1(psd(hh,5) - W17_1(hh,5))*180/pi;

mean_psd_W17_2_all = nanmean1(psd(dd,:) - W17_2(dd,:));
mean_psd_W17_2_all(jj) = nanmean1(psd(hh,jj) - W17_2(hh,jj));
mean_psd_W17_2_all(5) = nanmean1(psd(hh,5) - W17_2(hh,5))*180/pi;
med_psd_W17_2_all = nanmedian1(psd(dd,:) - W17_2(dd,:));
med_psd_W17_2_all(jj) = nanmedian1(psd(hh,jj) - W17_2(hh,jj));
med_psd_W17_2_all(5) = nanmedian1(psd(hh,5) - W17_2(hh,5))*180/pi;
std_psd_W17_2_all = nanstd1(psd(dd,:) - W17_2(dd,:));
std_psd_W17_2_all(jj) = nanstd1(psd(hh,jj) - W17_2(hh,jj));
std_psd_W17_2_all(5) = nanstd1(psd(hh,5) - W17_2(hh,5))*180/pi;

% psd vs STRATUS-16 stats for period near bouy
mean_psd_W17_1 = nanmean1(psd(aa,:) - W17_1(aa,:));
mean_psd_W17_1(jj) = nanmean1(psd(ff,jj) - W17_1(ff,jj));
mean_psd_W17_1(5) = nanmean1(psd(ff,5) - W17_1(ff,5))*180/pi;
med_psd_W17_1 = nanmedian1(psd(aa,:) - W17_1(aa,:));
med_psd_W17_1(jj) = nanmedian1(psd(ff,jj) - W17_1(ff,jj));
med_psd_W17_1(5) = nanmedian1(psd(ff,5) - W17_1(ff,5))*180/pi;
std_psd_W17_1 = nanstd1(psd(aa,:) - W17_1(aa,:));
std_psd_W17_1(jj) = nanstd1(psd(ff,jj) - W17_1(ff,jj));
std_psd_W17_1(5) = nanstd1(psd(ff,5) - W17_1(ff,5))*180/pi;

mean_psd_W17_2 = nanmean1(psd(aa,:) - W17_2(aa,:));
mean_psd_W17_2(jj) = nanmean1(psd(ff,jj) - W17_2(ff,jj));
mean_psd_W17_2(5) = nanmean1(psd(ff,5) - W17_2(ff,5))*180/pi;
med_psd_W17_2 = nanmedian1(psd(aa,:) - W17_2(aa,:));
med_psd_W17_2(jj) = nanmedian1(psd(ff,jj) - W17_2(ff,jj));
med_psd_W17_2(5) = nanmedian1(psd(ff,5) - W17_2(ff,5))*180/pi;
std_psd_W17_2 = nanstd1(psd(aa,:) - W17_2(aa,:));
std_psd_W17_2(jj) = nanstd1(psd(ff,jj) - W17_2(ff,jj));
std_psd_W17_2(5) = nanstd1(psd(ff,5) - W17_2(ff,5))*180/pi;
% 
% % psd vs STRATUS-18 stats for entire period at site
% mean_psd_W18_1_all = nanmean1(psd(dd,:) - W18_1(dd,:));
% mean_psd_W18_1_all(jj) = nanmean1(psd(hh,jj) - W18_1(hh,jj));
% mean_psd_W18_1_all(5) = nanmean1(psd(hh,5) - W18_1(hh,5))*180/pi;
% med_psd_W18_1_all = nanmedian1(psd(dd,:) - W18_1(dd,:));
% med_psd_W18_1_all(jj) = nanmedian1(psd(hh,jj) - W18_1(hh,jj));
% med_psd_W18_1_all(5) = nanmedian1(psd(hh,5) - W18_1(hh,5))*180/pi;
% std_psd_W18_1_all = nanstd1(psd(dd,:) - W18_1(dd,:));
% std_psd_W18_1_all(jj) = nanstd1(psd(hh,jj) - W18_1(hh,jj));
% std_psd_W18_1_all(5) = nanstd1(psd(hh,5) - W18_1(hh,5))*180/pi;
% 
% mean_psd_W18_2_all = nanmean1(psd(dd,:) - W18_2(dd,:));
% mean_psd_W18_2_all(jj) = nanmean1(psd(hh,jj) - W18_2(hh,jj));
% mean_psd_W18_2_all(5) = nanmean1(psd(hh,5) - W18_2(hh,5))*180/pi;
% med_psd_W18_2_all = nanmedian1(psd(dd,:) - W18_2(dd,:));
% med_psd_W18_2_all(jj) = nanmedian1(psd(hh,jj) - W18_2(hh,jj));
% med_psd_W18_2_all(5) = nanmedian1(psd(hh,5) - W18_2(hh,5))*180/pi;
% std_psd_W18_2_all = nanstd1(psd(dd,:) - W18_2(dd,:));
% std_psd_W18_2_all(jj) = nanstd1(psd(hh,jj) - W18_2(hh,jj));
% std_psd_W18_2_all(5) = nanstd1(psd(hh,5) - W18_2(hh,5))*180/pi;
% 
% % psd vs STRATUS-18 stats for period near bouy
% mean_psd_W18_1 = nanmean1(psd(aa,:) - W18_1(aa,:));
% mean_psd_W18_1(jj) = nanmean1(psd(gg,jj) - W18_1(gg,jj));
% mean_psd_W18_1(5) = nanmean1(psd(gg,5) - W18_1(gg,5))*180/pi;
% med_psd_W18_1 = nanmedian1(psd(aa,:) - W18_1(aa,:));
% med_psd_W18_1(jj) = nanmedian1(psd(gg,jj) - W18_1(gg,jj));
% med_psd_W18_1(5) = nanmedian1(psd(gg,5) - W18_1(gg,5))*180/pi;
% std_psd_W18_1 = nanstd1(psd(aa,:) - W18_1(aa,:));
% std_psd_W18_1(jj) = nanstd1(psd(gg,jj) - W18_1(gg,jj));
% std_psd_W18_1(5) = nanstd1(psd(gg,5) - W18_1(gg,5))*180/pi;
% 
% mean_psd_W18_2 = nanmean1(psd(aa,:) - W18_2(aa,:));
% mean_psd_W18_2(jj) = nanmean1(psd(gg,jj) - W18_2(gg,jj));
% mean_psd_W18_2(5) = nanmean1(psd(gg,5) - W18_2(gg,5))*180/pi;
% med_psd_W18_2 = nanmedian1(psd(aa,:) - W18_2(aa,:));
% med_psd_W18_2(jj) = nanmedian1(psd(gg,jj) - W18_2(gg,jj));
% med_psd_W18_2(5) = nanmedian1(psd(gg,5) - W18_2(gg,5))*180/pi;
% std_psd_W18_2 = nanstd1(psd(aa,:) - W18_2(aa,:));
% std_psd_W18_2(jj) = nanstd1(psd(gg,jj) - W18_2(gg,jj));
% std_psd_W18_2(5) = nanstd1(psd(gg,5) - W18_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_W17_1_all;med_psd_W17_1_all;std_psd_W17_1_all]; % psd - W17-1 for all
% Stats = [Stats;mean_psd_W17_2_all;med_psd_W17_2_all;std_psd_W17_2_all]; % psd - W17-2 for all
% Stats = [Stats;mean_psd_W18_1_all;med_psd_W18_1_all;std_psd_W18_1_all]; % psd - W18-1 for all
% Stats = [Stats;mean_psd_W18_2_all;med_psd_W18_2_all;std_psd_W18_2_all]; % psd - W18-2 for all
% Stats = [Stats;mean_psd_W17_1;med_psd_W17_1;std_psd_W17_1];             % psd - W17-1 near bouy
% Stats = [Stats;mean_psd_W17_2;med_psd_W17_2;std_psd_W17_2];             % psd - W17-2 near bouy
% Stats = [Stats;mean_psd_W18_1;med_psd_W18_1;std_psd_W18_1];             % psd - W18-1 near bouy
% Stats = [Stats;mean_psd_W18_2;med_psd_W18_2;std_psd_W18_2];             % psd - W18-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',lon17,lat17,'bx','MarkerSize',8,'LineWidth',1.5); grid;
    xlabel('Lon'); ylabel('Lat'); axis([lon17-1.5*tol,lon17+1.5*tol,lat17-1.5*tol,lat17+1.5*tol]);
    title('STRATUS-18 Intercomparison Points');
    plotit = true;
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_W17_Intercomp',graphformat]);
        print(graphdevice,ppath);
    end

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

    % wind speed
    %figure; plot(jd_ref_hr,spd1hr,'bo',jd_ref_hr,imspd1hr,'bx',jd_ref_hr,W17spd1hr_1,'r+');
    figure; 
    plot(jd_ref_hr,spd1hr,'bo',jd_ref_hr,W17spd1hr_1,'r+');
    %hold on; plot(jd_ref_hr,W17spd1hr_2,'ro',jd_ref_hr,W18spd1hr_1,'g+',jd_ref_hr,W18spd1hr_2,'go');
    hold on; plot(jd_ref_hr,W17spd1hr_2,'ro');
    
    grid; xlim([stjd endjd]); xlabel('DOY, 2018'); ylabel('U m/s @ z=3m');
    legend('PSD','SCS','W17-1','W17-2','W18-1','W18-2','location','northwest');
    title('STRATUS-2018 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),W17spd1hr_1(ff),'r+',spd1hr_3m(ff),W17spd1hr_2(ff),'ro');
    %hold on; plot(spd1hr_3m(gg),W18spd1hr_1(gg),'b+',spd1hr_3m(gg),W18spd1hr_2(gg),'bo');
    %hold on; plot(spd1hr_3m(gg),W18spd1hr_1(gg),'b+',spd1hr_3m(gg),W18spd1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); grid;
    legend('W17-1','W17-2','W18-1','W18-2','location','northwest');
    ylabel('U m/s, Buoy'); xlabel('U m/s @ 3m, PSD');
    title('STRATUS-2018 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,W17dir1hr_1,'r+');
   figure; plot(jd_ref_hr,dir1hr,'bx',jd_ref_hr,jd_ref_hr,W17dir1hr_1,'r+'); 
   % hold on; plot(jd_ref_hr,W17dir1hr_2,'ro',jd_ref_hr,W18dir1hr_1,'g+',jd_ref_hr,W18dir1hr_2,'go');
    hold on; plot(jd_ref_hr,W17dir1hr_1,'ro');
    grid; axis([stjd endjd 0 360]); ylabel('True Wind Direction'); xlabel('DOY, 2018');
    legend('PSD','SCS','W17-1','W17-2','W18-1','W18-2','location','northeast');
    title('STRATUS-2018 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),W17dir1hr_1(ff),'r+',dir1hr(ff),W17dir1hr_2(ff),'ro');
    hold on; plot(dir1hr(gg),W18dir1hr_1(gg),'b+',dir1hr(gg),W18dir1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); grid; axis([minx,maxx,minx,maxx]);
    legend('W17-1','W17-2','W18-1','W18-2','location','northwest');
    ylabel('Dir True, Buoy'); xlabel('Dir True, PSD');
    title('STRATUS-2018 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,W17ta1hr_1,'r+');
    figure; plot(jd_ref_hr,ta1hr,'bx',jd_ref_hr,W17ta1hr_1,'r+');
    hold on; plot(jd_ref_hr,W17ta1hr_2,'ro',jd_ref_hr,W18ta1hr_1,'g+',jd_ref_hr,W18ta1hr_2,'go');
    grid; axis([stjd,endjd,minx,maxx]); xlabel('DOY, 2018'); ylabel('Air Temp, C @ z=3m');
    legend('PSD','SCS','W17-1','W17-2','W18-1','W18-2','location','southwest');
    ylabel('T air, C'); xlabel('DOY, 2018');
    title('STRATUS-2018 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),W17ta1hr_1(ff),'r+',ta1hr_3m(ff),W17ta1hr_2(ff),'ro');
    hold on; plot(ta1hr_3m(gg),W18ta1hr_1(gg),'b+',ta1hr_3m(gg),W18ta1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); grid;
    legend('W17-1','W17-2','W18-1','W18-2','location','northwest');
    ylabel('T air C, Buoy'); xlabel('T air C @ 3m, PSD');
    title('STRATUS-2018 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,W17ts1hr_1,'r+');
    figure; plot(jd_ref_hr,ts1hr,'bx',jd_ref_hr,W17ts1hr_1,'r+');
    hold on; plot(jd_ref_hr,W17ts1hr_2,'ro',jd_ref_hr,W18ts1hr_1,'g+',jd_ref_hr,W18ts1hr_2,'go');
    grid; axis([stjd endjd minx maxx]); xlabel('DOY, 2018'); ylabel('Sea Temp, C');
    legend('tsnk','ts38','ts45','W17-1','W17-2','W18-1','W18-2','location','southeast');
    title('STRATUS-2018 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),W17ts1hr_1(aa),'r+',ts1hr(aa),W17ts1hr_2(aa),'ro');
    figure; plot(ts1hr(dd),ts1hr(aa),W17ts1hr_1(aa),'r+',ts1hr(aa),W17ts1hr_2(aa),'ro');
    hold on; plot(ts1hr(bb),W18ts1hr_1(bb),'b+',ts1hr(bb),W18ts1hr_2(bb),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('ts38','W17-1','W17-2','W18-1','W18-2','location','northwest');
    ylabel('SST C, Ship/Buoy'); xlabel('SST C, Sea Snake');
    title('STRATUS-2018 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,W17qa1hr_1,'r+',jd_ref_hr,W17qa1hr_2,'ro');
    plot(jd_ref_hr,W18qa1hr_1,'g+',jd_ref_hr,W18qa1hr_2,'go');
    grid; xlim([stjd endjd]); xlabel('DOY, 2018'); ylabel('Sp Hum g/kg');
    legend('PSD','W17-1','W17-2','W18-1','W18-2','location','south');
    title('STRATUS-2018 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),W17qa1hr_1(ff),'r+',qa1hr_3m(ff),W17qa1hr_2(ff),'ro');
    hold on; plot(qa1hr_3m(gg),W18qa1hr_1(gg),'b+',qa1hr_3m(gg),W18qa1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('W17-1','W17-2','W18-1','W18-2','location','northwest');
    ylabel('q g/kg, Buoy'); xlabel('q g/kg @ 3m, PSD');
    title('STRATUS-2018 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');
    figure; plot(jd_ref_hr,rh1hr,'bx',jd_ref_hr);
    hold on; plot(jd_ref_hr,W17rh1hr_1,'r+',jd_ref_hr,W17rh1hr_2,'ro');
    plot(jd_ref_hr,W18rh1hr_1,'go',jd_ref_hr,W18rh1hr_2,'g+');
    grid; xlim([stjd endjd]); xlabel('DOY, 2018'); ylabel('RH %');
    legend('PSD','SCS','W17-1','W17-2','W18-1','W18-2','location','northeast');
    title('STRATUS-2018 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),W17rh1hr_1(ff),'r+',rh1hr_3m(ff),W17rh1hr_2(ff),'ro');
    hold on; plot(rh1hr_3m(gg),W18rh1hr_1(gg),'b+',rh1hr_3m(gg),W18rh1hr_2(gg),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('W17-1','W17-2','W18-1','W18-2','location','northwest');
    ylabel('RH, Buoy'); xlabel('RH @ 3m, PSD');
    title('STRATUS-2018 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;
    figure; plot(jd_ref_hr,Pmb1hr,'bx',jd_ref_hr); hold on;
    plot(jd_ref_hr,W17Pmb1hr_1,'r+',jd_ref_hr,W17Pmb1hr_2,'ro');
    plot(jd_ref_hr,W18Pmb1hr_1,'g+',jd_ref_hr,W18Pmb1hr_2,'go');
    grid; axis([stjd,endjd,minx,maxx]); xlabel('DOY, 2018'); ylabel('P mb @ sea level');
    legend('PSD','SCS','W17-1','W17-2','W18-1','W18-2','location','north');
    title('STRATUS-2018 Pressure Comparison');
    if printit
        ppath = fullfile(path_plots,['STRATUS_2019_Pressure',graphformat]);
        print(graphdevice,ppath);
    end

    minx = 1014; maxx = 1021;
    figure; plot(Pmb1hr(ff),W17Pmb1hr_1(ff),'r+',Pmb1hr(ff),W17Pmb1hr_2(ff),'ro'); hold on;
    plot(Pmb1hr(gg),W18Pmb1hr_1(gg),'g+',Pmb1hr(gg),W18Pmb1hr_2(gg),'go');
    hold on; plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('W17-1','W17-2','W18-1','W18-2','location','northwest');
    ylabel('P, mb, Buoy'); xlabel('P, mb, PSD');
    title('STRATUS-2018 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,W17rs1hr_1,'r+',jd_ref_hr,W17rs1hr_2,'ro');
    plot(jd_ref_hr,W18rs1hr_1,'g+',jd_ref_hr,W18rs1hr_2,'go');
    grid; axis([stjd endjd minx maxx]); xlabel('DOY, 2018'); ylabel('Solar Downwelling, W/m^2');
    legend('PSD','W17-1','W17-2','W18-1','W18-2','location','northwest');
    title('STRATUS-2018 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); % W17
    zz = intersect(xx,bb); % W18
    minx = min(rs1hr(yy)); maxx = max(rs1hr(yy));
    figure; plot(rs1hr(yy),W17rs1hr_1(yy),'r+',rs1hr(yy),W17rs1hr_2(yy),'ro');
    hold on; plot(rs1hr(zz),W18rs1hr_1(zz),'b+',rs1hr(zz),W18rs1hr_2(zz),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('W17-1','W17-2','W18-1','W18-2','location','northwest');
    ylabel('PSP W/m^2, Buoy'); xlabel('PSP W/m^2, PSD');
    title('STRATUS-2018 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,W17rl1hr_1,'r+',jd_ref_hr,W17rl1hr_2,'ro');
    plot(jd_ref_hr,W18rl1hr_1,'g+',jd_ref_hr,W18rl1hr_2,'go');
    grid; axis([stjd endjd minx maxx]); xlabel('DOY, 2018'); ylabel('IR Downwelling, W/m^2');
    legend('PSD','W17-1','W17-2','W18-1','W18-2','location','southwest');
    title('STRATUS-2018 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),W17rl1hr_1(aa),'r+',rl1hr(aa),W17rl1hr_2(aa),'ro');
    hold on; plot(rl1hr(bb),W18rl1hr_1(bb),'b+',rl1hr(bb),W18rl1hr_2(bb),'bo');
    plot([minx,maxx],[minx,maxx],'k-'); axis([minx,maxx,minx,maxx]); grid;
    legend('W17-1','W17-2','W18-1','W18-2','location','northwest');
    ylabel('PIR W/m^2, Buoy'); xlabel('PIR W/m^2, PSD');
    title('STRATUS-2018 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,W17_1_diff_bin,~,~,~,~,~] = binave(tod(aa),W17rl1hr_1(aa)-rl1hr(aa),2,0,24);
    [~,W17_2_diff_bin,~,~,~,~,~] = binave(tod(aa),W17rl1hr_2(aa)-rl1hr(aa),2,0,24);
    [~,W18_1_diff_bin,~,~,~,~,~] = binave(tod(bb),W18rl1hr_1(bb)-rl1hr(bb),2,0,24);
    [~,W18_2_diff_bin,~,~,~,~,~] = binave(tod(bb),W18rl1hr_2(bb)-rl1hr(bb),2,0,24);
    tod_bin = tod_bin+1;

    figure; plot(tod_bin,W17_1_diff_bin,tod_bin,W17_2_diff_bin); hold on;
    plot(tod_bin,W18_1_diff_bin,tod_bin,W18_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('W17-1','W17-2','W18-1','W18-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),W17rl1hr_1(yy)-rl1hr(yy),'bx',rs1hr(yy),W17rl1hr_2(yy)-rl1hr(yy),'gx'); hold on;
    plot(rs1hr(zz),W18rl1hr_1(zz)-rl1hr(zz),'rx',rs1hr(zz),W18rl1hr_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('W17-1','W17-2','W18-1','W18-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

