%% README
% ======================================================================= %
% The file is constructed as follows:
% get_paths_and_types(local) retrieves the links to the files in the folder
% "HKN_12_to_15h_data", which is given as a .zip and NEEDS TO BE UNPACKED!
% Then, the user can select which cases to load. This is done via an ID
% given to each controller, below you find the list of avaiable controllers
% with their ID.
% The data of the selected codes is loaded and processed. The data is
% mainly composed of the power generated, from which we calculate the
% energy generated across the entire time, as well as a sliding window. The
% length of the sliding window can be set! In the publication 600s is used.
%
% User defined variables:
%   controller_ids          : Controllers to evaluate and load data from
%   time_window             : Sliding window length for energy and yaw
%                               travel
%   norm_id                 : Index of the controller data that should be
%                               used for normalization: -> 1 means the
%                               first controller in controller_ids is used.
%   controller_ids_tmp      : Used to only plot a selection of the loaded
%                               controllers.
%   iT                      : Index of the turbine to plot. This index
%                               starts at 1 and ends with 10, contrairy to
%                               the turbine names, thus to plot T4 -> iT=5
% ======================================================================= %
%% Evaluate long cases
local = true;   % Needs to be set to true 
cases = get_paths_and_types(local);
% Types
% % --------------------------------
% cases.names{1} =  'Baseline  2 deg 0.01 ki, no noise, 5s';       % 1 internal
% cases.names{3} =  'Baseline  5 deg 0.02 ki, no noise, 5s';       % 3 off
%
% cases.names{5} =  'Baseline  2 deg 0.01 ki, with noise, 60s';    % 5 internal
% cases.names{6} =  'Baseline  4 deg 0.01 ki, with noise, 60s';    % 6 internal
% % -------------------------------- internal LuT
% cases.names{11} =  'LuT int. 2 deg 0.01 ki, no noise, 5s';       % 11
% cases.names{12} =  'LuT int. 4 deg 0.01 ki, no noise, 5s';       % 12
%
% cases.names{15} =  'LuT int. 2 deg 0.01 ki, with noise, 60s';    % 15
% cases.names{16} =  'LuT int. 4 deg 0.01 ki, with noise, 60s';    % 16
% % --------------------------------
% cases.names{31} =  'CLC 2mr 100 ph, no noise, 60s';                % 31
% cases.names{32} =  'CLC 2mr 200 ph, no noise, 60s';                % 32
% cases.names{33} =  'CLC 2mr shifted no preview, no noise, 60s';    % 33
% cases.names{34} =  'CLC 2mr shifted with preview, no noise, 60s';  % 34
%
% cases.names{35} =  'CLC 2mr 100 ph, with noise, 60s';              % 35
% cases.names{36} =  'CLC 2mr 200 ph, with noise, 60s';              % 36
% cases.names{37} =  'CLC 2mr shifted no preview, with noise, 60s';  % 37
% % --------------------------------

% ==== Available control results (local) ====
% 12_to_15 = [1     3     5    6    11    12    15    16    31    32    33    34    35    36    37];
% =========================================
% %% Load Wind direction
% dir      = readmatrix("HKN_12_to_15h_dir.csv");
% dir(:,1) = (dir(:,1)-1200)/3600;
% dir      = dir(3:end,:);
%% Select case and controllers to investigate
% Initially more case studies were planned, but only this case study was eventually investigated
case_ids        = 3; % Do NOT change, needs to be 3!
controller_ids  = sort([1 5 11 12 15 16 32 36 31 35 33 37 34]);

% load power
cases = load_power(cases, case_ids, controller_ids);

% Calculate total energy
cases = calculate_energy(cases, case_ids, controller_ids);

% Calculate yaw QoI
cases = load_yaw(cases, case_ids, controller_ids);

% Calculate averaged energy
time_window = 600;
cases = calculate_window_energy(cases, case_ids, controller_ids, time_window);
cases = calculate_window_yaw_travel(cases, case_ids, controller_ids, time_window);
cases.t = cases.t/(60*60);
cases.t = cases.t - cases.t(1);



%% Plot total energy
norm_id = 1;
plot_norm_energy(case_ids, controller_ids, cases, norm_id);

%% Plot normed window energy - Figure 14
controller_ids_tmp = [1 11 12 31 32 33]; % Noise free
%controller_ids_tmp = [5 15 16 35 36 37]; % Disturbed
norm_id = 1;
plot_norm_window_energy(case_ids, controller_ids_tmp, cases, norm_id, time_window);

%% Plot absolute window energy - Figure 14
controller_ids_tmp = [1 11 12 31 32 33]; % Noise free
%controller_ids_tmp = [5 15 16 35 36 37]; % Noisy
plot_window_energy(case_ids, controller_ids_tmp, cases, time_window);


