OPT = struct(...
    'maindir', fullfile(fileparts(which(mfilename)), '..', 'simulations'),...
    'regexp_select', 'run(09).*xb_sum.nc',...
    'regexp_legend', 'run\d{2}.*?angle\d{3}(.*?rev\d{4}\D?)?',...
    'strrep_legend', {{{'\' '; '} {'s0' 's '} {'00' '00 m'}}},...
    'offset', 0,...
    'n_contours', 5,...
    'save', true,...
    'style', @styleThesis ...
    );

% OPT = setproperty(OPT, varargin);

%%
[~, ~, files] = dirlisting(OPT.maindir);

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

legtxt = cellfun(@(x) regexp(x, OPT.regexp_legend, 'match'), ncfiles);
for ilegrep = 1:length(OPT.strrep_legend)
    legtxt = cellfun(@(x) strrep(x, OPT.strrep_legend{ilegrep}{:}), legtxt,...
        'UniformOutput', false);
end
%%
if isscalar(ncfiles)
    fprintf('%s\n', ncfiles{:});
else
    fparts = cellfun(@(x) textscan(x, '%s', 'delimiter', repmat(filesep, 1, ispc+1)), ncfiles, 'UniformOutput', false);
    fparts = [fparts{:}];
    for i = 1:min(cellfun(@length, fparts))-1
        if ~isscalar(unique(cellfun(@(x) x{i}, fparts, 'UniformOutput', false)))
            fprintf('%s\n', fullfile(fparts{1}{1:i-1}));
            for j = 1:length(ncfiles)
                fprintf('\t%s\n', fullfile(' ', fparts{j}{i:end}));
            end
            break
        end
    end
end

%%
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,y) 270-repmat(str2double(x), 1, length(y)), [m{:}], Ero, 'uniformoutput', false);
angles_mean = cellfun(@mean, angles);
angles = [angles{:}];
Ero_std = cellfun(@std, Ero);
Ero_mean = cellfun(@mean, Ero);
Ero = [Ero{:}];

if OPT.save
    tag = '$HeadURL$ $Revision$';
    save('wave_angle_erosion_xbeach.mat', 'Ero', 'angles', 'Ero_mean', 'angles_mean', 'Ero_std', 'tag')
end