clc;
clear all;
close all

%% Wind rose
load('windRose.mat')
figure;
intervals = fliplr([[3.5:2:max(windSpeed_80m)-1],20]);
pax = polaraxes;
colores = parula(length(intervals)+2);
for i=1:length(intervals)-1
    name = strcat(num2str(intervals(i+1)),'-',num2str(intervals(i)));
    polarhistogram(deg2rad(windDirection_80m(windSpeed_80m<intervals(i))),deg2rad(-5:10:355),'FaceColor',colores(i+1,:),'FaceAlpha',1,'displayname',name)
    hold on
end
l = legend('Show');
h = polarhistogram(deg2rad(windDirection_80m(windDirection_80m<=260 & windDirection_80m>=100)),deg2rad(-5:10:355),'LineWidth',2);
h.DisplayStyle = 'stairs';
h.EdgeColor = [0 0 0];
pax.ThetaDir = 'clockwise';
pax.ThetaZeroLocation = 'top';
%l = legend('Show')
box(l,'off')
title(l,'Wind speed (m s^{-1})')
title('Wind Rose 80m')
rticks([0:800:2400]);
rticklabels({' ','2.5%','5%','7.5%'});
hold off;
set(gca,'FontSize',16)

%% Probability density functions
load('pdf_dirShear');

figure;
subplot(1,2,1)
plot(edgesDirShear/80,distr,'Color',[0 0 0],'LineWidth',1.5);
xlim([-0.5,1])
ylim([0,0.12])
xlabel('Direction shear (deg m^{-1})');
ylabel('Probability of occurrence');
box on;
set(gca,'FontSize',14);

subplot(1,2,2)
clear all;
load('pdf_alpha')
hold on;
plot([1/7,1/7],[0,0.08],'--k')
plot(edges_alpha_pdf,pdf_alpha,'Color',[0,0,0],'LineWidth',1.5)
ht = text(1/7 + 0.05,0.05,'\alpha = 1/7');
set(ht,'Rotation',90)
xlabel('Speed shear (\alpha)')
ylabel('Probability of occurrence')
box on
set(gca,'FontSize',14);

%% Shear evolution with height
clear all;
load('evolHeight_dirShear')
figure;
subplot(1,2,1)
hold on;
barh(heights-10,mean(veeringBacking(veeringBacking_120_40>0,:)/20),'group');
barh(heights-10,mean(veeringBacking(veeringBacking_120_40<0,:)/20),'group');
errorbar(mean(veeringBacking(veeringBacking_120_40<0,:)/20),heights-8,std(veeringBacking(veeringBacking_120_40<0,:)/20),'horizontal','Color',[0.2 0.2 0.2],'LineStyle','none')
errorbar(mean(veeringBacking(veeringBacking_120_40>0,:)/20),heights-12,std(veeringBacking(veeringBacking_120_40>0,:)/20),'horizontal','Color',[0 0 0],'LineStyle','none')
xlabel('{\Delta}Wind direction over 20-m layer depths (deg m^{-1})');
ylabel('Height (m)');
legend('Veering','Backing');
yticks([40:20:120])
ylim([38 122]);
box on;
set(gca,'FontSize',14);

clear all;
load('evolHeight_alpha')
alpha_heights = [zeros(length(alpha_60_40),1),alpha_60_40,alpha_80_60,alpha_100_80,alpha_120_100];
subplot(1,2,2)
hold on;
b1 = barh(heights-10,mean(alpha_heights(alpha_120_40>0,:)),'group');
b2 = barh(heights-10,mean(alpha_heights(alpha_120_40<0,:)),'group');
b1.FaceColor = [0.47,0.67,0.19];
b2.FaceColor = [0.93,0.69,0.13];
errorbar(mean(alpha_heights(alpha_120_40<0,:)),heights-8,std(alpha_heights(alpha_120_40<0,:)),'horizontal','Color',[0.2 0.2 0.2],'LineStyle','none')
errorbar(mean(alpha_heights(alpha_120_40>0,:)),heights-12,std(alpha_heights(alpha_120_40>0,:)),'horizontal','Color',[0 0 0],'LineStyle','none')
xlabel('Speed shear over 20-m layer depths (\alpha)');
ylabel('Height (m)');
legend('\alpha > 0','\alpha < 0');
yticks([40:20:120])
ylim([38 122]);
box on;
set(gca,'FontSize',14);

