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});

save('reference.mat', 'resultFORMdes', 'P_fdes', 'P_fround', 'AFdes')


%% calculation including wave angle
AngleFcn = @(Angle) interp1(50:-10:0, [0.2809 0.3011 0.2686 0.188 0.074 0], Angle);
% initialize U to zero
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('straightcoast.mat', 'resultFORMdes_angle', 'P_fdes_angle', 'P_fround_angle', 'AFdes_angle')

%% 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{:})