%%
OPT = struct(...
    'maindir_dewa', 'D:\cornelisdenheijer\sources\dune_erosion_wave_angle\simulations\00000000\xbeach\',...
    'regexp_select_dewa', 'run(09).*?0\\xb_sum\.nc' ...
    );

%% wave angle
[~, ~, files] = dirlisting(OPT.maindir_dewa);

ncfiles = files(~cellfun(@isempty, regexp(files, OPT.regexp_select_dewa)));

zb1 = cellfun(@(ncfile) nc_varget(ncfile, 'zb', [0 0 0], [2 -1 -1], [90 1 1]), ncfiles, 'uniformoutput', false);
[tdim ydim xdim] = cellfun(@(z) size(z), zb1);
zb0 = cellfun(@(z) repmat(z(1,:,:), [size(z,1) 1 1]), zb1, 'uniformoutput', false);
[x y t] = cellfun(@(ncfile) deal(nc_varget(ncfile, 'globalx'), nc_varget(ncfile, 'globaly'), nc_varget(ncfile, 'globaltime')), ncfiles, 'uniformoutput', false); 
t = cellfun(@(t) t - min(t(:)), t, 'uniformoutput', false);
g = cellfun(@xb_stagger, x, y);
A = cellfun(@(x,y,tdim) repmat(shiftdim(x'./y', -1), [tdim 1 1]), {g.dsdnz}, {g.dnz}, num2cell(tdim), 'uniformoutput', false);
Ero = cellfun(@(zb0,zb1,A) sum((max(zb1, 5) - max(zb0, 5)) .* A,3), zb0, zb1, A, 'uniformoutput', false);

idx_zonebounds = cellfun(@(y) interp1(y(:,1)/y(end,1), 1:size(y,1), [.3 .7], 'nearest'), y, 'uniformoutput', false);

Ero = cellfun(@(Ero,idx) -Ero(end,idx(1):idx(2)), Ero, idx_zonebounds, 'uniformoutput', false);

m = regexp(ncfiles, '(?<=angle)\d+', 'match');
angles = cellfun(@(x) str2double(x), [m{:}]);
phi_dewa = 270-angles;
ero_dewa = cellfun(@mean, Ero);

reference_erosion = mean(ero_dewa(phi_dewa==0));