%% Speed and direction shear evolution throughout the day
clear all;
load('dirShear_throughoutDay');
figure; 
hold on;
p1 = plot(edges,[meansFiltered(1:length(meansFiltered)-1),meansFiltered(1)],'k-.','LineWidth',1.75);
%e1 = errorbar(edges,[meansFiltered(1:length(meansFiltered)-1),meansFiltered(1)],[stdFiltered(1:length(stdFiltered)-1),stdFiltered(1)],'k-.');
plot(x_start,y_start/9,'Color','k','LineStyle','--');
plot(x_end,y_end/9,'Color','k','LineStyle','--');
xlim([0 23]);
ylim([0.025 0.45]);
xlabel('Local time');
ylabel('Direction shear (deg m^{-1})');
box on
% y = [[0.3,0.3],[0,0]]; % make closed patch
% x = [[5.5,9],[9,5.5]];
% patch(x,y,[0.8,0.8,0.8],'EdgeColor','none');
% alpha(0.2)
set(gca,'FontSize',14);

load('alpha_throughoutDay');
yyaxis right;
%e2 = errorbar(edges,[meansAlpha(1:length(meansAlpha)-1),meansAlpha(1)],[stdAlpha(1:length(stdAlpha)-1),stdAlpha(1)],'k:')
p2 = plot(edges,[meansAlpha(1:length(meansAlpha)-1),meansAlpha(1)],'k:','LineWidth',2.5);
% plot(x_start,y_start/9,'Color','k','LineStyle','--');
% plot(x_end,y_end/9,'Color','k','LineStyle','--');
% xlim([0 24]);
ylim([0 0.47]);
ylabel('Speed shear (\alpha)');
box on
y = [[0.466,0.466],[0,0]]; % make closed patch
x = [[6,9],[9,6]];
patch(x,y,[0.8,0.8,0.8],'EdgeColor','none');
alpha(0.2) 
% For NAWEA----
 y = [[0.466,0.466],[0,0]]; % make closed patch
 x = [[19.5,22.5],[22.5,19.5]];
 patch(x,y,[0.8,0.8,0.8],'EdgeColor','none');
 alpha(0.2)
%--------------
set(gca,'YColor',[0,0,0]);
legend([p1,p2],{'Direction shear','Speed shear'})
set(gca,'FontSize',14);

%% Speed and direction shear evolution with wind speed
clear all;
load('dirShear_windSpeed')
figure;
p1 = plot(edges(numLow:numHigh),medians_allTime(numLow:numHigh)/80,'LineWidth',1.95,'Color',[0.47,0.67,0.19]);
hold on;
%plot(edges(numLow:numHigh),means_allTime(numLow:numHigh)/80,'LineWidth',1.95,'LineStyle',':','Color',p1.Color);
p2 = plot(edges(numLow:numHigh),medians_daytime(numLow:numHigh)/80,'LineWidth',1.95,'Color',[0.93,0.69,0.13]);
%plot(edges(numLow:numHigh),means_daytime(numLow:numHigh)/80,'LineWidth',1.95,'LineStyle',':','Color',p2.Color);
p3 = plot(edges(numLow:numHigh),medians_nighttime(numLow:numHigh)/80,'LineWidth',1.95,'Color',[0.00,0.45,0.74]);
%plot(edges(numLow:numHigh),means_nighttime(numLow:numHigh)/80,'LineWidth',1.95,'LineStyle',':','Color',p3.Color);
%legend([p1 p2,p3],{'Median wind veer','Median daytime wind veer','Median nighttime wind veer'});
xlabel('Wind speed (m s^{-1})');
ylabel('Direction shear (deg m^{-1})');
%xlim([filterWSLow filterWSHigh]);
xlim([3.5,11])
xticks([filterWSLow:1:filterWSHigh]);
set(gca,'FontSize',14);

