% Analysis and figures for manuscript submitted to Atmospheric Measurement
% Techniques on the effects of surface winds on MASC measurements
%
% ---------------------------------
% Kyle E. Fitch, Maj, USAF (AFIT/CI)
% University of Utah
% Salt Lake City, Utah
% kyle.fitch@us.af.mil
% 7-Nov-2019  --  Version 0.2
% ---------------------------------

% Dependencies:
% 	tight_subplot.m
% 	clipLabEnds.m
% 	calcLocMin2.m
% 	calcLinFit.m
%  isMyMax.m
% 	oli_data_windeffects.mat


% Close all figures
clear;
close all;

% set path
mypwd = pwd;
if strcmp('C:\Users\Kyle',mypwd(1:13))
    myhomedir = 'C:\Users\Kyle\OneDrive\Documents\PhD\';
    myslash = '\';
elseif strcmp('C:\Users\fitc',mypwd(1:13))
    myhomedir = 'C:\Users\fitch\OneDrive\Documents\PhD\';
    myslash = '\';
elseif strcmp('/Users/kfitch/D',mypwd(1:15))
    myhomedir = '/Users/kfitch/Desktop/we_code/plots/';
    myslash = '/';
elseif strcmp('/Users/kfitch',mypwd(1:13))
    myhomedir = '/Users/kfitch/OneDrive/Documents/PhD/';
    myslash = '/';
end

% datasets compiled with loadstats_thingr.m
myvars_path = 'oli_data_windeffects'; % name of MASC & ARM .mat data file

rim_thresh = [1.35,2]; % complexity thresholds for determining particle type
sn_temp_thresh = 2; % temperature used to further filter out rain (degC)
fs_thresh = 10; % fallspeeds above this value are assumed to be erroneous & omitted (m/s)

% RGB colors for plotting
mybrown = [0.45 0.24 0.10];

% wind-restricted categories
wind_restr = [5,1.50]; % (m s^-1)
wind_restr2 = [1.00,0.50]; % (m s^-1)

% where to check for bimodal dist.
mylow = 0;
myhigh = 1.1; 

%% Load output
load(myvars_path);
   
%% set plot parameters
%output orientation and size
print_orientation = 'landscape'; % landscape or portrait

if strcmp(print_orientation,'landscape') == 1
        pagewidth = 11; pageheight = 8.5;
else
        pagewidth = 8.5; pageheight = 11;
end

% set figure directory path
fig_dir = [myhomedir,'figures'];

% standard font size for figures
text_size1 = 22;

mypts_fs = linspace(0,10,1e3); % points at which to evaluate fall speed PDF
mypts_theta = linspace(0,90,1e3); % points at which to evaluate orientation PDF
mypts_Usfc = linspace(0,19,1e3); % points at which to evaluate wind speed PDF

%% standardize variable names
% particle size
dmax = d_max; % maximum dimension (mm)

% particle shape
chi = oli_complexity; % complexity
alpha = aspect_ratio; % aspect ratio

% particle falling 
theta = orientation; % fall orientation (deg)
v = fallspeed; % fall speed (m s^-1)

% mean Doppler velocity (m s^-1) below cloud base, 10-min avg, offset by 5 min 
% to allow for hydrometeor settling time
mdv(mdv>0) = nan; % set positive (upward) Doppler velocities to nan
vd = -mdv; 

% surface variables 
T_sfc = met_temp; % surface temperature (C)
U_sfc = met_wspd; % surface wind speed (m s^-1)
p_sfc = met_press; % surface pressure (kPa)
T_sfc(T_sfc==-9999) = nan; % set missing values to nan
U_sfc(U_sfc==-9999) = nan; % set missing values to nan
p_sfc(p_sfc==-9999) = nan; % set missing values to nan

% date of wind shield installation at Oliktok Point
fence_start = [2016 08 22 00 00 00];

% datetime(s) of erroneous particles
bad_hr = [2017 06 20 23 00 00];

% find index ranges for three primary riming classes
aggregaterange = find(oli_complexity > rim_thresh(2));
rimedrange = find(oli_complexity <= rim_thresh(2) & oli_complexity > rim_thresh(1));
graupelrange = find(oli_complexity <= rim_thresh(1));

% date/time should be a double for use in calculations
masc_datetime2 = double(masc_datetime);

% index of particles within temp/fallspeed thresholds, avoiding bad
% datetimes
basic_inds = find(T_sfc < sn_temp_thresh & v < fs_thresh ...
    & (datenum(masc_datetime2) < datenum(bad_hr) | ...
    datenum(masc_datetime2) >= datenum(bad_hr)+1/24) );

% wind shield dates with non-nan Doppler velocity data
WS_inds = intersect(basic_inds,...
    find(datenum(masc_datetime2) > datenum(fence_start) & ~isnan(vd))); 

% non-wind-shield dates with non-nan Doppler velocity data
NWS_inds = intersect(basic_inds,...
    find(datenum(masc_datetime2) < datenum(fence_start) & ~isnan(vd))); 

%% Compare fall speeds & PDFs
ymax = 3; % max value for y-axis (to standardize)
xmax = 1.7; % max value for x-axis (to standardize)

h01 = figure;
set(h01,'units','inches','position',[0 0 pagewidth pageheight]);
h01a = tight_subplot(1,1,0,[0.105,0.02],[0.08,0.02]);
axes(h01a(1))
    % find values from kernel density calculation using pts specified above
    [N_v1,xloc_v1] = ksdensity(v(NWS_inds),mypts_fs);
    [N_vd1,xloc_vd1] = ksdensity(vd(NWS_inds),mypts_fs);
    [N_v4,xloc_v4] = ksdensity(v(WS_inds),mypts_fs);
    [N_vd4,xloc_vd4] = ksdensity(vd(WS_inds),mypts_fs);
    p01a = plot(xloc_v1,N_v1,'k--','LineWidth',2.5);
    hold on
    p02a = plot(xloc_vd1,N_vd1,'r--','LineWidth',2.5);
    p01b = plot(xloc_v4,N_v4,'k-','LineWidth',2.5);
    p02b = plot(xloc_vd4,N_vd4,'r-','LineWidth',2.5);
    max_v1 = isMyMax(xloc_v1,N_v1,0,10)
    max_v4 = isMyMax(xloc_v4,N_v4,0,10)
    max_v4_2 = isMyMax(xloc_v4,N_v4,0.5,1.5)
    v4_locmin = calcLocMin2(N_v4,xloc_v4,mylow,myhigh);
    yvals = linspace(0,ymax,100);
    xvals = v4_locmin.*ones(size(yvals));
    plot(xvals,yvals,'k.-','LineWidth',1)
    
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    set(gca,'FontSize',text_size1+2)
    set(gca,'XTick',0:0.5:xmax)
    ylabel('\textbf{PDF $\mathbf{[m\,s^{-1}]^{-1}}$}','FontSize',text_size1+2,...
        'FontWeight','bold','Interpreter','Latex')
    xlabel('\textbf{Fall Speed ($v_p$) [$\mathbf{m\,s^{-1}}$]}',...
        'FontSize',text_size1+2,'FontWeight','bold','Interpreter','Latex')
    text(0.12,2.5,'\textbf{MASC}','Color','k','FontSize',text_size1+2,...
        'Interpreter','Latex');
    text(0.8,2.0,'\textbf{KAZR}','Color','r','FontSize',text_size1+2,...
        'Interpreter','Latex');
    plot(linspace(1.10,1.23,10),2.5*ones(1,10),'LineStyle','--',...
        'Color',[0.5 0.5 0.5],'LineWidth',2.5);
    plot(linspace(1.10,1.23,10),2.3*ones(1,10),'LineStyle','-',...
        'Color',[0.5 0.5 0.5],'LineWidth',2.5);
    text(1.25,2.5,'\textbf{unshielded}','Color',[0.5 0.5 0.5],'FontSize',text_size1+2,...
        'Interpreter','Latex');
    text(1.25,2.3,'\textbf{shielded}','Color',[0.5 0.5 0.5],'FontSize',text_size1+2,...
        'Interpreter','Latex');
    text(v4_locmin+0.01,2.7,['$v_c\mathbf{=',num2str(v4_locmin,...
        '%4.2f'),'\,\rm{m\,s^{-1}}}$'],'FontSize',text_size1+2,'Interpreter','Latex');
    
orient(h01,'landscape');
print([fig_dir,myslash,'v_pdf.pdf'],'-dpdf','-painters','-fillpage');

%% Delta PDF
lab_xpos = 0.85;
lab_ypos = 0.90;

