function winds_scatter_wc_compare_regTA( flag_wc16_avg_exact , flag_v2_std_or_mc , height_examined_m ,  WC_nums , avg_pd_min , marker_area , ws_low_lim_m_s , ws_upp_lim_m_s )
%winds_scatter Scatter plots of wind measurements from Wind Cube LiDARs
%    Scatter plots of wind measurements from Wind Cube LiDARs.

% Author: Paul T. Quelet
%
% Last Update Noted: 26 May 2015
%
% Made with Matlab 2014b
% IMPORTANT NOTE: GRAPHICS WILL NOT WORK WITH EARLIER VERSIONS OF CODE
%
% Example function call:
% winds_scatter_wc_compare_regTA( 'exact' , 'mc' , 50 , 2 : 4 , 2 , 15.0 , 0.0 , 18.0 )
%
% =========================================================================
% INPUT                 DESCRIPTION                 TYPE        UNITS
% -------------------------------------------------------------------------
% flag_wc16_avg_exact   flag for WC16 avg or exact  string      ['avg' or 'exact']
% flag_v2_std_or_mc 'std' for standard files        string
%                   'mc'  for motion compensated files (normal)
% height_examined_m     flag for height(s) examined double      [m] 50, 100, 150
% WC_nums               Wind Cube Numbers           double      [unitless] 1 is 49, 2 is 61, 3 is 68, 4 is 16
% avg_pd_min            Avg. Period of WC STA files integer     [minutes]
% marker_area           area scatter marker         double      [unitless]
% ws_low_lim_m_s        Lower Limit Wind Speed      double      [m/sec]
% ws_upp_lim_m_s        Upper Limit Wind Speed      double      [m/sec]
% =========================================================================

%% Set Parameters

% Note font parameters:

fss                                             = 'FontSize';

font_size                                       = 14 ;

% Note the relevant heights for examination:

height_examined_m_BAO_lower                     = 50 ;

height_examined_m_BAO_middle                    = 100 ;

height_examined_m_BAO_upper                     = 150 ;

% Note the current experiement:

curr_exper_abbr                                 = 'XPIA' ;

% Other experiments string abbreviations include the following:
% 'BLADE' , 'NWTC' , 'CWEX11' , 'TWICS' , 'TODS' , 'CWEX13' , 'LUMEX'

% Date string format out:
% See http://www.mathworks.com/help/matlab/ref/datestr.html#bteqx6q for more formatting:

datestr_formatOut                               = 'dd-mmm-yyyy HH:MM:SS' ;

%% Load relvant data from WindCubes

% Set in directory for WindCube data:

inDir_WC                                        = 'E:\XPIA_Research\XPIA_Dataviewer_Matlab_Plot_Codes\WLS7\WC_regTA_output\' ;

% Loop over the WCs and load the data:

for j_wc = WC_nums
    
    % Assign information base upon Wind Cube Number:
    
    switch j_wc
        case 1
            WCnum                               = 49 ;
            WCnum_str                           = num2str( WCnum ) ;
        case 2
            WCnum                               = 61 ;
            WCnum_str                           = num2str( WCnum ) ;
        case 3
            WCnum                               = 68 ;
            WCnum_str                           = num2str( WCnum ) ;
        case 4
            WCnum                               = 16 ;
            WCnum_str                           = num2str( WCnum ) ;
    end
    
    load( strcat( inDir_WC , curr_exper_abbr , '_' , WCnum_str, '_' , num2str( avg_pd_min ) , 'min_' , flag_v2_std_or_mc , '_regTA.mat' ) )
    
end % (j_wc) End loop over WindCube Numbers

clearvars -except regTA alt* ws* wdir* fss font_size text_fs font_size curr_exper_abbr ...
    height_examined_m* avg_pd_min flag_wc16_avg_exact ...
    datestr_formatOut marker_area

%% Extract the dates relevant for the case study:

% Set the times in between to examine:

lower_datetime                                  = datetime( 2015, 03 , 02 , 00 , 00 , 00 ) ;
upper_datetime                                  = datetime( 2015, 06 , 01 , 00 , 00 , 00 ) ;

% Use the new "isbetween" function to get new true/false vectors to plug into the wind speed arrays:

specified_time_inds                             = isbetween( regTA , lower_datetime , upper_datetime ) ;

% Get the smaller WindCube matrices for the specified dates:

ws_case_61                                      = ws_regTA_61( specified_time_inds , : ) ;
ws_case_68                                      = ws_regTA_68( specified_time_inds , : ) ;
ws_case_16                                      = ws_regTA_16( specified_time_inds , : ) ;

%% Extract only the relevant WC Heights:

ws_WC61_40m_m_s                                 = ws_case_61( : , alt_61 == height_examined_m_BAO_lower - 10 ) ; % Index 1 is 40m
ws_WC68_40m_m_s                                 = ws_case_68( : , alt_68 == height_examined_m_BAO_lower - 10 ) ; % Index 1 is 40m
ws_WC16_40m_m_s                                 = ws_case_16( : , alt_16 == height_examined_m_BAO_lower - 10 ) ; % Index 1 is 40m

ws_WC61_60m_m_s                                 = ws_case_61( : , alt_61 == height_examined_m_BAO_lower + 10 ) ; % Index 2 is 60m
ws_WC68_60m_m_s                                 = ws_case_68( : , alt_68 == height_examined_m_BAO_lower + 10 ) ; % Index 2 is 60m
ws_WC16_60m_m_s                                 = ws_case_16( : , alt_16 == height_examined_m_BAO_lower + 10 ) ; % Index 3 is 60m

ws_WC16_50m_m_s                                 = ws_case_16( : , alt_16 == height_examined_m_BAO_lower ) ; % Index 2 is 50m, WC16 only

ws_WC61_100m_m_s                                = ws_case_61( : , alt_61 == height_examined_m_BAO_middle ) ; % Index 4 is 100m
ws_WC68_100m_m_s                                = ws_case_68( : , alt_68 == height_examined_m_BAO_middle ) ; % Index 4 is 100m
ws_WC16_100m_m_s                                = ws_case_16( : , alt_16 == height_examined_m_BAO_middle ) ; % Index 5 is 100m

ws_WC61_140m_m_s                                = ws_case_61( : , alt_61 == height_examined_m_BAO_upper - 10 ) ; % Index 6 is 140m
ws_WC68_140m_m_s                                = ws_case_68( : , alt_68 == height_examined_m_BAO_upper - 10 ) ; % Index 6 is 140m
ws_WC16_140m_m_s                                = ws_case_16( : , alt_16 == height_examined_m_BAO_upper - 10 ) ; % Index 6 is 140m

ws_WC61_160m_m_s                                = ws_case_61( : , alt_61 == height_examined_m_BAO_upper + 10 ) ; % Index 7 is 160m
ws_WC68_160m_m_s                                = ws_case_68( : , alt_68 == height_examined_m_BAO_upper + 10 ) ; % Index 7 is 160m
ws_WC16_160m_m_s                                = ws_case_16( : , alt_16 == height_examined_m_BAO_upper + 10 ) ; % Index 9 is 160m

ws_WC16_150m_m_s                                = ws_case_16( : , alt_16 == height_examined_m_BAO_upper ) ; % Index 2 is 150m, WC16 only

%% Horizontally concatenate WCs, then take mean (being restrictive) of heights above & below height of interest:

ws_WC61_40m60m_matrix_m_s                       = [ ws_WC61_40m_m_s ws_WC61_60m_m_s ] ;
ws_WC61_50m_m_s                                 = mean( ws_WC61_40m60m_matrix_m_s , 2 ) ;

ws_WC68_40m60m_matrix_m_s                       = [ ws_WC68_40m_m_s ws_WC68_60m_m_s ] ;
ws_WC68_50m_m_s                                 = mean( ws_WC68_40m60m_matrix_m_s , 2 ) ;

switch flag_wc16_avg_exact
    case 'avg'
        ws_WC16_40m60m_matrix_m_s               = [ ws_WC16_40m_m_s ws_WC16_60m_m_s ] ;
        ws_WC16_50m_m_s                         = mean( ws_WC16_40m60m_matrix_m_s , 2 ) ;
    case 'exact'
        % Use as is!
end

ws_WC61_140m160m_matrix_m_s                     = [ ws_WC61_140m_m_s ws_WC61_160m_m_s ] ;
ws_WC61_150m_m_s                                = mean( ws_WC61_140m160m_matrix_m_s , 2 ) ;

