close all; clear all;

%% input wind
% in sand engine simulation this approach wa snot followed?! Why? --> What
% is Joost doing?
wDir = load('winddir.dat');
wDur = load('windduration.dat');
wSp  = load('windspeed.dat');

figure;
subplot(211);
plot(wDur,wSp,'g-','LineWidth',2);
subplot(212);
plot(wDur,wDir,'g-','LineWidth',2);

%% bathymetry
zb = load('h_init_ext.dat'); % initial bathyemtry
zveg = load('h_vegh_ext.dat'); % height vegetation
zvegpos = load('h_vegpos_ext.dat'); % what is this?

dx = 5;
[x,y] = meshgrid([0:1:size(zb,1)-1]*dx,[0:1:size(zb,2)-1]*dx);

figure; surf(x,y,zb); shading flat; hold on;  
        colorbar; caxis([-5 15]);
        xlabel('x [m]'); ylabel('y [m]');
        
figure; plot(x(1,:),zb(1,:)+zveg(1,:),'g-','LineWidth',1.5); hold on;
        plot(x(1,:),zb(1,:)+zvegpos(1,:),'g--','LineWidth',1.5);
        plot(x(1,:),zb(1,:),'k-','LineWidth',1.5);

%% output
data = load('data1\time.dat');
time = data(:,2)/3600;
tout = 1;
time_out = tout:tout:time(end);


% estimate area that is always moprhologically updated
nm = floor(min(size(x))/sqrt(2)/2);
% find centre of model domain
m = ceil(size(x,1)/2);
n = ceil(size(x,2)/2);
% find indices for analysis
indm = m-nm:m+nm-1;
indn = n-nm:n+nm-1;

% read outputfiles --> Check what is in other output files
h1d =zeros(length(time_out),length(x));
h2d =zeros(length(time_out),length(x),length(y));
h1ne = h1d;
veg2d = h2d;
UX = h2d;
UY = h2d;
dhdt2d = h2d;
dhdt2dcum = h2d;

for i = 1:length(time_out);
    h = load(['data1\h_rot.',num2str(time_out(i),'%05.0f'),'.dat']);
    hne = load(['data1\h_nonerod_rot.',num2str(time_out(i),'%05.0f'),'.dat']);
    veg = load(['data1\veget_rot.',num2str(time_out(i),'%05.0f'),'.dat']);
    % u_x = load(['data1\shear_rot_x.',num2str(time_out(i),'%05.0f'),'.dat']);
    % u_y = load(['data1\shear_rot_y.',num2str(time_out(i),'%05.0f'),'.dat']);
    % dhdt= load(['data1\dhdt_rot.',num2str(time_out(i),'%05.0f'),'.dat']);
    % UX(i,:,:) = u_x(indm,indn);
    % UY(i,:,:) = u_y(indm,indn);
    % dhdt2d(i,:,:) = dhdt(indm,indn);
    h2d(i,:,:) = h;
    dhdt2dcum(i,:,:) = squeeze(h2d(i,:,:))-zb;
    h1ne(i,:) = hne(91,:);
    veg1d(i,:) = veg(91,:);
    h1d(i,:) = h(91,:);
end

% [t,xt] = meshgrid(time_out,x(1,:));
% figure; surf(t,xt,h1d'); shading flat;

view = 0;
if view == 1;
    indp = indm
else
    indp = 1:length(h);
end

figure; 
for i = 1:size(UX,1)
    surf(x(indp,indp),y(indp,indp),squeeze(h2d(i,indp,indp))); shading flat; colorbar;  hold on; %caxis([-2 2]);
    pause; hold off;
end


% figure; 
% Umag = sqrt(UX.^2+UY.^2);
% for i = 1:size(UX,1)
%     %     plot(x,y,'b.'); hold on;
%     pcolor(x(indm,indn),y(indm,indn),squeeze(Umag(i,:,:))); shading flat; colorbar; hold on;
%     quiver(x(indm,indn),y(indm,indn),squeeze(UX(i,:,:)),squeeze(UY(i,:,:)),'w');
%     pause; hold off; 
% end

figure;
for i = 1:length(time_out);
    plot(x(1,indp),zb(91,indp),'k-','LineWidth',1.5); hold on;
    plot(x(1,indp),veg1d(i,indp),'g-','LineWidth',1.5);
    plot(x(1,indp),h1d(i,indp),'r-','LineWidth',1.5);
    plot(x(1,indp),h1d(i,indp),'r-','LineWidth',1.5);
    plot(x(1,indp),h1ne(i,indp),'k-o','LineWidth',1.5);
    % surf(x,y,squeeze(h2d(i,:,:))); shading flat; colorbar;
    title(['bed level after ',num2str(time_out(i)),' hours']);
    xlabel('x [m]'); zlabel('zb [m]');
    %Ylim([-5 15]);
    pause(0.1); hold off;  
end

figure; surf(x,y,squeeze(h)); shading flat; colorbar; view(2);
figure; surf(x,y,squeeze(hne)); shading flat; colorbar; view(2);
        