load('alpha_windSpeed')
%figure;
yyaxis right;
ps1 = plot(edges(numLow:numHigh),medians_allTime(numLow:numHigh),'LineWidth',1.95,'LineStyle',':','Color',p1.Color);
hold on;
%plot(edges(numLow:numHigh),means_allTime(numLow:numHigh),'LineWidth',1.95,'LineStyle',':','Color',p1.Color);
ps2 = plot(edges(numLow:numHigh),medians_daytime(numLow:numHigh),'LineWidth',1.95,'LineStyle',':','Color',p2.Color);
%plot(edges(numLow:numHigh),means_daytime(numLow:numHigh),'LineWidth',1.95,'LineStyle',':','Color',p2.Color);
ps3 = plot(edges(numLow:numHigh),medians_nighttime(numLow:numHigh),'LineWidth',1.95,'LineStyle',':','Color',p3.Color);
%plot(edges(numLow:numHigh),means_nighttime(numLow:numHigh),'LineWidth',1.95,'LineStyle',':','Color',p3.Color);
%legend([p1 p2,p3],{'Median speed shear','Median daytime speed shear','Median nighttime speed shear'});
%xlabel('Wind speed (m s^{-1})');
ylabel('Speed shear (\alpha)');
%xlim([filterWSLow filterWSHigh]);
%xticks([filterWSLow:1:filterWSHigh]);
set(gca,'YColor',[0,0,0]);
legend([p1 p2,p3],{'Median whole-day','Median daytime','Median nighttime'});
set(gca,'FontSize',14);


%% Speed and direction shear distributions for wind speed bins
clear all;
load('dirShear_WSbins_pdf_cummPdf');
figure;
hold on;
colores = parula(12);
p1 = plot(distrIntervals,distr(:,edges==3.5),'LineWidth',1.5,'Color',colores(10,:));
p2 = plot(distrIntervals,distr(:,edges==5),'LineWidth',1.5,'Color',colores(8,:));
p3 = plot(distrIntervals,distr(:,edges==6.5),'LineWidth',1.5,'Color',colores(6,:));
p4 = plot(distrIntervals,distr(:,edges==8),'LineWidth',1.5,'Color',colores(4,:));
p5 = plot(distrIntervals,distr(:,edges==9.5),'LineWidth',1.5,'Color',colores(2,:));
hold off;
xlim([0,0.7]);
xlabel('Direction shear (deg m^{-1})');
ylabel('Probability of occurrence');
yyaxis right
ylabel('Cummulative probability');
hold on;
colores = parula(12);
p6 = plot(distrIntervals,cummDistr(:,edges==3.5),'LineWidth',1.5,'Color',colores(10,:),'LineStyle','--');
p7 = plot(distrIntervals,cummDistr(:,edges==5),'LineWidth',1.5,'Color',colores(8,:),'LineStyle','--');
p8 = plot(distrIntervals,cummDistr(:,edges==6.5),'LineWidth',1.5,'Color',colores(6,:),'LineStyle','--');
p9 = plot(distrIntervals,cummDistr(:,edges==8),'LineWidth',1.5,'Color',colores(4,:),'LineStyle','--');
p10 = plot(distrIntervals,cummDistr(:,edges==9.5),'LineWidth',1.5,'Color',colores(2,:),'LineStyle','--');
p6.Color(4) = 0.5;
p7.Color(4) = 0.5;
p8.Color(4) = 0.5;
p9.Color(4) = 0.5;
p10.Color(4) = 0.5;
plot([0.1875,0.1875],[0,1],':k','LineWidth',1.5);
ylim([0,1]);
box on;
l = legend([p1,p2,p3,p4,p5],{'3.5-5','5-6.5','6.5-8','8-9.5','9.5-11'});
title(l,'Wind speed (m s^{-1})');
set(gca,'FontSize',14);
set(gca,'YColor',[0,0,0]);