h0 = figure;
set(h0,'units','inches','position',[0 0 6.5 pageheight]);
h0a = tight_subplot(2,1,0,[0.13 0.06],[0.20 0.06]);
axes(h0a(1))
    [N1a,xloc1a] = ksdensity(U_sfc(intersect(WS_inds,find(v>v4_locmin))),mypts_Usfc);
    [N1b,xloc1b] = ksdensity(U_sfc(intersect(WS_inds,find(v<=v4_locmin))),mypts_Usfc);
    pa = plot(mypts_Usfc,N1a,'k-','LineWidth',2);
    xlim([0 20])
    hold on
    pb = plot(mypts_Usfc,N1b,'Color',[0.4 0.4 0.4],'LineWidth',2,'LineStyle','--');
    grid on
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:2:20)
    xticklabels('');
    ylabel('\textbf{PDF $\mathbf{[m\,s^{-1}]^{-1}}$}','FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex')
    clipLabEnds('YTickLabels','low');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(a)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    text(0.1,0.14,['$\mathbf{p(}v_p\mathbf{>',num2str(v4_locmin,'%4.2f'),...
        '\,\rm{m\,s^{-1}})}$'],'FontSize',text_size1,'Interpreter','Latex')
    text(10,0.11,['$\mathbf{p(}v_p\mathbf{\leq',num2str(v4_locmin,'%4.2f'),...
        '\,\rm{m\,s^{-1}})}$'],'FontSize',text_size1,'Interpreter','Latex',...
        'Color',[0.4 0.4 0.4])
axes(h0a(2))
    [Na,xloca] = ksdensity(U_sfc(intersect(WS_inds,find(v>v4_locmin))),mypts_Usfc);
    [Nb,xlocb] = ksdensity(U_sfc(intersect(WS_inds,find(v<=v4_locmin))),mypts_Usfc);
    N_diff = Na-Nb;
    mysep = find(N_diff<0 & mypts_Usfc > 1);
    pa1 = plot(mypts_Usfc(1:mysep(1)),N_diff(1:mysep(1)),'b-','LineWidth',2);
    hold on;
    pa2 = plot(mypts_Usfc(mysep(1):end),N_diff(mysep(1):end),'r-','LineWidth',2);
    xlim([0 20])
    ylim([-0.05 0.09])
    hold on
    myx = linspace(0,20,100);
    myy = zeros(size(myx));
    plot(myx,myy,'k--','LineWidth',2);
    grid on
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:2:20)
    ylabel('\textbf{$\mathbf{\Delta}$PDF $\mathbf{[m\,s^{-1}]^{-1}}$}',...
        'FontSize',text_size1,'FontWeight','bold','Interpreter','Latex')
    xlabel('\textbf{Wind Speed ($U_{sfc}$) $\mathbf{[m\,s^{-1}]}$}',...
        'FontSize',text_size1,'FontWeight','bold','Interpreter','Latex')
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(b)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    text(2.3,0.080,['$\mathbf{\Delta PDF=p(}v_p\mathbf{>',num2str(v4_locmin,'%4.2f'),...
        '\,\rm{m\,s^{-1}})}$'],'FontSize',text_size1-1,'Interpreter','Latex')
    text(6,0.070,['$\mathbf{-p(}v_p\mathbf{\leq',num2str(v4_locmin,'%4.2f'),...
        '\,\rm{m\,s^{-1}})}$'],'FontSize',text_size1-1,'Interpreter','Latex',...
        'Color',[0.4 0.4 0.4])
    text(5,0.040,['$v_p\mathbf{>',num2str(v4_locmin,'%4.2f'),...
        '\,\rm{m\,s^{-1}}}$'],'FontSize',text_size1-1,'Interpreter','Latex',...
        'Color','b')
    text(5,0.030,'\textbf{more probable}','FontSize',text_size1-1,'Interpreter','Latex',...
        'Color','b')
    text(12,-0.030,['$v_p\mathbf{\leq',num2str(v4_locmin,'%4.2f'),...
        '\,\rm{m\,s^{-1}}}$'],'FontSize',text_size1-1,'Interpreter','Latex',...
        'Color','r')
    text(12,-0.040,'\textbf{more probable}','FontSize',text_size1-1,'Interpreter','Latex',...
        'Color','r')
orient(h0,'portrait');
print([fig_dir,myslash,'delta_pdf.pdf'],'-dpdf','-painters','-fillpage');

%% Table values of counts & percentages

inds1 = intersect(NWS_inds,find(U_sfc>wind_restr(1)));
inds2 = intersect(NWS_inds,find(U_sfc<=wind_restr(1)));
inds3 = intersect(NWS_inds,find(U_sfc<=wind_restr(2)));
inds4 = intersect(WS_inds,find(U_sfc>wind_restr(1)));
inds5 = intersect(WS_inds,find(U_sfc<=wind_restr(1)));
inds6 = intersect(WS_inds,find(U_sfc<=wind_restr(2)));
inds7 = intersect(WS_inds,find(U_sfc<=wind_restr2(1)));
inds8 = intersect(WS_inds,find(U_sfc<=wind_restr2(2)));
inds07 = intersect(NWS_inds,find(U_sfc<=wind_restr2(1)));
inds08 = intersect(NWS_inds,find(U_sfc<=wind_restr2(2)));

inds1a = intersect(inds1,aggregaterange);
inds1b = intersect(inds1,rimedrange);
inds1c = intersect(inds1,graupelrange);

inds2a = intersect(inds2,aggregaterange);
inds2b = intersect(inds2,rimedrange);
inds2c = intersect(inds2,graupelrange);

inds3a = intersect(inds3,aggregaterange);
inds3b = intersect(inds3,rimedrange);
inds3c = intersect(inds3,graupelrange);

inds4a = intersect(inds4,aggregaterange);
inds4b = intersect(inds4,rimedrange);
inds4c = intersect(inds4,graupelrange);

inds5a = intersect(inds5,aggregaterange);
inds5b = intersect(inds5,rimedrange);
inds5c = intersect(inds5,graupelrange);

inds6a = intersect(inds6,aggregaterange);
inds6b = intersect(inds6,rimedrange);
inds6c = intersect(inds6,graupelrange);

inds7a = intersect(inds7,aggregaterange);
inds7b = intersect(inds7,rimedrange);
inds7c = intersect(inds7,graupelrange);

inds8a = intersect(inds8,aggregaterange);
inds8b = intersect(inds8,rimedrange);
inds8c = intersect(inds8,graupelrange);

inds07a = intersect(inds07,aggregaterange);
inds07b = intersect(inds07,rimedrange);
inds07c = intersect(inds07,graupelrange);

inds08a = intersect(inds08,aggregaterange);
inds08b = intersect(inds08,rimedrange);
inds08c = intersect(inds08,graupelrange);

disp(['NWS, U_{sfc}>',num2str(wind_restr(1)),':',num2str(length(inds1))])
disp(['NWS, U_{sfc}\leq',num2str(wind_restr(1)),':',num2str(length(inds2))])
disp(['NWS, U_{sfc}\leq',num2str(wind_restr(2)),':',num2str(length(inds3))])
disp(['NWS, U_{sfc}\leq',num2str(wind_restr2(1)),':',num2str(length(inds07))])
disp(['NWS, U_{sfc}\leq',num2str(wind_restr2(2)),':',num2str(length(inds08))])

disp(['NWS Agg, U_{sfc}>',num2str(wind_restr(1)),':',num2str(length(inds1a)),'(',...
    num2str(round(100*length(inds1a)/length(inds1))),'%)'])
disp(['NWS Agg, U_{sfc}\leq',num2str(wind_restr(1)),':',num2str(length(inds2a)),'(',...
    num2str(round(100*length(inds2a)/length(inds2))),'%)'])
disp(['NWS Agg, U_{sfc}\leq',num2str(wind_restr(2)),':',num2str(length(inds3a)),'(',...
    num2str(round(100*length(inds3a)/length(inds3))),'%)'])
disp(['NWS Agg, U_{sfc}\leq',num2str(wind_restr2(1)),':',num2str(length(inds07a)),'(',...
    num2str(round(100*length(inds07a)/length(inds07))),'%)'])
disp(['NWS Agg, U_{sfc}\leq',num2str(wind_restr2(2)),':',num2str(length(inds08a)),'(',...
    num2str(round(100*length(inds08a)/length(inds08))),'%)'])

disp(['NWS MR, U_{sfc}>',num2str(wind_restr(1)),':',num2str(length(inds1b)),'(',...
    num2str(round(100*length(inds1b)/length(inds1))),'%)'])
disp(['NWS MR, U_{sfc}\leq',num2str(wind_restr(1)),':',num2str(length(inds2b)),'(',...
    num2str(round(100*length(inds2b)/length(inds2))),'%)'])
disp(['NWS MR, U_{sfc}\leq',num2str(wind_restr(2)),':',num2str(length(inds3b)),'(',...
    num2str(round(100*length(inds3b)/length(inds3))),'%)'])
disp(['NWS MR, U_{sfc}\leq',num2str(wind_restr2(1)),':',num2str(length(inds07b)),'(',...
    num2str(round(100*length(inds07b)/length(inds07))),'%)'])
disp(['NWS MR, U_{sfc}\leq',num2str(wind_restr2(2)),':',num2str(length(inds08b)),'(',...
    num2str(round(100*length(inds08b)/length(inds08))),'%)'])

disp(['NWS R, U_{sfc}>',num2str(wind_restr(1)),':',num2str(length(inds1c)),'(',...
    num2str(round(100*length(inds1c)/length(inds1))),'%)'])
