clear all;clc;close all

%% load winddata
WINDFILE=csvread('wind2010v2.txt');
dir=WINDFILE(:,4)';
v_wind=WINDFILE(:,5)'/10;
clear WINDFILE

%% calc

rhoa = 1.225;

climate=zeros(21,37);

n=0;
for i=1:8760
    %windspeed=roundn(ustar(i),-2);
    direction=dir(i);
    
    if direction<900
    
    y=v_wind(i)+1;%y=round(windspeed*100+1);
    x=round(direction/10+1);
    
    climate(y,x)=climate(y,x)+1;

    n=n+1;
    else
    end
end

climate(:,1)=climate(:,1)+climate(:,37);
climate=climate(:,1:36);

duration = climate;


%% transform wind velocity field to wind shear stress field and transport
%% capacity field

v = 0:1:20;
dir = 0:10:350;

[V,DIR] = meshgrid(v,dir);

uster = V./((1/0.4)*log(10/10^-3));
tauster = rhoa*uster.^3;
tauster_w = tauster.*duration'; %transport cap

%% reduc wind conditions

A = 0.1;
d = 250e-6;
rhos  = 2650;
g = 9.81;
ucr = A*sqrt(d*g*(rhos-rhoa)/rhoa);
taucr = rhoa*ucr^3;

%set tauster to 0/NaN if smaller than taucr
tauster(tauster<=taucr)=0;
tauster(DIR<200 & DIR>20)=0;
tauster_w(tauster<=taucr)=0;
tauster_w(DIR<200 & DIR>20)=0;

%% plot results
figure('units','normalized','outerposition',[0 0 1 1]);

