%%
close all
id = desb_illustration_id;

%% duros+ results
durosnc = fullfile(fileparts(which(mfilename)), '..', 'simulations', sprintf('%08i', id), 'durosplus', 'run00', sprintf('%08i_pp.nc', id));

x = nc_varget(durosnc, 'cross_shore');
z = nc_varget(durosnc, 'altitude');

%%
ph = plot(x, squeeze(z(43,:,:)));
set(ph(1), 'DisplayName', 'Initial profile')
set(ph(2), 'DisplayName', 'DUROS+ erosion profile')

set(gca,...
    'XLim', [600 1400],...
    'YLim', [-6 12])

ise = zz(1,:) == zz(2,:);
x_z = x(...
    ise + ise([1 1:end-1]) + ise([2:end end]) == 2 |... % most landward and most seaward boundary
    ise + ise([1 1:end-1]) + ise([2:end end]) == 1 & ise ==1); % intermediate boundary

zh = vline(x_z);
set(zh, 'DisplayName', 'Zone boundaries')
set(zh(end), 'HandleVisibility', 'on')

title(sprintf('transect %i', id))
xlabel('cross-shore coordinate [m]')
ylabel('Altitude [m] w.r.t. NAP')

x_z2 = sort([x_z' get(gca, 'xlim')]);
xl = mean(x_z2([1:end-1; 2:end]));
yl = min(get(gca, 'ylim'));
for ix = 1:length(xl)
    text(xl(ix), yl, num2str(ix),...
        'verticalalignment', 'bottom')
end

legend show

%%
fh = figureStyle(@styleThesis,...
    {'figure' 'Position', [1 1 14 7]});
tightenfigure

[fpath fname] = fileparts(which(mfilename));
figurepath = fullfile(fpath, '..', 'figures');
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']))