%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIS CODE READS AND PLOTS THE OUTCOMES OF CLUSTERING ANALYSIS APPLIED
% TO HINDCAST DATA AT NODE 000300 
% (http://www3.dicca.unige.it/meteocean/hindcast.html)
% Cremonini G., De Leo F., Besio G., 2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clearvars
close all
clc

%--------------------------------------------------------------------------
% LOAD HISTORICAL SERIES OF METOCEAN PARAMETERS
%--------------------------------------------------------------------------
% date -> dates of the inital time series
% H -> significant wave height series
% Tm -> mean period series
% thm -> mean incoming direction series
% UW -> wind speed series
% thw -> wind direction series
load('hind_param.mat')

%--------------------------------------------------------------------------
% LOAD TIME-VARYING SCENARIOS
%--------------------------------------------------------------------------
% Xt -> matrix made of initial series reordered in time varying scenarios
% DATEt -> dates corresponding to the scenarios
% idDc -> rows corresponding to the selected model scenarios (centroids)
load('clusters.mat')

rdt = 6;    % time resolution of the case study [hrs]

%--------------------------------------------------------------------------
% PLOT DIFFERENT SIGNIFICANT SCENARIOS
%--------------------------------------------------------------------------

% circular colorscale ! makeColorMap function required
a = makeColorMap([0 0 1],[1 0 0],[1 1 0],16);
b = makeColorMap([1 1 0],[0 1 0],[0 0 1],16);
c = [a;b];

% INTENSE WAVES MOVING NW
IT = 24;    % arg index of the selected scenario
itime = find(datenum(date)>=DATEt(idDc(IT),1) &...
                datenum(date)<=DATEt(idDc(IT),end));
            
% wave data            
figure
p1=plot(1:rdt:24,X(idDc(IT),1:4),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),H(itime),100,thm(itime),'filled');hold on
p2=plot(1:numel(itime),H(itime),'k:');
text(0.02,0.95,'A)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta [°N]','FontSize',14)
xlim([1 19])
xticks([1 7 13 19])
ylim([0 3.5])
xlabel('Time [hr]')
ylabel('Hs [m]')
legend([p1 p2],'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% wind data
figure
p1 = plot(1:rdt:24,X(idDc(IT),13:16),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),UW(itime),100,thw(itime),'filled');hold on
p2 = plot(1:numel(itime),UW(itime),'k:');
text(0.02,0.95,'A)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta_w [°N]','FontSize',14)
xlabel('Time [hr]')
ylabel('Wind [m/s]')
xlim([1 19])
xticks([1 7 13 19])
ylim([0 20])
legend([p1 p2], 'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% MILD WAVES MOVING NW
IT = 14;    % arg index of the selected scenario
itime = find(datenum(date)>=DATEt(idDc(IT),1) &...
                    datenum(date)<=DATEt(idDc(IT),end));

% wave data
figure
p1=plot(1:rdt:24,X(idDc(IT),1:4),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),H(itime),100,thm(itime),'filled');hold on
p2=plot(1:numel(itime),H(itime),'k:');
text(0.02,0.95,'B)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta [°N]','FontSize',14)
xlim([1 19])
xticks([1 7 13 19])
ylim([0 3.5])
xlabel('Time [hr]')
ylabel('Hs [m]')
legend([p1 p2],'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% wind data
figure
p1 = plot(1:rdt:24,X(idDc(IT),13:16),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),UW(itime),100,thw(itime),'filled');hold on
p2 = plot(1:numel(itime),UW(itime),'k:');
text(0.02,0.95,'B)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta_w [°N]','FontSize',14)
xlabel('Time [hr]')
ylabel('Wind [m/s]')
xlim([1 19])
xticks([1 7 13 19])
ylim([0 20])
legend([p1 p2], 'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% INTENSE WAVES MOVING NE

IT = 4; % arg index of the selected scenario
itime = find(datenum(date)>=DATEt(idDc(IT),1) &...
                datenum(date)<=DATEt(idDc(IT),end));

% wave data
figure
p1=plot(1:rdt:24,X(idDc(IT),1:4),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),H(itime),100,thm(itime),'filled');hold on
p2=plot(1:numel(itime),H(itime),'k:');
text(0.02,0.95,'C)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta [°N]','FontSize',14)
xlim([1 19])
xticks([1 7 13 19])
ylim([0 3.5])
xlabel('Time [hr]')
ylabel('Hs [m]')
legend([p1 p2],'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% wind data
figure
p1 = plot(1:rdt:24,X(idDc(IT),13:16),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),UW(itime),100,thw(itime),'filled');hold on
p2 = plot(1:numel(itime),UW(itime),'k:');
text(0.02,0.95,'C)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta_w [°N]','FontSize',14)
xlabel('Time [hr]')
ylabel('Wind [m/s]')
xlim([1 19])
xticks([1 7 13 19])
ylim([0 20])
legend([p1 p2], 'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% MILD WAVES MOVING NE

IT = 31;    % arg index of the selected scenario
itime = find(datenum(date)>=DATEt(idDc(IT),1) &...
                    datenum(date)<=DATEt(idDc(IT),end));

% wave data
figure
p1=plot(1:rdt:24,X(idDc(IT),1:4),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),H(itime),100,thm(itime),'filled');hold on
p2=plot(1:numel(itime),H(itime),'k:');
text(0.02,0.95,'D)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta [°N]','FontSize',14)
xlim([1 19])
xticks([1 7 13 19])
ylim([0 3.5])
xlabel('Time [hr]')
ylabel('Hs [m]')
legend([p1 p2],'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% wind data
figure
p1 = plot(1:rdt:24,X(idDc(IT),13:16),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),UW(itime),100,thw(itime),'filled');hold on
p2 = plot(1:numel(itime),UW(itime),'k:');
text(0.02,0.95,'D)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta_w [°N]','FontSize',14)
xlabel('Time [hr]')
ylabel('Wind [m/s]')
xlim([1 19])
xticks([1 7 13 19])
ylim([0 20])
legend([p1 p2], 'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% INTENSE WAVES MOVING S

IT = 22;    % arg index of the selected scenario
itime = find(datenum(date)>=DATEt(idDc(IT),1) &...
                    datenum(date)<=DATEt(idDc(IT),end));

% wave data
figure
p1=plot(1:rdt:24,X(idDc(IT),1:4),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),H(itime),100,thm(itime),'filled');hold on
p2=plot(1:numel(itime),H(itime),'k:');
text(0.02,0.95,'E)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta [°N]','FontSize',14)
xlim([1 19])
xticks([1 7 13 19])
ylim([0 3.5])
xlabel('Time [hr]')
ylabel('Hs [m]')
legend([p1 p2], 'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% wind data
figure
p1 = plot(1:rdt:24,X(idDc(IT),13:16),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),UW(itime),100,thw(itime),'filled');hold on
p2 = plot(1:numel(itime),UW(itime),'k:');
text(0.02,0.95,'E)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta_w [°N]','FontSize',14)
xlabel('Time [hr]')
ylabel('Wind [m/s]')
xlim([1 19])
xticks([1 7 13 19])
ylim([0 20])
legend([p1 p2],'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% MILD WAVES MOVING S

IT = 43;    % arg index of the selected scenario
itime = find(datenum(date)>=DATEt(idDc(IT),1) &...
                datenum(date)<=DATEt(idDc(IT),end));

% wave data
figure
p1=plot(1:rdt:24,X(idDc(IT),1:4),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),H(itime),100,thm(itime),'filled');hold on
p2=plot(1:numel(itime),H(itime),'k:');
text(0.02,0.95,'F)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta [°N]','FontSize',14)
xlim([1 19])
xticks([1 7 13 19])
ylim([0 3.5])
xlabel('Time [hr]')
ylabel('Hs [m]')
legend([p1 p2],'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)

% wind data
figure
p1 = plot(1:rdt:24,X(idDc(IT),13:16),'k-','LineWidth',1.4);hold on
scatter(1:numel(itime),UW(itime),100,thw(itime),'filled');hold on
p2 = plot(1:numel(itime),UW(itime),'k:');
text(0.02,0.95,'F)','Units','Normalized','FontWeight','Bold','FontSize',24);hold on
colormap(c)
caxis([0 360])
h=colorbar;
set(h,'YTick',0:45:360)
ylabel(h, '\theta_w [°N]','FontSize',14)
xlabel('Time [hr]')
ylabel('Wind [m/s]')
xlim([1 19])
xticks([1 7 13 19])
ylim([0 20])
legend([p1 p2],'Centroid','Real')
legend boxoff
set(gca,'FontSize',16)