ws_WC68_140m160m_matrix_m_s                     = [ ws_WC68_140m_m_s ws_WC68_160m_m_s ] ;
ws_WC68_150m_m_s                                = mean( ws_WC68_140m160m_matrix_m_s , 2 ) ;

switch flag_wc16_avg_exact
    case 'avg'
        ws_WC16_140m160m_matrix_m_s             = [ ws_WC16_140m_m_s ws_WC16_160m_m_s ] ;
        ws_WC16_150m_m_s                        = mean( ws_WC16_140m160m_matrix_m_s , 2 ) ;
    case 'exact'
        % Use as is!
end

%% Cut off first values of WC time series to make same length:

ws_WC61_50m_m_s                                 = ws_WC61_50m_m_s(  2 : end ) ;
ws_WC68_50m_m_s                                 = ws_WC68_50m_m_s(  2 : end ) ;
ws_WC16_50m_m_s                                 = ws_WC16_50m_m_s(  2 : end ) ;

ws_WC61_100m_m_s                                = ws_WC61_100m_m_s( 2 : end ) ;
ws_WC68_100m_m_s                                = ws_WC68_100m_m_s( 2 : end ) ;
ws_WC16_100m_m_s                                = ws_WC16_100m_m_s( 2 : end ) ;

ws_WC61_150m_m_s                                = ws_WC61_150m_m_s( 2 : end ) ;
ws_WC68_150m_m_s                                = ws_WC68_150m_m_s( 2 : end ) ;
ws_WC16_150m_m_s                                = ws_WC16_150m_m_s( 2 : end ) ;

%% Create figure properties

% Create figure:
fig_name                                        = strcat( 'scatter_plot_WCs' ) ;
figure1                                         = figure( 'Name' , fig_name, 'units' , 'normalized' , 'outerposition', [ 0 0 1 1 ] );

% Create axes:
axes1                                           = axes( 'Parent' , figure1 ) ;

%% X-limits of the axes:

% Let the limits of the axes be the same so that scatterplot a square

xlim( axes1 , [ ws_low_lim_m_s ws_upp_lim_m_s ] ) ;

%% Y-limits of the axes:

ylim( axes1 , [ws_low_lim_m_s ws_upp_lim_m_s ] );

hold( axes1 , 'on' ) ;

%% Set Additional Plot properties

% Create xlabel:
x_label_name                                    = strcat( 'WC16 v2 (Leosphere) WS [m s^{-1}]' ) ;
xlabel( x_label_name , fss, font_size ) ;

% Create ylabel:
y_label_name                                    = strcat( 'WindCube v1 LiDARs WS [m s^{-1}]' ) ;
ylabel( y_label_name , fss , font_size ) ;

% Create title:
title_fontsize                                  = 17.5 ;
title_top_str                                   = ...
    [ 'WC16 v2 WS vs. WC61 & WC68 v1 WS Scatterplot, ' , num2str( height_examined_m ) , 'm' ]  ;
dates_str                               = [ datestr( lower_datetime , datestr_formatOut ) , ' to ' , datestr( upper_datetime , datestr_formatOut ) ] ;
title( { title_top_str ; dates_str } , fss , title_fontsize ) ;

%% Make Scatterplot

switch height_examined_m
    
    case 50
        % Hold the figure to continue plotting on it:
        hold on
        first_marker_color_rgb                  = 'r' ;
        scatter( ws_WC16_50m_m_s , ws_WC61_50m_m_s , marker_area , first_marker_color_rgb  , 'Marker' , '.' ) ;
                
        % Hold the figure to continue plotting on it:
        hold on
        second_marker_color_rgb                 = 'b' ;
        scatter( ws_WC16_50m_m_s , ws_WC68_50m_m_s , marker_area , second_marker_color_rgb  , 'Marker' , '.' ) ;
        
    case 100
        % Hold the figure to continue plotting on it:
        hold on
        first_marker_color_rgb                  = 'r' ;
        scatter( ws_WC16_100m_m_s , ws_WC61_100m_m_s , marker_area , first_marker_color_rgb  , 'Marker' , '.' ) ;

                
        % Hold the figure to continue plotting on it:
        hold on
        second_marker_color_rgb                 = 'b' ;
        scatter( ws_WC16_100m_m_s , ws_WC68_100m_m_s , marker_area , second_marker_color_rgb  , 'Marker' , '.' ) ;
        
    case 150
        % Hold the figure to continue plotting on it:
        hold on
        first_marker_color_rgb                  = 'r' ;
        scatter( ws_WC16_150m_m_s , ws_WC61_150m_m_s , marker_area , first_marker_color_rgb  , 'Marker' , '.' ) ;
        
        % Hold the figure to continue plotting on it:
        hold on
        second_marker_color_rgb                 = 'b' ;
        scatter( ws_WC16_150m_m_s , ws_WC68_150m_m_s , marker_area , second_marker_color_rgb  , 'Marker' , '.' ) ;
        
