%%plot Sue results
%Gundula Winter
%2013-02-14

% This script requires the Open Earth Toolbox

%% clean up
clc
clear
close all

%% input
xbrun='Sue_01_nonh_gladijs';

xbpath=['..\Sue_experiment\',xbrun,'\'];

%% read results

output=xb_read_netcdf([xbpath,'xboutput.nc'],'vars',{'DIMS','zs','zb'});

xb_peel(output);
xb_peel(DIMS);

%% load experimental data
load('Sue_SG2_IS4\eta.mat')
eta_exp=eta/1000;
load('Sue_SG2_IS4\time.mat')
time_exp=time;
load('Sue_SG2_IS4\x.mat');
x_exp=x/1000;

%% plot

% landslide block
% t=0:0.01:6;
a=1*ones(1,length(t'));
v=zeros(1,length(t'));
pos=zeros(1,length(t'));
v(t<=1.2)=(a(t<=1.2)).*(t(t<=1.2))';
pos(t<=1.2)=0.5*(a(t<=1.2)).*(t(t<=1.2))'.^2;
a(t>1.2)=-t(t>1.2);
v(t>1.2)=v(t==1.2)+a(t>1.2).*(t(t>1.2)-1.2)';
pos(t>1.2)=pos(t==1.2)+v(t==1.2).*(t(t>1.2)-1.2)'+0.5*(a(t>1.2)).*(t(t>1.2)-1.2)'.^2;
a(t>1.31)=-0.25;
v(t>1.31)=v(t==1.31)+a(t>1.31).*(t(t>1.31)-1.31)';
pos(t>1.31)=pos(t==1.31)+v(t==1.31).*(t(t>1.31)-1.31)'+0.5*(a(t>1.31)).*(t(t>1.31)-1.31)'.^2;
ind_pos=find((v>0));
pos(v<0)=pos(ind_pos(end));
pos=pos+0.4665;

% 2D contour plot of water surface profile time history
fc_exp=figure;
[xp_exp,yp_exp]=meshgrid(x_exp,time_exp);
% contourf(xp,yp,squeeze(zs(:,1,:)),-0.03:0.001:0.03);%shading flat;
pcolor(xp_exp,yp_exp,(eta_exp(:,:)));shading flat;
hold on
colorbar;caxis([-0.025 0.015])
xlabel('x (m)')
ylabel('t (s)')
xlim([0 10])
ylim([0 7])
scatter(pos,t,'k')

fc=figure;
[xp,yp]=meshgrid(globalx_DATA,t);
% contourf(xp,yp,squeeze(zs(:,1,:)),-0.03:0.001:0.03);%shading flat;
pcolor(xp,yp,squeeze(zs(:,1,:)));shading flat;
hold on
colorbar;caxis([-0.025 0.015])
xlabel('x (m)')
ylabel('t (s)')
xlim([0 10])
ylim([0 7])
scatter(pos,t,'k')

% Water surface time history
ind_x=[find(globalx_DATA==0.5) find(globalx_DATA==1.5) find(globalx_DATA==2.5) find(globalx_DATA==3.5) find(globalx_DATA==4.5) find(globalx_DATA==5.5)];
ind_t=[find(abs(t-0.6)==min(abs(t-0.6))) find(abs(t-1.6)==min(abs(t-1.6))) find(abs(t-2.6)==min(abs(t-2.6)))...
    find(abs(t-3.6)==min(abs(t-3.6))) find(abs(t-4.6)==min(abs(t-4.6))) find(abs(t-5.6)==min(abs(t-5.6)))];
for i=1:6
    temp_ind=find(abs(x_exp-i+0.5)<0.0005);
    ind_x_exp(i)=temp_ind(1);
end
for i=1:6
    temp_ind=find(abs(time_exp-i+0.4)<0.0005);
    ind_t_exp(i)=temp_ind(1);
end


% Time history and profile
ft=figure;
for k=1:length(ind_x)
    subplot(length(ind_x),1,k)
    plot(t,squeeze(zs(:,1,ind_x(k))),'-b')    
    hold on
    plot(time_exp,eta_exp(:,ind_x_exp(k)),'-r')    
    title(['x = ',num2str(globalx_DATA(ind_x(k)),'%2.1f'),' m'])
    ylabel('zs [m]')
    xlabel('time [s]')
    ylim([-0.04 0.04])
    grid on
    if k==1
        legend('XBeach','experiment','location','east')
    end
end

fx=figure;
for k=1:length(ind_t)
    subplot(length(ind_t),1,k)
    plot(globalx_DATA,squeeze(zs(ind_t(k),1,:)))
    hold on
    plot(x_exp,eta_exp(ind_t_exp(k),:),'-r')   
    scatter(pos(ind_x(k)),0,50,'or','filled')
    title(['t = ',num2str(t(ind_t(k)),'%2.1f'),' m'])
    ylabel('zs [m]')
    xlabel('x [m]')
    ylim([-0.04 0.04])
    xlim([0 10])
    grid on
    if k==1
        legend('XBeach','experiment','location','east')
    end
end

set(ft,'paperunits','centimeters', 'papersize',[15 25],'paperposition',[0, 0, 15, 25])
set(fx,'paperunits','centimeters', 'papersize',[15 25],'paperposition',[0, 0, 15, 25])

print(fc,'-dpng',[xbpath,'\characteristic_curve_model.png']);
print(fc_exp,'-dpng',[xbpath,'\characteristic_curve_experiment.png']);
print(ft,'-dpng',[xbpath,'\time_history.png']);
print(fx,'-dpng',[xbpath,'\profile.png']);


