OPT = struct(...
    'maindir_decc', 'd:\cornelisdenheijer\sources\dune_erosion_coastline_curvature\simulations\00000000\xbeach\',...
    'regexp_select_decc', 'run(40).*?0\\xb_sum\.nc' ...
    );

reference_erosion = 162.0591;
%% curvature
[~, ~, files] = dirlisting(OPT.maindir_decc);

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

for iradii = 1:length(ncfiles)
    ncfile = ncfiles{iradii};
    
    result(iradii) = decc_analyse(ncfile,...
        'philimits', [-90 90],...
        'phi_resolution', 1,...
        's', false,...
        'zb', true,...
        'contour', -2.2);
end

% radii_decc = cellfun(@(x,y) repmat(x, size(y)), {result.radius}, {result.Ero}, 'UniformOutput', false);
radii_decc = {result.radius};
Ero = cellfun(@(x) nanmax([x(2:end-1); fliplr(x(2:end-1))]), {result.Ero},...
    'uniformoutput', false);
phi_decc = cellfun(@(x) x(2:end-1), {result.phi},...
    'uniformoutput', false);
ero_decc = cellfun(@(Ero) Ero/reference_erosion-1, Ero, 'uniformoutput', false);

%% plot
x = vertcat(phi_decc{:})';
y = vertcat(ero_decc{:})';
id = x(:,1) >= 0;
[x2 y2] = deal(50:-10:0, [0.2809 0.3011 0.2686 0.188 0.074 0]);
ph = plot(x(id,:),y(id,:),x2,y2, 'x:');

legtxt = [cellfun(@(x) sprintf('radius %i m', x), radii_decc, 'uniformoutput', false) {'straight coast'}];

%
hsv1 = rgb2hsv([0 0 1]);
hsv2 = rgb2hsv([1 0 1]);
h = interp1(0:1, [hsv1(1) hsv2(1)], linspace(0,1,size(x,2)));
colors = cellfun(@(h) hsv2rgb([h .8 hsv1(3)]), num2cell(h),...
    'uniformoutput', false);

set(ph,...
    {'color'}, [colors'; {[0 0 0]}],...
    {'displayname'}, legtxt')

legend(gca,...
    'location', 'south')

xlabel('Angle w.r.t. shore normal [degrees]')
ylabel('Portion additional erosion volume')

figureStyle(@styleThesis,...
    {'figure' 'Position' [1 1 15 10]});
tightenfigure
[fpath fname] = fileparts(which(mfilename));
fname = 'decc_dewa_additional_erosion';
figurepath = fpath;%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']))

%%
Res = 70:10:100;
% initialize U to zero
setpref('FORM', 'finalU', 0)

% %% reference calculation
% reffun = @(Res) FORM(...
%     'x2zFunction', @x2z_DUROS,...
%     'startU', @() getpref('FORM', 'finalU'),...
%     'stochast', [exampleStochastVar struct('Name', {'Angle'}, 'Distr', {@deterministic}, 'Params', {{0}})],...
%     'x2zVariables', {'Resistance', Res});
% resultFORM = cellfun(reffun,...
%     num2cell(Res));
% 
% P_f = cellfun(@(x) x.P_f, {resultFORM.Output});
% P_fround = 10.^(ceil(min(log10(P_f))):floor(max(log10(P_f))));
% 
% AFdes = sort(round(interp1(log10(P_f), Res, log10(P_fround))));
% 
% resultFORMdes = cellfun(reffun,...
%     num2cell(AFdes));
% P_fdes = cellfun(@(x) x.P_f, {resultFORMdes.Output});

%% calculation including curvature and wave angle
% initialize U to zero
for iradii = 1:length(radii_decc)
    AngleFcn = @(Angle) interp1(phi_decc{iradii}, ero_decc{iradii}, Angle);
    setpref('FORM', 'finalU', 0)
    anglefun = @(Res) FORM(...
        'x2zFunction', @x2z_DUROS,...
        'startU', @() getpref('FORM', 'finalU'),...
        'stochast', [exampleStochastVar struct('Name', {'Angle'}, 'Distr', {@unif_inv}, 'Params', {{10 40}})],...
        'x2zVariables', {'Resistance', Res, 'AngleFcn', AngleFcn});
    resultFORM_angle = cellfun(anglefun,...
        num2cell(Res));
    
    P_f_angle = cellfun(@(x) x.P_f, {resultFORM_angle.Output});
    P_fround_angle = 10.^(ceil(min(log10(P_f_angle))):floor(max(log10(P_f_angle))));
    
    AFdes_angle = sort(round(interp1(log10(P_f_angle), Res, log10(P_fround_angle))));
    
    resultFORMdes_angle = cellfun(anglefun,...
        num2cell(AFdes_angle));
    P_fdes_angle = cellfun(@(x) x.P_f, {resultFORMdes_angle.Output});
    
    save(sprintf('radius%05i.mat', radii_decc{iradii}), 'resultFORM_angle', 'P_f_angle', 'P_fround_angle', 'AFdes_angle', 'resultFORMdes_angle', 'P_fdes_angle')
end

%% create a plot of the probability of failure (=exceedance) as function of
%% the retreat distance

semilogy(...
    Res, P_f, '-ob',...
    AFdes, P_fdes, 'ob',...
    Res, P_f_angle, '-xr',...
    AFdes_angle, P_fdes_angle, 'xr');
xlabel('Retreat distance [m]')
ylabel('Probability of exceedance [y^{-1}]')

%%
prob_FORM_designpoint(resultFORMdes_angle, 'rowdelimiter', '\\\\\n', 'columndelimiter', '&')
% varnames = cellfun(@(x) x, {resultFORMdes_angle(end).Input.Name}, 'uniformoutput', false);
% data = [varnames' num2cell([resultFORMdes(end).Output.alpha' resultFORMdes_angle(end).Output.alpha'].^2 *100)]';
% fprintf('%-10s%5.1f%5.1f\n', data{:})