end % End switch statement over heights examined

%% Plot One-to-One Line

% Hold the figure to continue plotting on it:

hold on

% Set the 1 to 1 line properties:

m_slope                                         = 1.0 ;
b_y_int                                         = 0.0 ;
h_one_to_one_line                               = refline( m_slope, b_y_int ) ;
h_one_to_one_line.Color                         = 'k' ; % Make the reference line black

%% Calculate slope and R^2 of best fit lines

% Set the order of polynomial to fit to:

poly_order                                      = 1.0 ;

switch height_examined_m
    
    case 50
        
        % Setup good (non-nan) value indicies to fit to the values:
        
        non_nan_inds_61                         = intersect( find( ~isnan( ws_WC16_50m_m_s ) ) , find( ~isnan( ws_WC61_50m_m_s ) ) ) ;
        non_nan_inds_68                         = intersect( find( ~isnan( ws_WC16_50m_m_s ) ) , find( ~isnan( ws_WC68_50m_m_s ) ) ) ;
        
        % Calculate the best fit line with a first order polynomial fit:
        
        lin_coeffs_61                           = polyfit( ws_WC16_50m_m_s( non_nan_inds_61 ) , ws_WC61_50m_m_s( non_nan_inds_61 ) , poly_order ) ;
        lin_coeffs_68                           = polyfit( ws_WC16_50m_m_s( non_nan_inds_68 ) , ws_WC68_50m_m_s( non_nan_inds_68 ) , poly_order ) ;
                
    case 100
        
        % Setup good (non-nan) value indicies to fit to the values:
        
        non_nan_inds_61                         = intersect( find( ~isnan( ws_WC16_100m_m_s ) ) , find( ~isnan( ws_WC61_100m_m_s ) ) ) ;
        non_nan_inds_68                         = intersect( find( ~isnan( ws_WC16_100m_m_s ) ) , find( ~isnan( ws_WC68_100m_m_s ) ) ) ;
        
        % Calculate the best fit line with a first order polynomial fit:
        
        lin_coeffs_61                           = polyfit( ws_WC16_100m_m_s( non_nan_inds_61 ) , ws_WC61_100m_m_s( non_nan_inds_61 ) , poly_order ) ;
        lin_coeffs_68                           = polyfit( ws_WC16_100m_m_s( non_nan_inds_68 ) , ws_WC68_100m_m_s( non_nan_inds_68 ) , poly_order ) ;
        
    case 150
        
        % Setup good (non-nan) value indicies to fit to the values:
        
        non_nan_inds_61                         = intersect( find( ~isnan( ws_WC16_150m_m_s ) ) , find( ~isnan( ws_WC61_150m_m_s ) ) ) ;
        non_nan_inds_68                         = intersect( find( ~isnan( ws_WC16_150m_m_s ) ) , find( ~isnan( ws_WC68_150m_m_s ) ) ) ;
        
        % Calculate the best fit line with a first order polynomial fit:
        
        lin_coeffs_61                           = polyfit( ws_WC16_150m_m_s( non_nan_inds_61 ) , ws_WC61_150m_m_s( non_nan_inds_61 ) , poly_order ) ;
        lin_coeffs_68                           = polyfit( ws_WC16_150m_m_s( non_nan_inds_68 ) , ws_WC68_150m_m_s( non_nan_inds_68 ) , poly_order ) ;
   
end % End switch statement over heights examined

% Get the slope and y-intercepts of the lines:

m_slope_61                                      = lin_coeffs_61( 1 ) ;
m_slope_68                                      = lin_coeffs_68( 1 ) ;

b_int_61                                        = lin_coeffs_61( 2 ) ;
b_int_68                                        = lin_coeffs_68( 2 ) ;

