close all; clear all;

% // If done, read the model output  
xb = xb_getdimensions('J:\WP4\t156_2D\');
xbo = xb_read_output('J:\WP4\t156_2D\', 'vars', {'zb'});
[zb] = xb_get(xbo,'zb'); 

%xbo = xb_read_output('J:\WP4\t156_2D\', 'vars', {'zb', 'zs', 'u', 'v', 'H'});
%[u v zs zb H] = xb_get(xbo,'u','v','zs','zb','H'); 

%vmag = sqrt(u.^2+v.^2);

%% animation
%figure(1)
%for i = 1:xb.nt;
    %subplot(2,2,1);
    %pcolor(xb.x,xb.y,squeeze(H(i,:,:))'); shading interp; caxis([0 6]); colorbar; hold on; 
    %xlabel('x [m]'); ylabel('y [m]'); title('H [m]'); axis equal;
    
    %subplot(2,2,2);
    %pcolor(xb.x,xb.y,squeeze(vmag(i,:,:))'); shading interp; caxis([-1 1]); colorbar; hold on;
    %quiver(xb.x,xb.y,squeeze(u(i,:,:))',squeeze(v(i,:,:))','w'); hold off;
    %xlabel('x [m]'); ylabel('y [m]'); title('|v| [m/s]'); axis equal;
    
    %subplot(2,2,3);
    %pcolor(xb.x,xb.y,squeeze(zs(i,:,:))'); shading interp; caxis([-1 1]); colorbar;
    %xlabel('x [m]'); ylabel('y [m]'); title('z_s [m]'); axis equal;
    
   % subplot(2,2,4);
   % pcolor(xb.x,xb.y,squeeze(zb(i,:,:)-zb(1,:,:))'); shading interp; caxis([-0.5 0.5]); colorbar;
    %xlabel('x [m]'); ylabel('y [m]'); title('\Delta z_b [m]'); axis equal;   
    %pause(1); hold off;
%end

%% compare measured and simulated bed level changes

% rea measured data from files
xyzpre = load('pre_grd.dat');
Zpre = griddata(xyzpre(:,1),xyzpre(:,2),xyzpre(:,3),xb.x,xb.y);
xyzpost = load('post_grd.dat');
Zpost = griddata(xyzpost(:,1),xyzpost(:,2),xyzpost(:,3),xb.x,xb.y);

%figure(2);
%subplot(121); pcolor(xb.x,xb.y,Zpost-Zpre); shading flat; axis equal;
%caxis ([-1.5 1.5])
%subplot(122); pcolor(xb.x,xb.y,squeeze(zb(end,:,:)-zb(1,:,:))'); shading flat; axis equal;
%caxis ([-1.5 1.5])

%% compare transects

% 4 transect
%xyt = [9639422.95797 4658381.84679
       %9640416.93795 4658571.97036];
       
 %18 transect
 %xyt = [9639604.74298 4657706.40745
        %9640550.35806 4657886.52444];
  
 % 20 transect
 %xyt = [9639678.12383 4657623.01976
        %9640563.70017 4657793.13042];
    
    
xori  = 9639388 ;
yori  = 4657484;
alpha = 190.525787;

[xr yr] = xb_grid_rotate(xb.x-xori, xb.y-yori, alpha-360);

figure (2); pcolor(xb.xc,xb.yc,squeeze(zb(1,:,:))'); shading flat; axis equal; hold on;
        plot(xr,yr,'g-o');

dx = (max(max(xr))-min(min(xr)))/150;
dy = (yr(2,2)-yr(1,1))/150;
xt = min(min(xr)):dx:max(max(xr));
yt = yr(1,1):dy:yr(2,2);
zbs0 = interp2(xb.xc',xb.yc',squeeze(zb(1,:,:)),xt',yt');
zbsend = interp2(xb.xc',xb.yc',squeeze(zb(end,:,:)),xt',yt');
zbmend = griddata(xyzpost(:,1),xyzpost(:,2),xyzpost(:,3),xt',yt');

%n = 85;
figure(3);
subplot(221); pcolor(xb.x,xb.y,Zpost-Zpre); shading flat; axis equal; hold on; 
caxis ([-1.5 1.5])
%plot(xb.x(:,n),xb.y(:,n),'r','LineWidth',1.5); 
plot([xyt(1,1) xyt(2,1)],[xyt(1,2) xyt(2,2)],'r-*','LineWidth',1.5); axis image

subplot(222); pcolor(xb.x,xb.y,squeeze(zb(end,:,:)-zb(1,:,:))');shading flat; hold on;
caxis ([-1.5 1.5])
%plot(xb.x(:,n),xb.y(:,n),'r','LineWidth',1.5); 
plot([xyt(1,1) xyt(2,1)],[xyt(1,2) xyt(2,2)],'r-*','LineWidth',1.5); axis image

subplot(2,2,[3:4]);
plot(xt,zbs0,'k--','LineWidth',1.5);hold on;
title('PROFILE 4','fontsize',12,'fontweight','b');hold on;
xlabel('Cross-shore distance, m','fontsize',10,'fontweight','b');hold on;
ylabel('Depth, m','fontsize',10,'fontweight','b');hold on;
set(gca,'ylim',[-20 5],'xlim',[-1400 -180])
plot(xt,zbsend,'r-','LineWidth',1.5);
set(gca,'ylim',[-20 5],'xlim',[-1400 -180])
plot(xt,zbmend,'k-','LineWidth',1.5);
set(gca,'ylim',[-20 5],'xlim',[-1400 -180])
%save('2D_pr20_pre','xt','zbs0','-ASCII')
%save('2D_pr20_calc','xt','zbsend','-ASCII')
%save('2D_pr20_post','xt','zbmend','-ASCII')