close all; clear all; clc;

% hourly winddata 2010 HvH(?)at 10 meter height
wdata=csvread('wind2010v2.txt');
dirt=wdata(:,4)';
uw=wdata(:,5)'/10;

% make wind rose uw
wind_rose(dirt,uw);

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

ustart=uw/((1/0.4)*log(10/10^-3)); % check this
ustar_crit=0.22; 

% make wind rose ustar and determine probabilities
[handles,reduced_data] = wind_rose(dirt,ustart,'n',36);

p = reduced_data/sum(sum(reduced_data));

[n,id] = size(reduced_data);
[dir, ustar] = meshgrid([0:1:n-1]*360/n, [0.5:1:id]*0.1);

% only onshore wind directions and wind intensities larger than ustar
cl_orientation = 40 % degrees Nautical with respect to north 
[j,i] = find(dir>cl_orientation-20 & dir<200 +cl_orientation);
[j2,i2] = find(ustar<ustar_crit);

p2 = p;
p2(i,j) = NaN;
p2(i2,j2) = NaN;

figure;
subplot(211);
pcolor(dir,ustar,p'); colorbar; caxis([0 0.02]);
xlabel('Direction [deg]'); ylabel('ustar [m/s]');
title('probability of wind conditions')
subplot(212);
pcolor(dir,ustar,p2'); colorbar; caxis([0 0.02]);
xlabel('Direction [deg]'); ylabel('ustar [m/s]');
title('probability of wind conditions that will cause transport')

%% now reduce wind conditions based on shear stress

addpath D:\documents\TuDelft\TU_CD2\2012\week4\UNIBESTtools\

tau = rhoa*ustar.^2;
tauw = tau.*p2';

nr_dir_bins = 8 % 15 degrees per bin
nr_ustar_bins = 6
% Number of grouped wave height bins; the grouping is based on an approx. equal weighted contribution

% set dir bins 
tauw(isnan(tauw))=0;
tauwD = cumsum(sum(tauw,1));
trh  = [0:1:nr_dir_bins]*tauwD(end)/nr_dir_bins;
for i = 2:nr_dir_bins+1
    dirbin{i-1} = find(tauwD <= trh(i) & tauwD>trh(i-1));
end

ddir = dir(1,2)-dir(1,1);
dustar = ustar(2,1)-ustar(1,1);
% and now reduce as in CDII
[Ustar_Dir_P] = wind_reduction(tauw,p',ustar(:,1)',dir(1,:),dirbin,nr_dir_bins,nr_ustar_bins,ddir,dustar); 

rn = randperm(length(Ustar_Dir_P));
Durt = 0;
for i = 1:length(rn)
    Durt = Durt + ceil(Ustar_Dir_P(rn(i),3)*365*24);
    Duration(i) = Durt;
    Directions(i) = 270-Ustar_Dir_P(rn(i),2);
    if Directions(i) > 180
        Directions(i) = Directions(i)-360
    elseif Directions(i) < -180
        Directions(i) = 
    Windspeeds(i) = Ustar_Dir_P(rn(i),1);
end
    
% save to files
save('windspeed.dat','Windspeeds','-ascii');
save('winddir.dat','Directions','-ascii');
save('windduration.dat','Duration','-ascii');