disp(['NWS R, U_{sfc}\leq',num2str(wind_restr(1)),':',num2str(length(inds2c)),'(',...
    num2str(round(100*length(inds2c)/length(inds2))),'%)'])
disp(['NWS R, U_{sfc}\leq',num2str(wind_restr(2)),':',num2str(length(inds3c)),'(',...
    num2str(round(100*length(inds3c)/length(inds3))),'%)'])
disp(['NWS R, U_{sfc}\leq',num2str(wind_restr2(1)),':',num2str(length(inds07c)),'(',...
    num2str(round(100*length(inds07c)/length(inds07))),'%)'])
disp(['NWS R, U_{sfc}\leq',num2str(wind_restr2(2)),':',num2str(length(inds08c)),'(',...
    num2str(round(100*length(inds08c)/length(inds08))),'%)'])

disp(['WS, U_{sfc}>',num2str(wind_restr(1)),':',num2str(length(inds4))])
disp(['WS, U_{sfc}\leq',num2str(wind_restr(1)),':',num2str(length(inds5))])
disp(['WS, U_{sfc}\leq',num2str(wind_restr(2)),':',num2str(length(inds6))])
disp(['WS, U_{sfc}\leq',num2str(wind_restr2(1)),':',num2str(length(inds7))])
disp(['WS, U_{sfc}\leq',num2str(wind_restr2(2)),':',num2str(length(inds8))])

disp(['WS Agg, U_{sfc}>',num2str(wind_restr(1)),':',num2str(length(inds4a)),'(',...
    num2str(round(100*length(inds4a)/length(inds4))),'%)'])
disp(['WS Agg, U_{sfc}\leq',num2str(wind_restr(1)),':',num2str(length(inds5a)),'(',...
    num2str(round(100*length(inds5a)/length(inds5))),'%)'])
disp(['WS Agg, U_{sfc}\leq',num2str(wind_restr(2)),':',num2str(length(inds6a)),'(',...
    num2str(round(100*length(inds6a)/length(inds6))),'%)'])
disp(['WS Agg, U_{sfc}\leq',num2str(wind_restr2(1)),':',num2str(length(inds7a)),'(',...
    num2str(round(100*length(inds7a)/length(inds7))),'%)'])
disp(['WS Agg, U_{sfc}\leq',num2str(wind_restr2(2)),':',num2str(length(inds8a)),'(',...
    num2str(round(100*length(inds8a)/length(inds8))),'%)'])

disp(['WS MR, U_{sfc}>',num2str(wind_restr(1)),':',num2str(length(inds4b)),'(',...
    num2str(round(100*length(inds4b)/length(inds4))),'%)'])
disp(['WS MR, U_{sfc}\leq',num2str(wind_restr(1)),':',num2str(length(inds5b)),'(',...
    num2str(round(100*length(inds5b)/length(inds5))),'%)'])
disp(['WS MR, U_{sfc}\leq',num2str(wind_restr(2)),':',num2str(length(inds6b)),'(',...
    num2str(round(100*length(inds6b)/length(inds6))),'%)'])
disp(['WS MR, U_{sfc}\leq',num2str(wind_restr2(1)),':',num2str(length(inds7b)),'(',...
    num2str(round(100*length(inds7b)/length(inds7))),'%)'])
disp(['WS MR, U_{sfc}\leq',num2str(wind_restr2(2)),':',num2str(length(inds8b)),'(',...
    num2str(round(100*length(inds8b)/length(inds8))),'%)'])

disp(['WS R, U_{sfc}>',num2str(wind_restr(1)),':',num2str(length(inds4c)),'(',...
    num2str(round(100*length(inds4c)/length(inds4))),'%)'])
disp(['WS R, U_{sfc}\leq',num2str(wind_restr(1)),':',num2str(length(inds5c)),'(',...
    num2str(round(100*length(inds5c)/length(inds5))),'%)'])
disp(['WS R, U_{sfc}\leq',num2str(wind_restr(2)),':',num2str(length(inds6c)),'(',...
    num2str(round(100*length(inds6c)/length(inds6))),'%)'])
disp(['WS R, U_{sfc}\leq',num2str(wind_restr2(1)),':',num2str(length(inds7c)),'(',...
    num2str(round(100*length(inds7c)/length(inds7))),'%)'])
disp(['WS R, U_{sfc}\leq',num2str(wind_restr2(2)),':',num2str(length(inds8c)),'(',...
    num2str(round(100*length(inds8c)/length(inds8))),'%)'])
    
%% Compare MASC fall speed to KAZR Doppler Fall Speed
ymax = 6; % max value for y-axis (to standardize)
xmax = 1.8; % max value for x-axis (to standardize)

lab_xpos = 0.80;
lab_ypos = 0.90;