clear all;
load('alpha_WSbins_pdf_cummPdf');
figure;
hold on;
colores = parula(12);
p1 = plot(distrIntervals,distr(:,edges==3.5),'LineWidth',1.5,'Color',colores(10,:));
p2 = plot(distrIntervals,distr(:,edges==5),'LineWidth',1.5,'Color',colores(8,:));
p3 = plot(distrIntervals,distr(:,edges==6.5),'LineWidth',1.5,'Color',colores(6,:));
p4 = plot(distrIntervals,distr(:,edges==8),'LineWidth',1.5,'Color',colores(4,:));
p5 = plot(distrIntervals,distr(:,edges==9.5),'LineWidth',1.5,'Color',colores(2,:));
hold off;
xlim([-0.3,0.8]);
ylim([0,0.4]);
xlabel('Speed shear (\alpha)');
ylabel('Probability of occurrence');
yyaxis right
ylabel('Cummulative probability');
hold on;
colores = parula(12);
p6 = plot(distrIntervals,cummDistr(:,edges==3.5),'LineWidth',1.5,'Color',colores(10,:),'LineStyle','--');
p7 = plot(distrIntervals,cummDistr(:,edges==5),'LineWidth',1.5,'Color',colores(8,:),'LineStyle','--');
p8 = plot(distrIntervals,cummDistr(:,edges==6.5),'LineWidth',1.5,'Color',colores(6,:),'LineStyle','--');
p9 = plot(distrIntervals,cummDistr(:,edges==8),'LineWidth',1.5,'Color',colores(4,:),'LineStyle','--');
p10 = plot(distrIntervals,cummDistr(:,edges==9.5),'LineWidth',1.5,'Color',colores(2,:),'LineStyle','--');
p6.Color(4) = 0.5;
p7.Color(4) = 0.5;
p8.Color(4) = 0.5;
p9.Color(4) = 0.5;
p10.Color(4) = 0.5;
plot([0.225,0.225],[0,1],':k','LineWidth',1.5);
ylim([0,1]);
box on;
l = legend([p1,p2,p3,p4,p5],{'3.5-5','5-6.5','6.5-8','8-9.5','9.5-11'});
title(l,'Wind speed (m s^{-1})');
disp('change of convexity of cummulative probaility distributions for speed shear')
set(gca,'YColor',[0,0,0]);
set(gca,'FontSize',14);

%% Relationship between speed and direction shear
clear all
load Rel_alpha_dirShear.mat

figure;
%subplot(2,1,1)
colores = parula(12);
hold on;
pp1 = plot(edges_alpha,mean_relDirShear_alphaColumns_3_5_WS_5_veer,'LineStyle','-','LineWidth',1.2,'Color',colores(10,:));
pp2 = plot(edges_alpha,mean_relDirShear_alphaColumns_5_WS_6_5_veer,'LineStyle','-','LineWidth',1.2,'Color',colores(8,:));
pp3 = plot(edges_alpha,mean_relDirShear_alphaColumns_6_5_WS_8_veer,'LineStyle','-','LineWidth',1.2,'Color',colores(6,:));
pp4 = plot(edges_alpha,mean_relDirShear_alphaColumns_8_WS_9_5_veer,'LineStyle','-','LineWidth',1.2,'Color',colores(4,:));
pp5 = plot(edges_alpha,mean_relDirShear_alphaColumns_9_5_WS_11_veer,'LineStyle','-','LineWidth',1.2,'Color',colores(2,:));

plot(edges_alpha,mean_relDirShear_alphaColumns_3_5_WS_5_back,'Color',pp1.Color,'LineStyle','--','LineWidth',1.2)
plot(edges_alpha,mean_relDirShear_alphaColumns_5_WS_6_5_back,'Color',pp2.Color,'LineStyle','--','LineWidth',1.2)
plot(edges_alpha,mean_relDirShear_alphaColumns_6_5_WS_8_back,'Color',pp3.Color,'LineStyle','--','LineWidth',1.2)
plot(edges_alpha,mean_relDirShear_alphaColumns_8_WS_9_5_back,'Color',pp4.Color,'LineStyle','--','LineWidth',1.2)
plot(edges_alpha,mean_relDirShear_alphaColumns_9_5_WS_11_back,'Color',pp5.Color,'LineStyle','--','LineWidth',1.2)

p1 = plot(edges_alpha(occurrences_veering>30),mean_alpha_dirShear_veer(occurrences_veering>30),'LineWidth',1.5,'Color',[0 0 0]);
p2 = plot(edges_alpha(occurrences_backing>30),mean_alpha_dirShear_back(occurrences_backing>30),'LineStyle','--','LineWidth',1.5,'Color',[0 0 0]);
plot([0,0],[-0.16,0.43],'--','Color',[0.65,0.65,0.65],'LineWidth',1.5);
legend([p1,p2],{'Veering','Backing'});
xlim([-0.3,0.8])
ylim([-0.16,0.43])
xlabel('Speed shear (\alpha)');
ylabel('Direction shear (deg m^{-1})');
box on;

l = legend([pp1,pp2,pp3,pp4,pp5,p1],{'3.5-5','5-6.5','6.5-8','8-9.5','9.5-11','3.5-11'});
title(l,'Wind speed (m s^{-1})')
box on;
set(gca,'FontSize',14);
