function Ningaloo_Statistics(info, dirs)
%% NINGALOO REEF: SKILL STATISTICS
% This script calculates and plots skill statistics for the simulation
% against the measured data.

%% LOAD DATA
Meas = load([dirs.data filesep '1DStormData']);
Mod = load([dirs.data filesep 'data_LOWE'],...
    'Hsig_swell',...
    'Hsig_IG',...
    'MWL');

%% EXTRACT DATA OF INTEREST
Meas.Hsig_swell = Meas.Hsig_SW(:,7:115);
Meas.Hsig_IG = Meas.Hsig_IG(:,7:115);
Meas.MWL = Meas.MWL(:,7:115);

Mod.Hsig_swell = Mod.Hsig_swell(:,133:241);
Mod.Hsig_IG = Mod.Hsig_IG(:,133:241);
Mod.MWL = Mod.MWL(:,133:241);

%% OUTPUT VARIABLES
RMSE = ones(5,3)*NaN;
BIAS = ones(5,3)*NaN;

%% INSTRUMENT C1
[~, RMSE(1,1), ~, BIAS(1,1)] = CET_SkillStatistics(Mod.Hsig_swell(1,:), Meas.Hsig_swell(1,:));
[~, RMSE(1,2), ~, BIAS(1,2)] = CET_SkillStatistics(Mod.Hsig_IG(1,:), Meas.Hsig_IG(1,:));
[~, RMSE(1,3), ~, BIAS(1,3)] = CET_SkillStatistics(Mod.MWL(1,:), Meas.MWL(1,:));

%% INSTRUMENT C3
[~, RMSE(2,1), ~, BIAS(2,1)] = CET_SkillStatistics(Mod.Hsig_swell(2,:), Meas.Hsig_swell(2,:));
[~, RMSE(2,2), ~, BIAS(2,2)] = CET_SkillStatistics(Mod.Hsig_IG(2,:), Meas.Hsig_IG(2,:));
[~, RMSE(2,3), ~, BIAS(2,3)] = CET_SkillStatistics(Mod.MWL(2,:), Meas.MWL(2,:));

%% INSTRUMENT C4
[~, RMSE(3,1), ~, BIAS(3,1)] = CET_SkillStatistics(Mod.Hsig_swell(3,:), Meas.Hsig_swell(3,:));
[~, RMSE(3,2), ~, BIAS(3,2)] = CET_SkillStatistics(Mod.Hsig_IG(3,:), Meas.Hsig_IG(3,:));
[~, RMSE(3,3), ~, BIAS(3,3)] = CET_SkillStatistics(Mod.MWL(3,:), Meas.MWL(3,:));

%% INSTRUMENT C5
[~, RMSE(4,1), ~, BIAS(4,1)] = CET_SkillStatistics(Mod.Hsig_swell(4,:), Meas.Hsig_swell(4,:));
[~, RMSE(4,2), ~, BIAS(4,2)] = CET_SkillStatistics(Mod.Hsig_IG(4,:), Meas.Hsig_IG(4,:));
[~, RMSE(4,3), ~, BIAS(4,3)] = CET_SkillStatistics(Mod.MWL(4,:), Meas.MWL(4,:));

%% INSTRUMENT C6
C6_Mod_Hsig_swell   = Mod.Hsig_swell(5,2:2:end);
C6_Mod_Hsig_IG      = Mod.Hsig_IG(5,2:2:end);
C6_Mod_MWL          = Mod.MWL(5,2:2:end);
C6_Meas_Hsig_swell  = Meas.Hsig_swell(5,2:2:end);
C6_Meas_Hsig_IG     = Meas.Hsig_IG(5,2:2:end);
C6_Meas_MWL         = Meas.MWL(5,2:2:end);

[~, RMSE(5,1), ~, BIAS(5,1)] = CET_SkillStatistics(C6_Mod_Hsig_swell, C6_Meas_Hsig_swell);
[~, RMSE(5,2), ~, BIAS(5,2)] = CET_SkillStatistics(C6_Mod_Hsig_IG, C6_Meas_Hsig_IG);
[~, RMSE(5,3), ~, BIAS(5,3)] = CET_SkillStatistics(C6_Mod_MWL, C6_Meas_MWL);

%% PLOT STATISICS RESULTS (BAR PLOT)
figJGR = figure;
colormap([1 1 1; 0 0 0; .5 .5 .5]);
% BIAS PLOT
Axes1 = subplot (1,2,1);
hold(Axes1);
bar(BIAS);
set(Axes1,'xlim',[0 6],...
    'xtick',[1 2 3 4 5],...
    'xticklabel', ['C1';'C3';'C4';'C5';'C6'],...
    'ylim',[-0.15 0.15],...
    'ytick',[-0.15:0.05:0.15]);
xlabel('Instruments');
ylabel('Bias (m)');
text(0.25,0.13,'(a)');

% RMSE PLOT
Axes2 = subplot (1,2,2);
hold(Axes2);
bar(RMSE);
set(Axes2,'xlim',[0 6],...
    'xtick',[1 2 3 4 5],...
    'xticklabel', ['C1';'C3';'C4';'C5';'C6'],...
    'ylim',[0 0.15],...
    'ytick',[0:0.05:0.15]);
xlabel('Instruments');
ylabel('RMSE (m)');
legend('Swell','IG','z_{s}','Location','SouthEastOutside');
text(0.25,0.14,'(b)');

% Save as Fig File
saveas(figJGR,'JGR_Statistics');

%% FORMAT FOR JGR JOURNAL AND MAKE TIFF FILE
display 'Formatting for JGR Journal and saving file'
JGR_FigFormat

% RESIZE PLOTS
set(0, 'DefaultAxesXGrid','off');
set(0, 'DefaultAxesYGrid','on');

set(figJGR, 'PaperSize', [16 5]);
set(figJGR, 'PaperPosition', [0 0 16 5]);

set(Axes1,'Position',[3 1 3.5 3.5]);
set(Axes2,'Position',[8 1 3.5 3.5]);

annotation(figJGR,'line',[1 1],[1 1],...
    'Units','Centimeters',...
    'Position',[13.2 1.45 0.2 0]);

xb_write_plot(figJGR, dirs.output, 'JGR_Statistics');
end