switch height_examined_m
    
    case 50
        
        % Calculate the correlation coefficient and then square it:
        
        r_sq_61                                 = corr( ws_WC16_50m_m_s( non_nan_inds_61 ) , ws_WC61_50m_m_s( non_nan_inds_61 ) ) .^ 2.0 ;
        r_sq_68                                 = corr( ws_WC16_50m_m_s( non_nan_inds_68 ) , ws_WC68_50m_m_s( non_nan_inds_68 ) ) .^ 2.0 ;
        
        % Also store the number of points:
        
        num_pts_61                              = length( ws_WC61_50m_m_s( non_nan_inds_61 ) ) ;
        num_pts_68                              = length( ws_WC68_50m_m_s( non_nan_inds_68 ) ) ;
        
    case 100
        
        % Calculate the correlation coefficient and then square it:
        
        r_sq_61                                 = corr( ws_WC16_100m_m_s( non_nan_inds_61 ) , ws_WC61_100m_m_s( non_nan_inds_61 ) ) .^ 2.0 ;
        r_sq_68                                 = corr( ws_WC16_100m_m_s( non_nan_inds_68 ) , ws_WC68_100m_m_s( non_nan_inds_68 ) ) .^ 2.0 ;
        
        % Also store the number of points:
        
        num_pts_61                              = length( ws_WC61_100m_m_s( non_nan_inds_61 ) ) ;
        num_pts_68                              = length( ws_WC68_100m_m_s( non_nan_inds_68 ) ) ;
        
    case 150
        
        % Calculate the correlation coefficient and then square it:
        
        r_sq_61                                 = corr( ws_WC16_150m_m_s( non_nan_inds_61 ) , ws_WC61_150m_m_s( non_nan_inds_61 ) ) .^ 2.0 ;
        r_sq_68                                 = corr( ws_WC16_150m_m_s( non_nan_inds_68 ) , ws_WC68_150m_m_s( non_nan_inds_68 ) ) .^ 2.0 ;
        
        % Also store the number of points:
        
        num_pts_61                              = length( ws_WC61_150m_m_s( non_nan_inds_61 ) ) ;
        num_pts_68                              = length( ws_WC68_150m_m_s( non_nan_inds_68 ) ) ;
        
end % End switch statement over height examined

%% Display Coefficients on the chart

right_margin_displace                           = 6.00 ;

OT                                              = 1000.0 ; % One Thousand, used to round decimals to the thousandths place

eq_str                                          = 'y = mx + b' ;
text( ws_upp_lim_m_s - right_margin_displace , ws_low_lim_m_s + 4.50 , eq_str       , 'Color' , 'black' ) ;

r_sq_str_61                                     = ['WC61, R^2 = ' , num2str( round( r_sq_61 * OT ) / OT ), ', N = ' , num2str( num_pts_61 ) ] ;
text( ws_upp_lim_m_s - right_margin_displace , ws_low_lim_m_s + 3.50 , r_sq_str_61 , 'Color' , 'red') ;
lin_eq_str_61                                   = [ '(NCAR), m = ' , num2str( round( m_slope_61 * OT ) / OT ) , ', b = ' , num2str( round( b_int_61  * OT ) / OT )   ] ;
text( ws_upp_lim_m_s - right_margin_displace , ws_low_lim_m_s + 2.70 , lin_eq_str_61 , 'Color' , 'red' ) ;

r_sq_str_68                                     = ['WC68, R^2 = ' , num2str( round( r_sq_68 * OT ) / OT ), ', N = ' , num2str( num_pts_68 ) ] ;
text( ws_upp_lim_m_s - right_margin_displace , ws_low_lim_m_s + 1.50 , r_sq_str_68 , 'Color' , 'blue' ) ;
lin_eq_str_68                                   = [ '(CU), m = ' , num2str( round( m_slope_68 * OT ) / OT ) , ', b = ' , num2str( round( b_int_68  * OT ) / OT ) ] ;
text( ws_upp_lim_m_s - right_margin_displace , ws_low_lim_m_s + 0.70 , lin_eq_str_68 , 'Color' , 'blue' ) ;

%% Save Scatterplot

% Set the output plot directory:

outDir                                          = 'E:\XPIA_Research\BAO_Tower_Analysis_long_dataset\plots\' ;

% Save the plot created:

saveas( gcf , strcat( outDir , 'WC_v2_v1_scatter_STA_' , num2str( height_examined_m ) , 'm' ) , 'png' ) ;

end % End of the winds_scatter_wc_compare_regTA() function