h1 = figure;
set(h1,'units','inches','position',[0 0 pagewidth pageheight]);
h1a = tight_subplot(2,3,0,[0.10 0.08],[0.16 0]);
axes(h1a(1))
    p11 = plot(xloc_v1,N_v1,'k-','LineWidth',2);
    hold on
    p12 = plot(xloc_vd1,N_vd1,'r-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    text(0.12,2.5,'\textbf{MASC}','Color','k','FontSize',text_size1,...
        'Interpreter','Latex');
    text(0.8,2.0,'\textbf{KAZR}','Color','r','FontSize',text_size1,...
        'Interpreter','Latex');
    yl1 = ylabel({'\textbf{Unshielded}','(Nov 2015 -- Aug 2016)'},'FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex');
    title({['$\mathbf{>',num2str(wind_restr(1),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:0.5:xmax)
    set(gca,'XTickLabels','')
    clipLabEnds('YTickLabels','both');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(a)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h1a(2))
    [N_v2,xloc_v2] = ksdensity(v(inds2),mypts_fs);
    [N_vd2,xloc_vd2] = ksdensity(vd(inds2),mypts_fs);
    p21 = plot(xloc_v2,N_v2,'k-','LineWidth',2);
    hold on
    p22 = plot(xloc_vd2,N_vd2,'r-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    yl2 = ylabel('\textbf{PDF $\mathbf{[m\,s^{-1}]^{-1}}$}','FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex');
    title({'$U_{sfc}\mathbf{\,[m\,s^{-1}]}$',['$\mathbf{\leq',num2str(wind_restr(1),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:0.5:xmax)
    set(gca,'XTickLabels','')
    set(gca,'YTickLabels','')
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(b)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h1a(3))
    [N_v3,xloc_v3] = ksdensity(v(inds3),mypts_fs);
    [N_vd3,xloc_vd3] = ksdensity(vd(inds3),mypts_fs);
    p31 = plot(xloc_v3,N_v3,'k-','LineWidth',2);
    hold on
    p32 = plot(xloc_vd3,N_vd3,'r-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    v3_locmin = calcLocMin2(N_v3,xloc_v3,mylow,myhigh);
    yvals3 = linspace(0,ymax,100);
    xvals3 = v3_locmin.*ones(size(yvals3));
    plot(xvals3,yvals3,'k.-','LineWidth',1)
    text(v3_locmin+0.01,myylim(1)+lab_ypos*(myylim(2)-myylim(1)),...
        '$v_c\mathbf{=}$','FontSize',text_size1+2,'Interpreter','latex')
    text(v3_locmin+0.01,myylim(1)+lab_ypos*(myylim(2)-myylim(1))-0.50,...
        ['$\mathbf{',num2str(v3_locmin,...
        '%4.2f'),'\,\rm{m\,s^{-1}}}$'],'FontSize',text_size1+2,'Interpreter','Latex')
    title({['$\mathbf{\leq',num2str(wind_restr(2),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:0.5:xmax)
    set(gca,'XTickLabels','')
    set(gca,'YTickLabels','')
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(c)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h1a(4))
    p41 = plot(xloc_v4,N_v4,'k-','LineWidth',2);
    hold on
    yvals4 = linspace(0,ymax,100);
    xvals4 = v4_locmin.*ones(size(yvals4));
    plot(xvals4,yvals4,'k.-','LineWidth',1)
    p42 = plot(xloc_vd4,N_vd4,'r-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    text(v4_locmin+0.01,3.5,...
        '$v_c\mathbf{=}$','FontSize',text_size1+2,'Interpreter','latex')
    text(v4_locmin+0.01,3.5-0.5,['$\mathbf{',num2str(v4_locmin,...
        '%4.2f'),'\,\rm{m\,s^{-1}}}$'],'FontSize',text_size1+2,'Interpreter','latex')
    xlabel('')
    yl4 = ylabel({'\textbf{Shielded}','(Sep 2016 -- Aug 2018)'},'FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex');
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:0.5:xmax)
    clipLabEnds('XTickLabels','low');
    clipLabEnds('YTickLabels','both');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(d)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h1a(5))
    [N_v5,xloc_v5] = ksdensity(v(inds5),mypts_fs);
    [N_vd5,xloc_vd5] = ksdensity(vd(inds5),mypts_fs);
    p51 = plot(xloc_v5,N_v5,'k-','LineWidth',2);
    hold on
    v5_locmin = calcLocMin2(N_v5,xloc_v5,mylow,myhigh);
    yvals5 = linspace(0,ymax,100);
    xvals5 = v5_locmin.*ones(size(yvals5));
    plot(xvals5,yvals5,'k.-','LineWidth',1)
    p52 = plot(xloc_vd5,N_vd5,'r-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    text(v5_locmin+0.01,3.5,...
        '$v_c\mathbf{=}$','FontSize',text_size1+2,'Interpreter','latex')
    text(v5_locmin+0.01,3.5-0.5,['$\mathbf{',num2str(v5_locmin,...
        '%4.2f'),'\,\rm{m\,s^{-1}}}$'],'FontSize',text_size1+2,'Interpreter','latex')
    xlabel('\textbf{Fall Speed $(v_p)\mathbf{\,[m\,s^{-1}]}$}',...
        'FontSize',text_size1,'FontWeight','bold','Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:0.5:xmax)
    set(gca,'YTickLabels','')
    clipLabEnds('XTickLabels','low');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(e)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h1a(6))
    [N_v6,xloc_v6] = ksdensity(v(inds6),mypts_fs);
    [N_vd6,xloc_vd6] = ksdensity(vd(inds6),mypts_fs);
    p61 = plot(xloc_v6,N_v6,'k-','LineWidth',2);
    hold on
    v6_locmin = calcLocMin2(N_v6,xloc_v6,mylow,myhigh);
    yvals6 = linspace(0,ymax,100);
    xvals6 = v6_locmin.*ones(size(yvals6));
    plot(xvals6,yvals6,'k.-','LineWidth',1)
    p62 = plot(xloc_vd6,N_vd6,'r-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    text(v6_locmin+0.01,3.5,...
        '$v_c\mathbf{=}$','FontSize',text_size1+2,'Interpreter','latex')
    text(v6_locmin+0.01,3.5-0.5,['$\mathbf{',num2str(v6_locmin,...
        '%4.2f'),'\,\rm{m\,s^{-1}}}$'],'FontSize',text_size1+2,'Interpreter','latex')
    xlabel('')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:0.5:xmax)
    set(gca,'YTickLabels','')
    clipLabEnds('XTickLabels','low');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(f)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    pct_below = 100*length(find(v(inds6)<v6_locmin))/length(inds6);
    
% center y-label
myyl2_x = yl2.Position(1);
yl2.Position(2) = 0;
yl2.Position(1) = myyl2_x - xmax - 0.8;

myyl1_y = yl1.Position(2);
yl1.Position(2) = myyl1_y + 0.2;

myyl4_y = yl4.Position(2);
yl4.Position(2) = myyl4_y - 0.2;
    
orient(h1,'landscape');
print([fig_dir,myslash,'v_vd_comp.pdf'],'-dpdf','-painters','-fillpage');

%% Calculate v_cutoff(U_sfc)

my_incr = 0.50;

wind_array = v4_locmin:my_incr:10+my_incr;

locmin = nan(size(wind_array));
locmin_agg = nan(size(wind_array));
locmin_mrim = nan(size(wind_array));
locmin_rim = nan(size(wind_array));
frac_agg = nan(size(wind_array));
frac_mrim = nan(size(wind_array));
frac_rim = nan(size(wind_array));

for m = 1:length(wind_array)
    myinds = intersect(WS_inds,find(U_sfc<=wind_array(m)));
    myinds_agg = intersect(myinds,aggregaterange);
    myinds_mrim = intersect(myinds,rimedrange);
    myinds_rim = intersect(myinds,graupelrange);
    frac_agg(m) = length(myinds_agg)/length(myinds);
    frac_mrim(m) = length(myinds_mrim)/length(myinds);
    frac_rim(m) = length(myinds_rim)/length(myinds);
    [N,xloc] = ksdensity(v(myinds),mypts_fs);
    [N_agg,xloc_agg] = ksdensity(v(myinds_agg),mypts_fs);
    [N_mrim,xloc_mrim] = ksdensity(v(myinds_mrim),mypts_fs);
    [N_rim,xloc_rim] = ksdensity(v(myinds_rim),mypts_fs);
    locmin(m) = calcLocMin2(N,xloc,mylow,myhigh);
    if isnan(locmin(m))
        disp(['iter #',num2str(m),' for all is not bimodal'])
    end
    locmin_agg(m) = calcLocMin2(N_agg,xloc_agg,mylow,myhigh);
    if isnan(locmin_agg(m))
        disp(['iter #',num2str(m),' for agg. is not bimodal'])
    end
    locmin_mrim(m) = calcLocMin2(N_mrim,xloc_mrim,mylow,myhigh);
    if isnan(locmin_mrim(m))
        disp(['iter #',num2str(m),' for mod. rimed is not bimodal'])
    end
    locmin_rim(m) = calcLocMin2(N_rim,xloc_rim,mylow,myhigh);
    if isnan(locmin_rim(m))
        disp(['iter #',num2str(m),' for graupel is not bimodal'])
    end
end

%
[slope,yint,yfit,yhi,ylo,gof] = calcLinFit(wind_array,locmin);
locmin_agg(locmin_agg<0.02) = nan;
locmin_mrim(locmin_mrim<0.02) = nan;
locmin_rim(locmin_rim<0.02) = nan;
[slope_agg,yint_agg,yfit_agg,yhi_agg,ylo_agg,gof_agg] = calcLinFit(wind_array,locmin_agg);
[slope_mrim,yint_mrim,yfit_mrim,yhi_mrim,ylo_mrim,gof_mrim] = calcLinFit(wind_array,locmin_mrim);
[slope_rim,yint_rim,yfit_rim,yhi_rim,ylo_rim,gof_rim] = calcLinFit(wind_array,locmin_rim);

%% plot v_cutoff(U_sfc)
ymin = 0.16;
ymax = 0.59;

lab_xpos = 0.02;
lab_ypos = 0.95;

h02 = figure;
set(h02,'units','inches','position',[0 0 pagewidth pageheight]);
h02a = tight_subplot(2,2,[0.01 0],[0.10 0.025],[0.08 0.02]);
axes(h02a(1))
    myconf = fill([wind_array,fliplr(wind_array)],[yhi',fliplr(ylo')],mybrown);
    myconf.FaceColor = [0.8 0.6 0.4];      
    myconf.EdgeColor = 'none';     
    hold on
    plot(wind_array,locmin,'o','Color',mybrown)
    grid on
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:2:10);
    xlabel('');
    set(gca,'XTickLabels','');
    ylabel('');
    text(1,0.5,['$v_{c}\mathbf{=',num2str(slope,'%4.3f'),'}U_{sfc}\mathbf{+',...
        num2str(yint,'%4.2f'),'}$'],'FontSize',text_size1,'Color',mybrown,'Interpreter','Latex');
    text(8,0.30,['$\mathbf{\tau = ',num2str(corr(wind_array',locmin','Type','Kendall',...
        'rows','complete'),'%4.2f'),'}$'],'FontSize',text_size1,...
        'Color',mybrown,'Interpreter','Latex');
    text(8,0.27,'$\mathbf{\rm{all}}$','FontSize',text_size1,...
        'Color',mybrown,'Interpreter','Latex');
    xlim([0 11])
    ylim([ymin ymax])
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(a)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    plot(linspace(0,12,100),linspace(0,12,100).*slope+yint,'Color',mybrown,...
        'LineWidth',2)
axes(h02a(2))
    myconf_agg = fill([wind_array(~isnan(locmin_agg)),fliplr(wind_array(~isnan(locmin_agg)))],...
        [yhi_agg',fliplr(ylo_agg')],'k');
    myconf_agg.FaceColor = [0.8 0.8 0.8];      
    myconf_agg.EdgeColor = 'none';     
    hold on
    plot(wind_array,locmin_agg,'ko')
    grid on
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:2:10);
    xlabel('');
    set(gca,'XTickLabels','');
    ylabel('');
    set(gca,'YTickLabels','');
    text(1,0.5,['$v_{c}\mathbf{=',num2str(slope_agg,'%4.3f'),'}U_{sfc}\mathbf{+',...
        num2str(yint_agg,'%4.2f'),'}$'],'FontSize',text_size1,'Interpreter','Latex');
    text(7.9,0.25,['$\mathbf{\tau = ',num2str(corr(wind_array',locmin_agg','Type',...
        'Kendall','rows','complete'),'%4.2f'),'}$'],'FontSize',text_size1,'Interpreter','Latex');
    text(7.9,0.22,'$\mathbf{\rm{aggregates}}$','FontSize',text_size1,...
        'Color','k','Interpreter','Latex');
    xlim([0 11])
    ylim([ymin ymax])
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(b)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    plot(linspace(0,12,100),linspace(0,12,100).*slope_agg+yint_agg,'Color','k',...
        'LineWidth',2)
axes(h02a(3))
    myconf_mrim = fill([wind_array(~isnan(locmin_mrim)),fliplr(wind_array(~isnan(locmin_mrim)))],...
        [yhi_mrim',fliplr(ylo_mrim')],'r');
    myconf_mrim.FaceColor = [1 0.8 0.8];      
    myconf_mrim.EdgeColor = 'none';     
    hold on
    plot(wind_array,locmin_mrim,'ro')
    grid on
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:2:10);
    xlabel('');
    y1 = ylabel('\textbf{Cutoff Fall Speed $(v_c)\mathbf{\,[m\,s^{-1}]}$}',...
        'FontSize',text_size1,'FontWeight','bold','Interpreter','Latex');
    y1.Position(2) = ymax;
    text(1,0.5,['$v_{c}\mathbf{=',num2str(slope_mrim,'%4.3f'),'}U_{sfc}\mathbf{+',...
        num2str(yint_mrim,'%4.2f'),'}$'],'FontSize',text_size1,'Color','r','Interpreter','Latex');
    text(7.3,0.30,['$\mathbf{\tau = ',num2str(corr(wind_array',locmin_mrim','Type',...
        'Kendall','rows','complete'),'%4.2f'),'}$'],'FontSize',text_size1,...
        'Color','r','Interpreter','Latex');
    text(7.3,0.27,'$\mathbf{\rm{moderately}}$','FontSize',text_size1,...
        'Color','r','Interpreter','Latex');
    text(7.3,0.24,'$\mathbf{\rm{rimed}}$','FontSize',text_size1,...
        'Color','r','Interpreter','Latex');
    xlim([0 11])
    ylim([ymin ymax])
    clipLabEnds('XTickLabels','low');
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(c)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    plot(linspace(0,12,100),linspace(0,12,100).*slope_mrim+yint_mrim,'Color','r',...
        'LineWidth',2)
