rundir = fullfile(fileparts(which(mfilename)), '..\simulations\00000000\xbeach\run09\');

D = dir2(...
    'basepath', rundir,...
    'dir_incl', '.*270$',...
    'file_incl', 'xb_sum.nc$',...
    'no_dirs', true);

ncfiles = cellfun(@(p,f) [p f], {D.pathname}, {D.name},...
    'uniformoutput', false);

t = cellfun(@(x) nc_varget(x, 'globaltime'), ncfiles,...
    'uniformoutput', false);
[tmin tmax] = cellfun(@(x) deal(min(x), max(x)), t);
% assign tmax as the tend of the shortest simulation
tmax = min(tmax-tmin) + tmin;
for i = 1:length(t)
    if t{i}(end) > tmax(i)
        fprintf('WARNING: simulation only analysed upto t = %.0f (of %.0f)\n', tmax(i), t{i}(end))
    end
end
t = cellfun(@(x,y) x(x<=y), t, num2cell(tmax), 'uniformoutput', false);
tdim = cellfun(@length, t);

tm = cellfun(@(x) nc_varget(x, 'meantime'), ncfiles,...
    'uniformoutput', false);
tm = cellfun(@(x,y) x(x<=y), tm, num2cell(tmax), 'uniformoutput', false);
tmdim = cellfun(@length, tm);

for icalc = 1:length(ncfiles)
    ncfile = ncfiles{icalc};
    [x(icalc,:) zb0(icalc,:) zb1(icalc,:) ccg(icalc,:) Svsg(icalc,:) ve(icalc,:) tmpEro tmpSed zs(icalc,:)] = duwa_analyse(ncfile,...
        'tdim', tdim(icalc),...
        'tmdim', tmdim(icalc));
    Ero{icalc} = tmpEro;
    Sed{icalc} = tmpSed;
    Bal{icalc} = tmpSed+tmpEro;
end

[zwl ixwl] = unique(mean(zb0,1));
xwl = interp1(zwl, x(1,ixwl), 5) - x(1);

angledirs = dir(fullfile(fileparts(ncfile), '..', 'angle*'));
% anglestrings = {angledirs.name};
angles = 270 - cellfun(@str2double, regexp({angledirs.name}, '\d+', 'match'));
anglestrings = cellfun(@(x) sprintf('%2i\\circ', x), num2cell(angles),...
    'uniformoutput', false);
upstreambound = ceil(tand(angles)*xwl/20);
nplots = 5;

fh = figure;

subplot(nplots,1,1);
plot(x', zb0', ':')
hold on
plot(x', zb1')
ylabel('zb')
legend(anglestrings{:})

subplot(nplots,1,2);
plot(x', Svsg')
ylabel('Svsg')

subplot(nplots,1,3);
plot(x', ccg'*2650)
ylabel('ccg')

subplot(nplots,1,4);
plot(x', ve')
ylabel('ve')

subplot(nplots,1,5);
plot(x', zs')
ylabel('zs')

linkaxes(findobj(fh, 'type', 'axes', '-and', '-not', 'tag', 'legend'), 'x')

%%
c = 'rgbcymk';
figure;
% marg = 0:.05:.4;
for j = 1:3
    for i = 1:length(Ero)
        subplot(3,1,j)
        xp{i} = linspace(0, 1, length(Ero{i}));
%         plot(xp{i}, Ero{i}, c(i))
switch j
    case 1
        plot(Ero{i}, c(i))
    case 2
        plot(Sed{i}, c(i))
    case 3
        plot(Bal{i}, c(i))
end

        hold on
        vline([upstreambound(i) length(Ero{i}) - upstreambound(i)], c(i))
%         mean_id = xp{i} > marg(imarg) & xp{i} < 1-marg(imarg);
%         mean_Ero(imarg,i) = mean(Ero{i}(mean_id));
    end
end
legend(anglestrings{:})

%%
A = cellfun(@(x,y) mean(x(max(y,1):end-y)), Ero, num2cell(upstreambound));
save(fullfile(fileparts(which(mfilename)), 'wave_angle_erosion.mat'), 'A', 'angles')