%% Plot box chart - Figure 15
norm_id = 1;
plot_norm_window_energy_box(case_ids, controller_ids, cases, norm_id, time_window);

%% Normed energy over normed yaw travel (pareto curve)  - Figure 16
norm_id = 1;
plot_final_energy_over_yaw_travel(case_ids, controller_ids, cases, norm_id);


%% Bar Chart of efficiency - not used
norm_id = 1;
final_energy_eff_bar_chart(case_ids, controller_ids, cases, norm_id)


%% Turbine Bar Chart of efficiency - Figure 13 + unused
norm_id = 1;
controller_ids_tmp = [1 11 12 31 32 33]; % Noise free
%controller_ids_tmp = [5 15 16 35 36 37]; % Noisy
energy_eff_bar_chart_turbine(case_ids, controller_ids_tmp, cases, norm_id)


%% yaw travel  over time window - not used
plot_window_yaw_travel(case_ids, controller_ids, cases, time_window);

%% Effectiveness of turbines
iT = 9;
controller_ids_tmp = [1 11 12 31 32 33];
norm_id = 1;
plot_turbine_norm_window_energy(case_ids, controller_ids_tmp, cases, iT, norm_id, time_window);
% xlim([140, 170]./60)
% xticks([140:10:170]./60)
% xticklabels({'140-150','150-160','160-170','170-180'})
%% Norm power generated turbine
iT = 9;
norm_id = 1;
plot_turbine_norm_power(case_ids, controller_ids_tmp, cases, iT, norm_id);
xlim([140, 170]./60)


%% Effectiveness of turbines
iT = 5;
%controller_ids_tmp = [1 11 32 36];
plot_turbine_window_energy(case_ids, controller_ids_tmp, cases, iT, time_window);

%% Window energy turbines combined - not used
iT = [5,9];
controller_ids_tmp = [1 11 12 31 32 33 34]; % Noise free
%controller_ids_tmp = [5 15 16 35 36 37]; % Noisy
plot_turbine_window_energy_combined(case_ids, controller_ids_tmp, cases, iT, time_window)

% Time selection for Fig 12
% xlim([140, 170]./60)
% xticks([140:10:170]./60)
% xticklabels({'140-150','150-160','160-170','170-180'})


%% Normed window energy turbines combined - Figure 12
iT = [5,9];
norm_id = 1;
controller_ids_tmp = [1 11 12 31 32 33 34]; % Noise free
%controller_ids_tmp = [5 15 16 35 36 37]; % Noisy
plot_turbine_norm_window_energy_combined(case_ids, controller_ids_tmp, cases, iT, time_window, norm_id)

% Time selection for Fig 12
% xlim([140, 170]./60)
% xticks([140:10:170]./60)
% xticklabels({'140-150','150-160','160-170','170-180'})


%% Turbine individual energy normed
iT = 9;
controller_ids_tmp = [1 11 32 33 35];
norm_id = 1;
plot_turbine_diff_window_energy(case_ids, controller_ids_tmp, cases, iT, norm_id, time_window);


%% Plot yaw over time - not used
% 3D Plot! 
iT      = 3;
norm_id = 1;
controller_ids_tmp = [1 11 12 31 32 33]; % Noise free
%controller_ids_tmp = [5 15 16 35 36 37]; % Noisy

plot_yaw_turbine(case_ids, controller_ids_tmp, cases, iT, cases.dir_domain{case_ids})


%% Plot yaw over time - Figure 12 (a)
% 3D Plot! 
iT      = 5;
norm_id = 1;
controller_ids_tmp = [1 11 12 31 32 33]; % Noise free
%controller_ids_tmp = [5 15 16 35 36 37]; % Noisy

plot_ori_turbine(case_ids, controller_ids_tmp, cases, iT, cases.dir_domain{case_ids})

%% Absolute yaw misalignment CDF histogram - not used
controller_ids_tmp = [1 11 12 31 32 33]; % Noise free
%controller_ids_tmp = [5 15 16 35 36 37]; % Noisy
plot_yaw_histogram(case_ids, controller_ids_tmp, cases,cases.dir_domain{case_ids})


%% Plot LuT - Figure 7
% Load LuT
lut = readmatrix('HKN_LuT_orientation.csv');

% Convert into misalignment instead of orientation
lut = (1:2:359)' - lut;

figure
imagesc(1:2:359,0:9,lut')
colormap(flipud(RdBu(61)));
clim([-30.5,30.5])
c = colorbar(gca,"southoutside");
c.Label.String = 'Yaw misalignment [deg]';
xlabel('Wind Direction [deg]')
ylabel('Turbine iD')