subplot(2,2,1);pcolor(DIR,V,climate');title('wind climate');
subplot(2,2,2);pcolor(DIR,V,tauster_w);title('reduced effective transport');
subplot(2,2,3);plot(dir,sum(tauster_w,2),'r');title('total eff. transport/dir');

tausterS = sum(tauster_w,2);

cum_tauster=cumsum(tausterS',2);

transport=sum(tausterS',2);

%number of direction bins
nrbins=10;

y=0:transport/nrbins:transport

%x location

xbin(1)=0

for n=2:11
i=1;    
while cum_tauster(i)<y(n)
    i=i+1
end
xbin(n)=(10/(cum_tauster(i)-cum_tauster(i-1)))*(y(n)-cum_tauster(i-1))+(i-2)*10
end
    
subplot(2,2,4); plot(dir,cum_tauster,'r');title('cum. eff. transport & bands');
hold on;
for i=1:11
    line([0 350], [y(i) y(i)],'linestyle',':','color','k')
end
for i=1:11
    line([xbin(i) xbin(i)], [0 transport],'linestyle','-','color','k')
end    

%% shore relative dist
 clearvars -except climate rhoa

v = [0:1:20];
dir = [0:10:170];

duration(:,1:14)=climate(:,22:35);
duration(:,15:18)=climate(:,1:4);

[V,DIR] = meshgrid(v,dir);

uster = V./((1/0.4)*log(10/10^-3));
tauster = rhoa*uster.^3;
tauster_w = tauster.*duration'; %transport cap

%% reduc wind conditions

A = 0.1;
d = 250e-6;
rhos  = 2650;
g = 9.81;
ucr = A*sqrt(d*g*(rhos-rhoa)/rhoa);
taucr = rhoa*ucr^3;

%set tauster to 0/NaN if smaller than taucr
tauster(tauster<=taucr)=0;
tauster_w(tauster<=taucr)=0;


%% plot results

figure('units','normalized','outerposition',[0 0 1 1]);
subplot(5,4,[1,5]);pcolor(DIR,V,duration');title('wind climate', 'FontSize', 11 );

ylabel( 'Wind Speed (m/s)', 'FontSize', 10, 'Rotation', 90 ,'FontName','Arial');
xlabel( 'Wind Dir', 'FontSize', 10, 'FontName','Arial');
colorbar('YTickLabel',{'0 hours','10 hours','20 hours','30 hours','40 hours','50 hours'});


subplot(5,4,[2,6]);pcolor(DIR,V,tauster_w);title('reduced onshore effective transport', 'FontSize', 11 );
ylabel( 'Wind Speed (m/s)', 'FontSize', 10, 'Rotation', 90 ,'FontName','Arial');
xlabel( 'Wind Dir', 'FontSize', 10, 'FontName','Arial');
colorbar

figure('units','normalized','outerposition',[0 0 1 1]);

subplot(5,4,[1,5]):plot(dir,sum(tauster_w,2),'r');title('total eff. onshore transport/dir', 'FontSize', 11 );

ylabel( 'transport', 'FontSize', 10, 'Rotation', 90 ,'FontName','Arial');
xlabel( 'Wind Dir', 'FontSize', 10, 'FontName','Arial');

tausterS = sum(tauster_w,2);

cum_tauster=cumsum(tausterS',2);

transport=sum(tausterS',2);

%number of direction bins
nrbins=10

y=0:transport/nrbins:transport

%x location

xbin(1)=0

for n=2:11
i=1;    
while cum_tauster(i)<y(n)
    i=i+1
end
xbin(n)=(10/(cum_tauster(i)-cum_tauster(i-1)))*(y(n)-cum_tauster(i-1))+(i-2)*10
end
    
subplot(5,4,[2,6]); plot(dir,cum_tauster,'r');title('cum. eff. transport & onshore bands','FontSize', 11);

ylabel( 'transport', 'FontSize', 10, 'Rotation', 90 ,'FontName','Arial');
xlabel( 'Wind Dir', 'FontSize', 10, 'FontName','Arial');

hold on;
for i=1:11
    line([0 170], [y(i) y(i)],'linestyle',':','color','k')
end
for i=1:11
    line([xbin(i) xbin(i)], [0 transport],'linestyle','-','color','k')
end    

ybin=1:4:20

figure;pcolor(DIR,V,duration');shading flat;title('wind climate & bands');
ylabel( 'wind speed', 'FontSize', 10, 'Rotation', 90 ,'FontName','Arial');
xlabel( 'Wind Dir', 'FontSize', 10, 'FontName','Arial');
hold on;
for i=1:5
    line([0 170], [ybin(i) ybin(i)],'color','w')
end
for i=1:11
    line([xbin(i) xbin(i)], [0 transport],'linestyle','-','color','w')
end

%% merge bins info

ybin=10:10:160

x=sort([ybin(1:end),xbin(1:end)])

a=xbin;
b=ybin;

for i=1:26
    dx(i)=(x(i+1)-x(i));
end
dx=dx./10;


nrbands=[2 2 3 3 3 3 3 2 2 2]


clear dur_w avrg avrgbin dur avrgdur avrgspeed

v1=[1,5,9,12,17]
v2=[4,8,12,16,20]

for j=1:5
v=v1(j):v2(j)
n=1;
for i=1:10
    avrg(v,i)=0;
    dur(v,i)=0
    for ii=1:nrbands(i);
        dur_w(v,n)=duration(v,ii).*dx(n);
        avrg(v,i)=avrg(v,i)+dur_w(v,n).*v';
        dur(v,i)=dur(v,i)+dur_w(v,n)
        n=n+1
        

    end
    avrgbin(j,i)=sum(avrg(v,i));
    avrgdur(j,i)=sum(dur(v,i));
    avrgspeed(j,i)=avrgbin(j,i)/avrgdur(j,i);

end

end

%is avrg dur de benodigde duration?

plotavrg=zeros(5,11);plotavrg(1:5,1:10)=avrgbin(:,:)
plotspeed=zeros(5,11);plotspeed(1:5,1:10)=avrgspeed(:,:)
figure
subplot(2,2,1);pcolor(avrgdur)
subplot(2,2,2);pcolor(plotspeed)


%list

n=1;
for i=1:10
     for j=1:5
        climate_list(n,1)=(xbin(i+1)-xbin(i))/2+xbin(i);
        climate_list(n,2)=avrgspeed(j,i);
        climate_list(n,3)=avrgdur(j,i);
        n=n+1;
    end
end

%randomize list

for i=1:50
    randseq=randperm(50);
    climate_list_rand(i,:)=climate_list(randseq(i),:);
end    

%shift directions
climate_list_shift=climate_list_rand(:,1)-90

%fix windspeed -> u*
climate_list_rand_uster = climate_list_rand(:,2)./((1/0.4)*log(10/10^-3));

% cumsum duration
climate_list_rand2=cumsum(climate_list_rand(:,3));
climate_list_rand3(1,1)=0;
climate_list_rand3(2:50,1)=climate_list_rand2(1:49);

delete('winddir.dat');delete('windspeed.dat');delete('windduration.dat');

fid = fopen('winddir.dat','w+');
fprintf(fid, '%6.0f', climate_list_shift);
fclose(fid);

fid = fopen('windspeed.dat','w');
fprintf(fid, '%6.2f', climate_list_rand_uster);
fclose(fid);

fid = fopen('windduration.dat','w');
fprintf(fid, '%6.0f', climate_list_rand3);
fclose(fid);

%% avrg in bins

% y=1
% x=1
% for y=1:3
% 
% dirbin=x
% i=ybin(y):ybin(y+1)-1
% 
% (1*sum(duration(i,dirbin).*i')/(sum(duration(i,dirbin)))+...
% (sum(duration(i,dirbin+1).*i')/(sum(duration(i,dirbin+1)))*(xbin(x+1)/(dirbin+1*10)-1) ))/...
% (xbin(x+1)/(dirbin+1*10))    
% 
% end
% 
% 
% 
% 
% 
% avrg=duration(i,dirbin).*i'+duration(i,dirbin+1).*i'
% avrgclimate(y,x)=sum(avrg)/(sum(duration(i,1))+xbin(x+1)/(dirbin+1*10)*(sum(duration(i+1,1))))
% 
% avrg=duration(i,dirbin+1).*i'
% avrgclimate(y,x)=avrgclimate(y,x)+(sum(avrg)/sum(duration(i,1)))*xbin(x+1)/(dirbin+1*10)



%% plot climate

climate(:,37)=0; %dummycell
climate(22,:)=0 %dummycell

figure('units','normalized','outerposition',[0 0 1 1]);
subplot(12,1,1:7)
pcolor(climate);shading flat;axis tight
set(gca, 'TickDir', 'out', 'XTick', [1.5:2:370], 'YTick', [1.5:1:22],'XTickLabel',{[0:20:360]},'YTickLabel',{[0:1:20]},'fontsize',8,'FontName','arial')

%xlabel( 'Direction (deg)', 'FontSize', 10 ,'FontName','Arial');
ylabel( 'Wind Speed (m/s)', 'FontSize', 10, 'Rotation', 90 ,'FontName','Arial');
colorbar('YTickLabel',{'0 hours','10 hours','20 hours','30 hours','40 hours','50 hours'});
title('Windspeed vs direction', 'FontSize', 11 )


%% reduction

for i=1:18
    red_dir_climate(:,i)=climate(:,2*i-1)+climate(:,2*i);
end

for i=1:5
    try red_climate(i,:)=red_dir_climate(4*i-2,:)+red_dir_climate(4*i-1,:)+red_dir_climate(4*i-0,:)+red_dir_climate(4*i+1,:);end
    i
end

red_climate(:,19)=0 %dummycell
red_climate(6,:)=0 %dummycell

subplot(12,1,9:12)
pcolor(red_climate);shading flat;axis tight;

xlabel( 'Direction (deg)', 'FontSize', 10 ,'FontName','Arial');
ylabel( 'Wind Speed (m/s)', 'FontSize', 10, 'Rotation', 90 ,'FontName','Arial');
set(gca, 'TickDir', 'out', 'XTick', [1.5:1:370], 'YTick', [1.5:1:21],'XTickLabel',{[5:20:360]},'YTickLabel',{[2.5:4:20]},'fontsize',8,'FontName','arial')
colorbar('YTickLabel',{' 0 hours',' 50 hours','100 hours','150 hours','200 hours','250 hours','300 hours'});
title('Windspeed vs direction reduced', 'FontSize', 11 )


%%
% n=1;
% direction=5;
% v=2.5;
% 
% for i=1:10
%     v=2.5;
%     for j=1:5
%         climate_list(n,1)=direction;
%         climate_list(n,2)=v;
%         climate_list(n,3)=red_climate(j,i);
%         v=v+4;
%         n=n+1;
%     end
%     direction=direction+20;
%     
% end
% 
% %randomize list
% 
% for i=1:50
%     randseq=randperm(50);
%     climate_list_rand(i,:)=climate_list(randseq(i),:);
% end

%% save data to file

% fid = fopen('windfield.dat','w');
% 
% fprintf(fid, '%6.0f %6.0f %6.0f\r\n', climate_list_rand)
% 
% fclose(fid);


%% save

print -depsc -tiff -r0 plot_waveclimate.eps
print -r0 -dtiff plot_waveclimate.tif


       