axes(h02a(4))
    myconf_rim = fill([wind_array(~isnan(locmin_rim)),fliplr(wind_array(~isnan(locmin_rim)))],...
        [yhi_rim',fliplr(ylo_rim')],'r');
    myconf_rim.FaceColor = [0.8 0.8 1];      
    myconf_rim.EdgeColor = 'none';     
    hold on
    plot(wind_array,locmin_rim,'bo')
    grid on
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:2:10);
    x1 = xlabel('\textbf{Wind Speed $(U_{sfc})\mathbf{\,[m\,s^{-1}]}$}',...
        'FontSize',text_size1,'FontWeight','bold','Interpreter','Latex');
    x1.Position(1) = 0; x1.Position(2) = 0.11;
    ylabel('');
    set(gca,'YTickLabels','');
    text(1,0.5,['$v_{c}\mathbf{=',num2str(slope_rim,'%4.3f'),'}U_{sfc}\mathbf{+',...
        num2str(yint_rim,'%4.2f'),'}$'],'FontSize',text_size1,'Color','b','Interpreter','Latex');
    text(8,0.30,['$\mathbf{\tau = ',num2str(corr(wind_array',locmin_rim','Type',...
        'Kendall','rows','complete'),'%4.2f'),'}$'],'FontSize',text_size1,...
        'Color','b','Interpreter','Latex');
    text(8,0.27,'$\mathbf{\rm{rimed}}$','FontSize',text_size1,...
        'Color','b','Interpreter','Latex');
    xlim([0 11])
    ylim([ymin ymax])
    clipLabEnds('XTickLabels','low');
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(d)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    plot(linspace(0,12,100),linspace(0,12,100).*slope_rim+yint_rim,'Color','b',...
        'LineWidth',2)
    
orient(h02,'landscape');
print([fig_dir,myslash,'v_cutoff.pdf'],'-dpdf','-painters','-fillpage');

%% Compare MASC theta for different WF & wind speeds
ymax = 3e-2; % max value for y-axis (to standardize)
xmax = 90; % max value for x-axis (to standardize)

lab_xpos = 0.02;
lab_ypos = 0.90;
 
