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

D = dir2(...
    'basepath', rundir,...
    'dir_incl', '.*',...
    '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);


result = dewa_analyse(ncfiles);%,...
%     'tdim', tdim(icalc),...
%     'tmdim', tmdim(icalc));
x = cell2mat({result.globalx}');
zb0 = cell2mat(cellfun(@(zb) zb(1,:), {result.zb}, 'uniformoutput', false)');
zb1 = cell2mat(cellfun(@(zb) zb(end,:), {result.zb}, 'uniformoutput', false)');
Svsg = cell2mat({result.Svsg_mean}');
Susg = cell2mat({result.Susg_mean}');
ve = cell2mat({result.ve_mean}');
ue = cell2mat({result.ue_mean}');
ccg = cell2mat({result.ccg_mean}');
zs = cell2mat({result.zs_mean}');

Ero = sum((max(zb1, 5) - max(zb0, 5)) .* diff(x(:,[1:end end]),1,2), 2);

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

angledirs = dir(fullfile(fileparts(ncfiles{1}), '..', '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);
tmp =  regexp(rundir, 'run\d{2}', 'match');
runname = tmp{1};
keys = fieldnames(result);
save([runname '.mat'], 'result', 'keys', 'angles')
nplots = 3;
%%
fh = figure;

xlim = [-400 max(x(:))];

xid = x>=min(xlim) & x<=max(xlim);

subplot(nplots,1,1);
plot(x', zb1')
hold on
ph = plot(x', zb0', ':');
set(ph, 'handlevisibility', 'off')
ylabel('Altitude [m]')
ylimitszb = minmax(zb0(xid)');
ylimitszb = mean(ylimitszb) + (diff(ylimitszb)/2  * [-1.1 1.1]);
set(gca, 'ylim', ylimitszb)
legend(anglestrings{:})

subplot(nplots,1,2);
plot(x', Svsg')
ylabel('LS transport [m^3/m/s]')
ylimitsSvsg = minmax(Svsg(:)');
ylimitsSusg = minmax(Susg(:)');
ylimits = max(cellfun(@diff, {ylimitsSvsg ylimitsSusg}))/2 * [-1.1 1.1];
set(gca, 'Ylim', mean(ylimitsSvsg) + ylimits)

subplot(nplots,1,3);
plot(x', Susg')
ylabel('CS transport [m^3/m/s]')
xlabel('Cross-shore coordinate [m]')
set(gca, 'Ylim', mean(ylimitsSusg) + ylimits)

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

figureStyle(@stylePresentation,...
    {'figure', 'Position', [1 1 18 21]});
[fpath fname] = fileparts(which(mfilename));
fname = 'dewa_plot_transport';
figurepath = fullfile(fpath);
set(gcf, 'tag', '$HeadURL$ $Revision$')
print(gcf, '-depsc2', fullfile(figurepath, [fname '.eps']))
if ispc
    saveas(gcf, fullfile(figurepath, [fname '.emf']))
end
saveas(gcf, fullfile(figurepath, [fname '.fig']))

%%
fh(2) = figure;
nplots = 3;

subplot(nplots,1,1);
plot(x', ccg'*2650)
ylabel('Concentration [kg/m^3]')
ylimitsccg = minmax(ccg(:)'*2650);
ylimitsccg = mean(ylimitsccg) + (diff(ylimitsccg)/2  * [-1.1 1.1]);
set(gca, 'ylim', ylimitsccg)


subplot(nplots,1,2);
plot(x', ve')
ylabel('LS velocity [m/s]')
ylimitsve = minmax(ve(:)');
ylimitsue = minmax(ue(:)');
ylimits = max(cellfun(@diff, {ylimitsve ylimitsue}))/2 * [-1.1 1.1];
set(gca, 'Ylim', mean(ylimitsve) + ylimits)

subplot(nplots,1,3);
plot(x', ue')
xlabel('Cross-shore coordinate [m]')
ylabel('CS velocity [m/s]')
set(gca, 'Ylim', mean(ylimitsue) + ylimits)

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

linkaxes(findobj(fh(2), 'type', 'axes', '-and', '-not', 'tag', 'legend'), 'x')
set(gca, 'xlim', xlim)

figureStyle(@stylePresentation,...
    {'figure', 'Position', [1 1 18 21]});
[fpath fname] = fileparts(which(mfilename));
fname = 'dewa_plot_concentration_velocity';
figurepath = fullfile(fpath);
set(gcf, 'tag', '$HeadURL$ $Revision$')
print(gcf, '-depsc2', fullfile(figurepath, [fname '.eps']))
if ispc
    saveas(gcf, fullfile(figurepath, [fname '.emf']))
end
saveas(gcf, fullfile(figurepath, [fname '.fig']))


%%
% 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'), 'Ero', 'angles')