function Weibull_Distribution_WCs( flag_wc16_avg_exact , flag_v2_std_or_mc , WC_nums , avg_pd_min )
%Weibull_Distribution_WCs Create a Weibull distribution of tower wind speeds for each wind farm.
%   Create a Weibull distribution for the tower using the built in Matlab packages.

% Example function call:
% Weibull_Distribution_WCs( 'exact' , 'mc' , 2 : 4 , 2 )

%% Set parameters

fss                                     = 'FontSize';

font_size                               = 14 ;

% Note the relevant heights for examination:

height_m_BAO_lower                      = 50 ;

height_m_BAO_middle                     = 100 ;

height_m_BAO_upper                      = 150 ;

heights                                 = [height_m_BAO_lower ; height_m_BAO_middle ; height_m_BAO_upper ] ;

% 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 relevant data from WindCubes

% Set in directory for WindCube data:

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

% Preallocate the WCs array info:

WCnum                                   = NaN(  length( WC_nums ) , 1 ) ; 
WCnum_str                               = cell( length( WC_nums ) , 1 ) ; 

% Loop over the WCs and load the data, and store the information:

for j_wc = WC_nums
    
    % Assign information base upon Wind Cube Number:
    
    switch j_wc
        case 1
            WCnum( j_wc )               = 49 ;
            WCnum_str{ j_wc }           = [ 'WC' , num2str( WCnum( j_wc ) ) , ' v1 (CU)' ] ;
        case 2
            WCnum( j_wc )               = 61 ;
            WCnum_str{ j_wc }           = [ 'WC' , num2str( WCnum( j_wc ) ) , ' v1 (NCAR)' ] ;
        case 3
            WCnum( j_wc )               = 68 ;
            WCnum_str{ j_wc }           = [ 'WC' , num2str( WCnum( j_wc ) ) , ' v1 (CU)' ] ;
        case 4
            WCnum( j_wc )               = 16 ;
            WCnum_str{ j_wc }           = [ 'WC' , num2str( WCnum( j_wc ) ) , ' v2 (NRG)' ] ;
    end
        
	load( strcat( inDir_WC , curr_exper_abbr , '_' , num2str( WCnum( j_wc ) ) , '_' , 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 curr_exper_abbr ...
    height* heights avg_pd_min flag_wc16_avg_exact ...
    fss text_fs font_size datestr_formatOut...
    WCnum WC_nums WCnum_str

%% 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_m_BAO_lower - 10 ) ; % Index 1 is 40m
ws_WC68_40m_m_s                         = ws_case_68( : , alt_68 == height_m_BAO_lower - 10 ) ; % Index 1 is 40m
ws_WC16_40m_m_s                         = ws_case_16( : , alt_16 == height_m_BAO_lower - 10 ) ; % Index 1 is 40m

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

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

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

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

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

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

%% Horizontally concatenate WCs, then take mean (being restrictive) of heights above & below height of interest if needed to compare against BAO Tower:

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 as the tower timeseries:

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 ) ;

%% Setup the matrices of what was read in:

% The matricies of wind speed are dimensions time x height analyzed for each WC device:

ws_WC61_m_sec                           = [ ws_WC61_50m_m_s ws_WC61_100m_m_s ws_WC61_150m_m_s ] ;
ws_WC68_m_sec                           = [ ws_WC68_50m_m_s ws_WC68_100m_m_s ws_WC68_150m_m_s ] ;
ws_WC16_m_sec                           = [ ws_WC16_50m_m_s ws_WC16_100m_m_s ws_WC16_150m_m_s ] ;

% Clear out the excess memory:

clear *50m_m_sec *100m_m_sec *150m_m_sec

%% Loop over all tower booms to make separate plots:

for j_wc = WC_nums
    %% Loop over all the heights to make separate plots:
    
    if     j_wc == 2
        ws_m_sec                        = ws_WC61_m_sec ;
    elseif j_wc == 3
        ws_m_sec                        = ws_WC68_m_sec ;
    elseif j_wc == 4
        ws_m_sec                        = ws_WC16_m_sec ;
    end % End if statement on choice of boom
    
    for i_height = 1 : length( heights )
        
        % Set wind speed to the current wind speed:
        
        wind_speed_m_s                  = ws_m_sec( : , i_height ) ;
        
        %% Weibull Code from Dr. Matt Aitken!
        
        % This is a lazy copy so we do not have to change the u in the code below:
        % (Also does not take up any more RAM since u points to the wind speed variable that came in)
        
        u                               = wind_speed_m_s ;
        
        %% Weibull Calculations
        
        % Throw out bad measurements, only taking non-NaN and nonzero measurements (a la Dr. Matt Aitken):
        u                               = u( ~isnan( u ) & u ~= 0  ) ;
        
        % Setup the bins and the bin centers
        uStep                           = 0.5 ; % Step size between the bins
        uBin                            = uStep / 2 : uStep : ceil( max( u ) ) ; % ceil( max( u ) ) rounds the maxium upward
        
        % relative frequency per uStep interval at bin centers (will make the area under the Weibull equal to 1)
        nU                              = hist( u , uBin ) / ( length( u ) * uStep ) ; % normalizes the relative frequency by the bin interval (Note: length(u)* uStep gives area)
        
        % Maximum likelihood estimates for Weibull parameters from wblfit()
        p                               = wblfit( u ) ; % p(1)--> Scale parameter ("a" or "lambda"), p(2) --> Shape parameter ("b" or "k")
        
        % Calculate the best-fit Weibull pdf given the scale and shape parameters:
        n_pts_calc_fit                  = 1000 ;
        xFit                            = linspace( 0 , max( u ) , n_pts_calc_fit ) ; % Evenly spaced x values to calculate a fit curve
        yFit                            = wblpdf( xFit , p( 1 ) , p( 2 ) ) ; % This is the best Weibull fit curve given the scale and shape parameters
        
        %% Plot the wind speed histogram
        
        % Create figure
        fig_name                        = strcat( 'BAO_Tower_Weibull_Distribution' ) ;
        
        figure( 'Name' , fig_name, 'units' , 'normalized' , 'outerposition', [ 0 0 1 1 ] ) ;
        
        % Set the Weibull plot settings:
        lws                             = 'LineWidth';
        lw                              = 2;
        
        % Set the horizontal and vertical axes maxima:
        
        ws_horiz_lim_max                = 18 ;
        prob_vert_lim_max               = 0.18 ;
        
        switch j_wc
            case 2 % WC61
                switch i_height
                    case 1 % 50m
                        bar( uBin, nU, 'FaceColor', 'r' , 'EdgeColor', 'w' ); % Draw the unique color histogram as a bar graph with the bins and relative frequency of the bins
                    case 2 % 100m
                        bar( uBin, nU, 'FaceColor', 'r' , 'EdgeColor', 'w' ); % Draw the unique color histogram as a bar graph with the bins and relative frequency of the bins
                    case 3 % 150m
                        bar( uBin, nU, 'FaceColor', 'r' , 'EdgeColor', 'w' ); % Draw the unique color histogram as a bar graph with the bins and relative frequency of the bins
                end
            case 3 % WC68
                switch i_height
                    case 1 % 50m
                        bar( uBin, nU, 'FaceColor', 'b' , 'EdgeColor', 'w' ); % Draw the unique color histogram as a bar graph with the bins and relative frequency of the bins
                    case 2 % 100m
                        bar( uBin, nU, 'FaceColor', 'b' , 'EdgeColor', 'w' ); % Draw the unique color histogram as a bar graph with the bins and relative frequency of the bins
                    case 3 % 150m
                        bar( uBin, nU, 'FaceColor', 'b' , 'EdgeColor', 'w' ); % Draw the unique color histogram as a bar graph with the bins and relative frequency of the bins
                end
            case 4 % WC16
                switch i_height
                    case 1 % 50m
                        bar( uBin, nU, 'FaceColor', 'g' , 'EdgeColor', 'w' ); % Draw the unique color histogram as a bar graph with the bins and relative frequency of the bins
                    case 2 % 100m
                        bar( uBin, nU, 'FaceColor', 'g' , 'EdgeColor', 'w' ); % Draw the unique color histogram as a bar graph with the bins and relative frequency of the bins
                    case 3 % 150m
                        bar( uBin, nU, 'FaceColor', 'g' , 'EdgeColor', 'w' ); % Draw the unique color histogram as a bar graph with the bins and relative frequency of the bins
                end
                
        end
        
        %% Plot the Weibull fit curve
        
        hold on ; % Keep the figure there, and overplot the fit line
        
        plot( xFit , yFit , 'k' , lws , lw ) ; % Overplot the Weibull fit
        
        set(gca,'XLim',[0 ws_horiz_lim_max],'YLim',[0 prob_vert_lim_max], fss, font_size, 'FontWeight','bold');
        
        xlabel('Wind Speed [m s^{-1}]', 'FontWeight', 'bold' );
        
        ylabel(['relative freq. per ' num2str(uStep) ' m s^{-1} interval'], 'FontWeight', 'bold');
        
        %% Add the legend, noting the extra spacing:
        
        legend( 'Observed' , [ 'Weibull fit:            '; '\lambda = ' num2str( p( 1 ) , '%4.2f' ) ', k = ' num2str( p( 2 ) , '%4.2f' ) ] , 'Location' , 'Northeast' ) ;
        
        %% Set the title depending upon the WC device
        
        title_fontsize              = 17.5 ;
        
        title_top_str               = strcat( {'2-min Avg. Winds Histogram & Weibull: '} , char( WCnum_str{ j_wc } ) , {', '} , num2str( heights( i_height ) ) , {'m Hgt.'} ) ;
        
        dates_str                   = [ datestr( lower_datetime , datestr_formatOut ) , ' to ' , datestr( upper_datetime , datestr_formatOut ) ] ;
        title( { char( title_top_str ) ; dates_str } , fss , title_fontsize , 'FontWeight' , 'bold' ) ; % Create the unique plot title:
        
        %% Save the figure:
        
        % Set the output plot directory:
        
        outDir                                  = 'E:\XPIA_Research\BAO_Tower_Analysis_long_dataset\plots\' ;
        
        % Save the plot created:
        saveas( gcf , strcat( outDir , 'WC_Hist_and_Weibull_dist_WC' ,  num2str( ( WCnum( j_wc) ) ) , '_' ,  num2str( heights( i_height ) ) , 'm' ) , 'png' ) ;
        
    end % End loop over heights (i_height)
    
end % End loop over WC Numbers (j_wc)

end % End of Weibull_Distribution_WCs() function