h2 = figure;
set(h2,'units','inches','position',[0 0 pagewidth pageheight]);
h2a = tight_subplot(2,3,0,[0.10 0.08],[0.12 0]);
axes(h2a(1))
    [N_theta1,xloc_theta1] = ksdensity(theta(inds1),mypts_theta);
    p1 = plot(xloc_theta1,N_theta1,'k-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    ylabel('\textbf{Unshielded}','FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex')
    title({['$\mathbf{>',num2str(wind_restr(1),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:15:90)
    set(gca,'XTickLabels','')
    clipLabEnds('YTickLabels','low');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(a)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    xq_theta1 = linspace(xloc_theta1(1),xloc_theta1(end),1e5);
    max_theta1 = xq_theta1(islocalmax(interp1(xloc_theta1,N_theta1,xq_theta1,...
        'spline'))==1)
axes(h2a(2))
    [N_theta2,xloc_theta2] = ksdensity(theta(inds2),mypts_theta);
    p2 = plot(xloc_theta2,N_theta2,'k-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    title({'$U_{sfc}\mathbf{\,[m\,s^{-1}]}$',['$\mathbf{\leq',num2str(wind_restr(1),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    yl2 = ylabel('\textbf{PDF $\mathbf{[^{\circ}]^{-1}}$}','FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex');
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:15:90)
    set(gca,'XTickLabels','')
    set(gca,'YTickLabels','')
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(b)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h2a(3))
    [N_theta3,xloc_theta3] = ksdensity(theta(inds3),mypts_theta);
    p3 = plot(xloc_theta3,N_theta3,'k-','LineWidth',2);    
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    title({['$\mathbf{\leq',num2str(wind_restr(2),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:15:90)
    set(gca,'XTickLabels','')
    set(gca,'YTickLabels','')
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(c)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    xq_theta3 = linspace(xloc_theta3(1),xloc_theta3(end),1e5);
    max_theta3 = xq_theta3(islocalmax(interp1(xloc_theta3,N_theta3,xq_theta3,...
        'spline'))==1)
axes(h2a(4))
    [N_theta4,xloc_theta4] = ksdensity(theta(inds4),mypts_theta);
    p4 = plot(xloc_theta4,N_theta4,'k-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    xlabel('')
    ylabel('\textbf{Shielded}','FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:15:90)
    clipLabEnds('XTickLabels','both');
    clipLabEnds('YTickLabels','low');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(d)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h2a(5))
    [N_theta5,xloc_theta5] = ksdensity(theta(inds5),mypts_theta);
    p5 = plot(xloc_theta5,N_theta5,'k-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    xlabel('\textbf{Orientation $(\theta)\mathbf{\,[^{\circ}]}$}',...
        'FontSize',text_size1,'FontWeight','bold','Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'YTickLabels','')
    set(gca,'XTick',0:15:90)
    clipLabEnds('XTickLabels','both');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(e)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h2a(6))
    [N_theta6,xloc_theta6] = ksdensity(theta(inds6),mypts_theta);
    p6 = plot(xloc_theta6,N_theta6,'k-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    grid on
    xlabel('')
    set(gca,'FontSize',text_size1)
    set(gca,'YTickLabels','')
    set(gca,'XTick',0:15:90)
    clipLabEnds('XTickLabels','both');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(f)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    xq_theta6 = linspace(xloc_theta6(1),xloc_theta6(end),1e5);
    max_theta6 = xq_theta6(islocalmax(interp1(xloc_theta6,N_theta6,xq_theta6,...
        'spline'))==1)
    
% center y-label
myyl2_x = yl2.Position(1);
yl2.Position(2) = 0;
yl2.Position(1) = myyl2_x - xmax - 30;
    
orient(h2,'landscape');
print([fig_dir,myslash,'theta_wf_usfc.pdf'],'-dpdf','-painters','-fillpage');    
    
%% Compare MASC D_max for different WF & wind speeds
ymin = 10^(-3.5);
ymax = 10^1; % max value for y-axis (to standardize)
xmax = 9.2; % max value for x-axis (to standardize)

lab_xpos = 0.05;
slope_thresh = 10^-3;

mypts_test = linspace(0,6,100);
 
h3 = figure;
set(h3,'units','inches','position',[0 0 pagewidth pageheight]);
h3a = tight_subplot(2,3,0,[0.10 0.08],[0.12 0]);
axes(h3a(1))
    [N_dmax1,xloc_dmax1] = ksdensity(dmax(inds1),mypts_test);
    test11 = N_dmax1 > slope_thresh;
    test12 = test11(2:end)-test11(1:end-1);
    test13 = find(test12~=0);
    N_dmax1new = N_dmax1(test13(1)+1:test13(2));
    xloc_dmax1new = xloc_dmax1(test13(1)+1:test13(2));
    p1a = plot(xloc_dmax1new,N_dmax1new,'k-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope1 = calcSlope(N_dmax1new,xloc_dmax1new,'r',2);
    text(2,10^(0),strcat('$\mathbf{\lambda=\,',num2str(round(slope1,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','r');
    xlim([0 xmax])
    ylim([ymin ymax])
    grid on
    ylabel('\textbf{Unshielded}','FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex')
    title({['$\mathbf{>',num2str(wind_restr(1),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTickLabels','')
    set(gca,'XTick',0:1:10)
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)*(myylim(2)/myylim(1))^lab_ypos,'\textbf{(a)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h3a(2))
    [N_dmax2,xloc_dmax2] = ksdensity(dmax(inds2),mypts_test);
    test21 = N_dmax2 > slope_thresh;
    test22 = test21(2:end)-test21(1:end-1);
    test23 = find(test22~=0);
    if length(test23)>1
        N_dmax2new = N_dmax2(test23(1)+1:test23(2));
        xloc_dmax2new = xloc_dmax2(test23(1)+1:test23(2));
    end
    if length(test23)>1
        p2a = plot(xloc_dmax2new,N_dmax2new,'k-','LineWidth',2);
    else
        p2a = plot(xloc_dmax2,N_dmax2,'k-','LineWidth',2);
    end
    set(gca,'YScale','log')
     hold on
    
    if length(test23)>1
        slope2 = calcSlope(N_dmax2new,xloc_dmax2new,'r',2);
    else
        slope2 = calcSlope(N_dmax2,xloc_dmax2,'r',2);
    end
    text(2,10^(0),strcat('$\mathbf{\lambda=\,',num2str(round(slope2,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','r');
    xlim([0 xmax])
    ylim([ymin ymax])
    grid on
    title({'$U_{sfc}\mathbf{\,[m\,s^{-1}]}$',['$\mathbf{\leq',num2str(wind_restr(1),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    yl2 = ylabel('\textbf{PDF $\mathbf{[mm]^{-1}}$}','FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex');
    set(gca,'FontSize',text_size1)
    set(gca,'XTickLabels','')
    set(gca,'YTickLabels','')
    set(gca,'XTick',0:1:10)
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)*(myylim(2)/myylim(1))^lab_ypos,'\textbf{(b)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h3a(3))
    [N_dmax3,xloc_dmax3] = ksdensity(dmax(inds3),mypts_test);
    test31 = N_dmax3 > slope_thresh;
    test32 = test31(2:end)-test31(1:end-1);
    test33 = find(test32~=0);
    N_dmax3new = N_dmax3(test33(1)+1:test33(2));
    xloc_dmax3new = xloc_dmax3(test33(1)+1:test33(2));
    p3a = plot(xloc_dmax3new,N_dmax3new,'k-','LineWidth',2);
    set(gca,'YScale','log')
     hold on
    slope3 = calcSlope(N_dmax3new,xloc_dmax3new,'r',2);
    text(2,10^(0),strcat('$\mathbf{\lambda=\,',num2str(round(slope3,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','r');
    xlim([0 xmax])
    ylim([ymin ymax])
    grid on
    title({['$\mathbf{\leq',num2str(wind_restr(2),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTickLabels','')
    set(gca,'YTickLabels','')
    set(gca,'XTick',0:1:10)
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)*(myylim(2)/myylim(1))^lab_ypos,'\textbf{(c)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h3a(4))
    [~,xloc_dmax4] = ksdensity(dmax(inds4));
    new_pts4 = linspace(xloc_dmax4(1),xloc_dmax4(end),2000);
    [N_dmax4,xloc_dmax4] = ksdensity(dmax(inds4),new_pts4);
    test41 = N_dmax4 > slope_thresh;
    test42 = test41(2:end)-test41(1:end-1);
    test43 = find(test42~=0);
    N_dmax4new = N_dmax4(test43(1)+1:test43(2));
    xloc_dmax4new = xloc_dmax4(test43(1)+1:test43(2));
    plot(xloc_dmax4new,N_dmax4new,'k-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope4 = calcSlope(N_dmax4new,xloc_dmax4new,'r',2);
    text(2,10^(0),strcat('$\mathbf{\lambda=\,',num2str(round(slope4,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','r');
    xlim([0 xmax])
    ylim([ymin ymax])
    grid on
    xlabel('')
    ylabel('\textbf{Shielded}','FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:1:10)
    clipLabEnds('XTickLabels','both');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)*(myylim(2)/myylim(1))^lab_ypos,'\textbf{(d)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h3a(5))
    [~,xloc_dmax5] = ksdensity(dmax(inds5));
    new_pts5 = linspace(xloc_dmax5(1),xloc_dmax5(end),2000);
    [N_dmax5,xloc_dmax5] = ksdensity(dmax(inds5),new_pts5);
    test51 = N_dmax5 > slope_thresh;
    test52 = test51(2:end)-test51(1:end-1);
    test53 = find(test52~=0);
    N_dmax5new = N_dmax5(test53(1)+1:test53(2));
    xloc_dmax5new = xloc_dmax5(test53(1)+1:test53(2));
    plot(xloc_dmax5new,N_dmax5new,'k-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope5 = calcSlope(N_dmax5new,xloc_dmax5new,'r',2);
    text(2,10^(0),strcat('$\mathbf{\lambda=\,',num2str(round(slope5,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','r');
    xlim([0 xmax])
    ylim([ymin ymax])
    grid on
    xlabel('\textbf{Maximum Dimension $(D_{max})\mathbf{\,[mm]}$}',...
        'FontSize',text_size1,'FontWeight','bold','Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'YTickLabels','')
    set(gca,'XTick',0:1:10)
    clipLabEnds('XTickLabels','both');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)*(myylim(2)/myylim(1))^lab_ypos,'\textbf{(e)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h3a(6))
    [~,xloc_dmax6] = ksdensity(dmax(inds6));
    new_pts6 = linspace(xloc_dmax6(1),xloc_dmax6(end),5000);
    [N_dmax6,xloc_dmax6] = ksdensity(dmax(inds6),new_pts6);
    test61 = N_dmax6 > slope_thresh;
    test62 = test61(2:end)-test61(1:end-1);
    test63 = find(test62~=0);
    N_dmax6new = N_dmax6(test63(1)+1:test63(2));
    xloc_dmax6new = xloc_dmax6(test63(1)+1:test63(2));
    plot(xloc_dmax6new,N_dmax6new,'k-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope6 = calcSlope(N_dmax6new,xloc_dmax6new,'r',2);
    text(2,10^(0),strcat('$\mathbf{\lambda=\,',num2str(round(slope6,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','r');
    xlim([0 xmax])
    ylim([ymin ymax])
    grid on
    xlabel('')
    set(gca,'FontSize',text_size1)
    set(gca,'YTickLabels','')
    set(gca,'XTick',0:1:10)
    clipLabEnds('XTickLabels','both');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)*(myylim(2)/myylim(1))^lab_ypos,'\textbf{(f)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    
% center y-label
myyl2_x = yl2.Position(1);
yl2.Position(2) = ymin;
yl2.Position(1) = myyl2_x - xmax - 3;
    
orient(h3,'landscape');
print([fig_dir,myslash,'dmax_wf_usfc.pdf'],'-dpdf','-painters','-fillpage');


%% MASC dmax (by type,#2)
ymin = 2e-3;
ymax = 10^1.3; % max value for y-axis (to standardize)
xmax = 8.2; % max value for x-axis (to standardize)

slope_thresh = 2e-3;
 
h7 = figure;
set(h7,'units','inches','position',[0 0 pagewidth 6]);
h7a = tight_subplot(1,3,0,[0.12 0.10],[0.12 0.02]);
axes(h7a(1))
        
    [~,xloc_dmax4a] = ksdensity(dmax(inds4a));
    new_pts4a = linspace(xloc_dmax4a(1),xloc_dmax4a(end),1000);
    [N_dmax4a,xloc_dmax4a] = ksdensity(dmax(inds4a),new_pts4a);
    test41a = N_dmax4a > slope_thresh;
    test42a = test41a(2:end)-test41a(1:end-1);
    test43a = find(test42a~=0);
    N_dmax4anew = N_dmax4a(test43a(1)+1:test43a(2));
    xloc_dmax4anew = xloc_dmax4a(test43a(1)+1:test43a(2));
    p4a = plot(xloc_dmax4anew,N_dmax4anew,'k-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope4a = calcSlope(N_dmax4anew,xloc_dmax4anew,'k',1.5);
    
    [~,xloc_dmax4b] = ksdensity(dmax(inds4b));
    new_pts4b = linspace(xloc_dmax4b(1),xloc_dmax4b(end),1000);
    [N_dmax4b,xloc_dmax4b] = ksdensity(dmax(inds4b),new_pts4b);
    test41b = N_dmax4b > slope_thresh;
    test42b = test41b(2:end)-test41b(1:end-1);
    test43b = find(test42b~=0);
    N_dmax4bnew = N_dmax4b(test43b(1)+1:test43b(2));
    xloc_dmax4bnew = xloc_dmax4b(test43b(1)+1:test43b(2));
    p4b = plot(xloc_dmax4bnew,N_dmax4bnew,'r-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope4b = calcSlope(N_dmax4bnew,xloc_dmax4bnew,'r',1.5);
    
    [~,xloc_dmax4c] = ksdensity(dmax(inds4c));
    new_pts4c = linspace(xloc_dmax4c(1),xloc_dmax4c(end),1000);
    [N_dmax4c,xloc_dmax4c] = ksdensity(dmax(inds4c),new_pts4c);
    test41c = N_dmax4c > slope_thresh;
    test42c = test41c(2:end)-test41c(1:end-1);
    test43c = find(test42c~=0);
    N_dmax4cnew = N_dmax4c(test43c(1)+1:test43c(2));
    xloc_dmax4cnew = xloc_dmax4c(test43c(1)+1:test43c(2));
    p4c = plot(xloc_dmax4cnew,N_dmax4cnew,'b-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope4c = calcSlope(N_dmax4cnew,xloc_dmax4cnew,'b',1.5);
    
    text(2,10^(1.10),strcat('$\mathbf{\lambda=\,',num2str(round(slope4a,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','k');
    text(2,10^(0.75),strcat('$\mathbf{\lambda=\,',num2str(round(slope4b,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','r');
    text(2,10^(0.40),strcat('$\mathbf{\lambda=\,',num2str(round(slope4c,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','b');
    xlim([0 xmax])
    ylim([ymin ymax])
    grid on
    xlabel('')
    ylabel('\textbf{PDF $\mathbf{[mm]^{-1}}$ (Shielded)}','FontSize',text_size1,...
        'FontWeight','bold','Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:1:10)
    clipLabEnds('XTickLabels','both');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)*(myylim(2)/myylim(1))^lab_ypos,'\textbf{(a)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    text(3,10^0,'\textbf{aggregates}','FontSize',...
        text_size1-2,'FontWeight','bold','Interpreter','Latex')
    text(3,10^-0.3,'\textbf{mod. rimed}','FontSize',text_size1-2,'Color','r',...
        'FontWeight','bold','Interpreter','Latex')
    text(3,10^-0.6,'\textbf{rimed}','FontSize',text_size1-2,'Color','b',...
        'FontWeight','bold','Interpreter','Latex')
    title({['$\mathbf{>',num2str(wind_restr(1),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
axes(h7a(2))
            
    [~,xloc_dmax5a] = ksdensity(dmax(inds5a));
    new_pts5a = linspace(xloc_dmax5a(1),xloc_dmax5a(end),1000);
    [N_dmax5a,xloc_dmax5a] = ksdensity(dmax(inds5a),new_pts5a);
    test51a = N_dmax5a > slope_thresh;
    test52a = test51a(2:end)-test51a(1:end-1);
    test53a = find(test52a~=0);
    N_dmax5anew = N_dmax5a(test53a(1)+1:test53a(2));
    xloc_dmax5anew = xloc_dmax5a(test53a(1)+1:test53a(2));
    p5a = plot(xloc_dmax5anew,N_dmax5anew,'k-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope5a = calcSlope(N_dmax5anew,xloc_dmax5anew,'k',1.5);
    
    [~,xloc_dmax5b] = ksdensity(dmax(inds5b));
    new_pts5b = linspace(xloc_dmax5b(1),xloc_dmax5b(end),1000);
    [N_dmax5b,xloc_dmax5b] = ksdensity(dmax(inds5b),new_pts5b);
    test51b = N_dmax5b > slope_thresh;
    test52b = test51b(2:end)-test51b(1:end-1);
    test53b = find(test52b~=0);
    N_dmax5bnew = N_dmax5b(test53b(1)+1:test53b(2));
    xloc_dmax5bnew = xloc_dmax5b(test53b(1)+1:test53b(2));
    p5b = plot(xloc_dmax5bnew,N_dmax5bnew,'r-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope5b = calcSlope(N_dmax5bnew,xloc_dmax5bnew,'r',1.5);
    
    [~,xloc_dmax5c] = ksdensity(dmax(inds5c));
    new_pts5c = linspace(xloc_dmax5c(1),xloc_dmax5c(end),1000);
    [N_dmax5c,xloc_dmax5c] = ksdensity(dmax(inds5c),new_pts5c);
    test51c = N_dmax5c > slope_thresh;
    test52c = test51c(2:end)-test51c(1:end-1);
    test53c = find(test52c~=0);
    N_dmax5cnew = N_dmax5c(test53c(1)+1:test53c(2));
    xloc_dmax5cnew = xloc_dmax5c(test53c(1)+1:test53c(2));
    p5c = plot(xloc_dmax5cnew,N_dmax5cnew,'b-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope5c = calcSlope(N_dmax5cnew,xloc_dmax5cnew,'b',1.5);
    
    text(2,10^(1.10),strcat('$\mathbf{\lambda=\,',num2str(round(slope5a,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','k');
    text(2,10^(0.75),strcat('$\mathbf{\lambda=\,',num2str(round(slope5b,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','r');
    text(2,10^(0.40),strcat('$\mathbf{\lambda=\,',num2str(round(slope5c,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','b');
    xlim([0 xmax])
    ylim([ymin ymax])
    grid on
    xlabel('\textbf{Maximum Dimension $(D_{max})\mathbf{\,[mm]}$}',...
        'FontSize',text_size1,'FontWeight','bold','Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'YTickLabels','')
    set(gca,'XTick',0:1:10)
    clipLabEnds('XTickLabels','both');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)*(myylim(2)/myylim(1))^lab_ypos,'\textbf{(b)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    title({'$U_{sfc}\mathbf{\,[m\,s^{-1}]}$',['$\mathbf{\leq',num2str(wind_restr(1),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
axes(h7a(3))
            
    [~,xloc_dmax6a] = ksdensity(dmax(inds6a));
    new_pts6a = linspace(xloc_dmax6a(1),xloc_dmax6a(end),500);
    [N_dmax6a,xloc_dmax6a] = ksdensity(dmax(inds6a),new_pts6a);
    test61a = N_dmax6a > slope_thresh;
    test62a = test61a(2:end)-test61a(1:end-1);
    test63a = find(test62a~=0);
    N_dmax6anew = N_dmax6a(test63a(1)+1:test63a(2));
    xloc_dmax6anew = xloc_dmax6a(test63a(1)+1:test63a(2));
    plot(xloc_dmax6anew,N_dmax6anew,'k-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope6a = calcSlope(N_dmax6anew,xloc_dmax6anew,'k',1.5);
    
    [N_dmax6b,xloc_dmax6b] = ksdensity(dmax(inds6b));
    test61b = N_dmax6b > slope_thresh;
    test62b = test61b(2:end)-test61b(1:end-1);
    test63b = find(test62b~=0);
    N_dmax6bnew = N_dmax6b(test63b(1)+1:test63b(2));
    xloc_dmax6bnew = xloc_dmax6b(test63b(1)+1:test63b(2));
    plot(xloc_dmax6bnew,N_dmax6bnew,'r-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope6b = calcSlope(N_dmax6bnew,xloc_dmax6bnew,'r',1.5);
    
    [N_dmax6c,xloc_dmax6c] = ksdensity(dmax(inds6c));
    test61c = N_dmax6c > slope_thresh;
    test62c = test61c(2:end)-test61c(1:end-1);
    test63c = find(test62c~=0);
    N_dmax6cnew = N_dmax6c(test63c(1)+1:test63c(2));
    xloc_dmax6cnew = xloc_dmax6c(test63c(1)+1:test63c(2));
    plot(xloc_dmax6cnew,N_dmax6cnew,'b-','LineWidth',2);
    set(gca,'YScale','log')
    hold on
    slope6c = calcSlope(N_dmax6cnew,xloc_dmax6cnew,'b',1.5);
    
    text(2,10^(1.10),strcat('$\mathbf{\lambda=\,',num2str(round(slope6a,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','k');
    text(2,10^(0.75),strcat('$\mathbf{\lambda=\,',num2str(round(slope6b,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','r');
    text(2,10^(0.40),strcat('$\mathbf{\lambda=\,',num2str(round(slope6c,2)),...
        '\,\rm{mm^{-1}}}$'),'Interpreter','Latex','FontSize',text_size1+2,'Color','b');
    xlim([0 xmax])
    ylim([ymin ymax])
    grid on
    xlabel('')
    set(gca,'FontSize',text_size1)
    set(gca,'YTickLabels','')
    set(gca,'XTick',0:1:10)
    clipLabEnds('XTickLabels','both');
    myxlim = xlim;
    myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)*(myylim(2)/myylim(1))^lab_ypos,'\textbf{(c)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    title({['$\mathbf{\leq',num2str(wind_restr(2),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
  
set(h7,'papersize',[6 11])            
orient(h7,'landscape');
print([fig_dir,myslash,'dmax_chi2.pdf'],'-dpdf','-painters','-fillpage');

%% MASC v (by type,#2)
ymax = 3; % max value for y-axis (to standardize)
xmax = 1.8; % max value for x-axis (to standardize)
 
h7v = figure;
set(h7v,'units','inches','position',[0 0 pagewidth 6]);
h7va = tight_subplot(1,3,0,[0.12 0.10],[0.12 0.02]);
axes(h7va(1))        
    [N_v6a,xloc_v6a] = ksdensity(v(inds6a),mypts_fs);
    plot(xloc_v6a,N_v6a,'k-','LineWidth',2);
    hold on
    [N_v6b,xloc_v6b] = ksdensity(v(inds6b),mypts_fs);
    plot(xloc_v6b,N_v6b,'r-','LineWidth',2);
    [N_v6c,xloc_v6c] = ksdensity(v(inds6c),mypts_fs);
    plot(xloc_v6c,N_v6c,'b-','LineWidth',2);
    max_v6a = isMyMax(xloc_v6a,N_v6a,0.5,1.5)
    max_v6b = isMyMax(xloc_v6b,N_v6b,0.5,1.5)
    max_v6c = isMyMax(xloc_v6c,N_v6c,0.5,1.5)
    xlim([0 xmax])
    ylim([0 ymax])
    ylabel('\textbf{PDF $\mathbf{[m\,s^{-1}]^{-1}}$ (Shielded)}','FontSize',...
        text_size1,'FontWeight','bold','Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:0.5:xmax)
    grid on;
    clipLabEnds('XTickLabels','low');
    text(0.25,2.50,'\textbf{aggregates}','FontSize',...
        text_size1-2,'FontWeight','bold','Interpreter','Latex')
    text(0.25,2.25,'\textbf{mod. rimed}','FontSize',text_size1-2,'Color','r',...
        'FontWeight','bold','Interpreter','Latex')
    text(0.25,2.00,'\textbf{rimed}','FontSize',text_size1-2,'Color','b',...
        'FontWeight','bold','Interpreter','Latex')
    title({['$\mathbf{\leq',num2str(wind_restr(2),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    myxlim = xlim; myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(a)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h7va(2))        
    [N_v7a,xloc_v7a] = ksdensity(v(inds7a),mypts_fs);
    plot(xloc_v7a,N_v7a,'k-','LineWidth',2);
    hold on
    [N_v7b,xloc_v7b] = ksdensity(v(inds7b),mypts_fs);
    plot(xloc_v7b,N_v7b,'r-','LineWidth',2);
    [N_v7c,xloc_v7c] = ksdensity(v(inds7c),mypts_fs);
    plot(xloc_v7c,N_v7c,'b-','LineWidth',2);
    max_v7a = isMyMax(xloc_v7a,N_v7a,0.5,1.5)
    max_v7b = isMyMax(xloc_v7b,N_v7b,0.5,1.5)
    max_v7c = isMyMax(xloc_v7c,N_v7c,0.5,1.5)
    xlim([0 xmax])
    ylim([0 ymax])
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:0.5:xmax)
    grid on;
    clipLabEnds('XTickLabels','low');
    set(gca,'YTickLabels','');
    xlabel('\textbf{Fall Speed $(v_p)\mathbf{\,[m\,s^{-1}]}$}',...
        'FontSize',text_size1,'FontWeight','bold','Interpreter','Latex')
    title({'$U_{sfc}\mathbf{\,[m\,s^{-1}]}$',['$\mathbf{\leq',num2str(wind_restr2(1),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(b)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h7va(3))        
    [N_v8a,xloc_v8a] = ksdensity(v(inds8a),mypts_fs);
    plot(xloc_v8a,N_v8a,'k-','LineWidth',2);
    hold on
    [N_v8b,xloc_v8b] = ksdensity(v(inds8b),mypts_fs);
    plot(xloc_v8b,N_v8b,'r-','LineWidth',2);
    [N_v8c,xloc_v8c] = ksdensity(v(inds8c),mypts_fs);
    plot(xloc_v8c,N_v8c,'b-','LineWidth',2);
    max_v8a = isMyMax(xloc_v8a,N_v8a,0.5,1.5)
    max_v8b = isMyMax(xloc_v8b,N_v8b,0.5,1.5)
    max_v8c = isMyMax(xloc_v8c,N_v8c,0.5,1.5)
    xlim([0 xmax])
    ylim([0 ymax])
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:0.5:xmax)
    grid on;
    clipLabEnds('XTickLabels','low');
    set(gca,'YTickLabels','');
    title({['$\mathbf{\leq',num2str(wind_restr2(2),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(c)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    
set(h7v,'papersize',[6 11])            
orient(h7v,'landscape');
print([fig_dir,myslash,'v_chi2.pdf'],'-dpdf','-painters','-fillpage');


%% MASC theta (by type,#2)
ymax = 0.035; % max value for y-axis (to standardize)
xmax = 90; % max value for x-axis (to standardize)
lab_xpos = 0.80;
lab_ypos = 0.90;
 
h7t = figure;
set(h7t,'units','inches','position',[0 0 pagewidth 6]);
h7ta = tight_subplot(1,3,0,[0.12 0.10],[0.12 0.02]);
axes(h7ta(1))        
    [N_t6a,xloc_t6a] = ksdensity(theta(inds6a),mypts_theta);
    p6a = plot(xloc_t6a,N_t6a,'k-','LineWidth',2);
    hold on
    [N_t6b,xloc_t6b] = ksdensity(theta(inds6b),mypts_theta);
    p6b = plot(xloc_t6b,N_t6b,'r-','LineWidth',2);
    [N_t6c,xloc_t6c] = ksdensity(theta(inds6c),mypts_theta);
    p6c= plot(xloc_t6c,N_t6c,'b-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    ylabel('\textbf{PDF $\mathbf{\,[^{\circ}]^{-1}}$ (Shielded)}','FontSize',...
        text_size1,'FontWeight','bold','Interpreter','Latex')
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:15:xmax)
    grid on;
    clipLabEnds('XTickLabels','both');
    text(30,0.028,'\textbf{aggregates}','FontSize',...
        text_size1-2,'FontWeight','bold','Interpreter','Latex')
    text(30,0.025,'\textbf{mod. rimed}','FontSize',text_size1-2,'Color','r',...
        'FontWeight','bold','Interpreter','Latex')
    text(30,0.022,'\textbf{rimed}','FontSize',text_size1-2,'Color','b',...
        'FontWeight','bold','Interpreter','Latex')
    title({['$\mathbf{\leq',num2str(wind_restr(2),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    myxlim = xlim; myylim = ylim;
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(a)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h7ta(2))        
    [N_t7a,xloc_t7a] = ksdensity(theta(inds7a),mypts_theta);
    p7a = plot(xloc_t7a,N_t7a,'k-','LineWidth',2);
    hold on
    [N_t7b,xloc_t7b] = ksdensity(theta(inds7b),mypts_theta);
    p7b = plot(xloc_t7b,N_t7b,'r-','LineWidth',2);
    [N_t7c,xloc_t7c] = ksdensity(theta(inds7c),mypts_theta);
    p7c= plot(xloc_t7c,N_t7c,'b-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:15:xmax)
    grid on;
    clipLabEnds('XTickLabels','both');
    set(gca,'YTickLabels','');
    xlabel('\textbf{Orientation $(\theta)\mathbf{\,[^{\circ}]}$}',...
        'FontSize',text_size1,'FontWeight','bold','Interpreter','Latex')
    title({'$U_{sfc}\mathbf{\,[m\,s^{-1}]}$',['$\mathbf{\leq',num2str(wind_restr2(1),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(b)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
axes(h7ta(3))        
    [N_t8a,xloc_t8a] = ksdensity(theta(inds8a),mypts_theta);
    p8a = plot(xloc_t8a,N_t8a,'k-','LineWidth',2);
    hold on
    [N_t8b,xloc_t8b] = ksdensity(theta(inds8b),mypts_theta);
    p8b = plot(xloc_t8b,N_t8b,'r-','LineWidth',2);
    [N_t8c,xloc_t8c] = ksdensity(theta(inds8c),mypts_theta);
    p8c= plot(xloc_t8c,N_t8c,'b-','LineWidth',2);
    xlim([0 xmax])
    ylim([0 ymax])
    set(gca,'FontSize',text_size1)
    set(gca,'XTick',0:15:xmax)
    grid on;
    clipLabEnds('XTickLabels','both');
    set(gca,'YTickLabels','');
    title({['$\mathbf{\leq',num2str(wind_restr2(2),'%4.1f'),...
        '}$']},'FontSize',text_size1,'Interpreter','Latex')
    text(myxlim(1)+lab_xpos*(myxlim(2)-myxlim(1)),...
        myylim(1)+lab_ypos*(myylim(2)-myylim(1)),'\textbf{(c)}','FontSize',...
        text_size1+2,'FontWeight','bold','Interpreter','Latex')
    
set(h7t,'papersize',[6 11])            
orient(h7t,'landscape');
print([fig_dir,myslash,'theta_chi2.pdf'],'-dpdf','-painters','-fillpage');

%% Reynolds numbers
Re2 = calcRe2(dmax,dmax.*alpha,theta,v,p_sfc.*1e3,T_sfc+273.15);
Re_filt = Re2(intersect(inds6,find(v>v6_locmin)));
my_incr2 = 0.05/2;
myquant = quantile(Re_filt,my_incr2:my_incr2:1-my_incr2);
myquant(1)
myquant(end)