close all; clear all;

dir = [pwd,'\'];

nam = [{'zb'};{'H'};{'zs'};{'zs0'};{'u'};];
%% dimensions
fid = fopen([dir,'dims.dat'],'r'); 
temp = fread(fid,[3,1],'double'); 
par.nt = temp(1);
par.nx = temp(2)+1;
par.ny = temp(3)+1;
fclose(fid);

%% read grid coordinates 
fid = fopen([dir,'xy.dat'],'r');
d.xw = fread(fid,[par.nx,par.ny],'double'); 
d.yw = fread(fid,[par.nx,par.ny],'double');
d.x = fread(fid,[par.nx,par.ny],'double'); 
d.y = fread(fid,[par.nx,par.ny],'double');
fclose(fid);

%% read XBeach output

for j = 1:length(nam)
    temp = zeros(par.nx,par.ny,par.nt);
    fid = fopen([dir,nam{j},'.dat'],'r');
    for i = 1:par.nt
        temp(:,:,i) = fread(fid,[par.nx,par.ny],'double');  % all data
    end
    fclose(fid);
    d.(nam{j}) = zeros(par.nx,par.nt);
    d.(nam{j}) = squeeze(temp(:,2,:));
end

%% movieplot
par.xstart=0;
par.zs0 =0;

ind = 2;

hfac = 8/1000/9.81;
fig = figure;
set(fig,'DoubleBuffer','on');
set(gca,'FontSize',12,'nextplot','replace','Visible','off');
for i = 1:par.nt
    % water surface elevation & short wave height
    plot(d.x(:,ind)+par.xstart,d.zs(:,i)-par.zs0,'b','LineWidth',1.5); hold on; 
    % flows
    % profiles
    plot(d.x(:,ind)+par.xstart,d.zb(:,1)-par.zs0,'k--'); hold on; 
    plot(d.x(:,ind)+par.xstart,d.zb(:,i)-par.zs0,'k','LineWidth',2);
    % figure settings
    title([num2str(i*2/60,'%1.1f'),' hours']); legend('water surface (including long waves)','short wave height','Eulerian flow velocity','initial profile','simulated profile',4);
    xlabel('x [m]'); ylabel('z_b [m] & \eta [m] & H_{rms} [m] & u_E [m/s]');
    axis([1000 2000 -10 10])
    pause(0.1); 
    pause(0.1); hold off;
end