%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIS CODE COMPUTES THE ECs FOR POINT_005284 AND POINT_018830 ACCORDING
% TO VARYING LEVELS OF t AND Tr
% De Leo F. - 10/2020
% francesco.deleo@edu.unige.it / fdeleo@calpoly.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clearvars 
close all
clc

%% USER EDITABLE PARAMETERS
iloc = 1;   % 1/2 -> Point_005284 / Point_010830
t    = 40;  % reference time t [years]
Tr   = 100; % return period [years]

%% DEFINE GEV AND N-GEV PARAMETERS
% stationary parameters
gev  = [4.77 0.91 -0.18;
        5.92 0.99 -0.25];
locST = gev(iloc,1);
scST  = gev(iloc,2);
shST  = gev(iloc,3);

% non-stationary parameters
ngev = [4.39  0.020	 0.93 -0.25;
        6.71  -0.034 1.02 -0.40];
locNS  = ngev(iloc,1);
loctNS = ngev(iloc,2);
scNS   = ngev(iloc,3);
shNS   = ngev(iloc,4);

%% DEFINE a AND b coefficients FOR IFORM MODELLING
a = [0.56 1.18 0.26;
     0.20 1.53 0.20];
    
b = [0.017 0.32 -0.39;
     0.045 0.17 -0.52];

a1 = a(iloc,1);
a2 = a(iloc,2);
a3 = a(iloc,3);

b1 = b(iloc,1);
b2 = b(iloc,2);
b3 = b(iloc,3);

%% I-FORM ENVIRONMENTAL CONTOUR COMPUTATION    
F = 1 - 1./(Tr);          % evaluate F(Tr)

% analysis in the Gaussian space
u = norminv(F);
u1 = -u:0.0001:u;                               % Hs
u2 = [sqrt(u.^2-u1.^2) -sqrt(u.^2-u1.^2)];      % Tp

fu1 = normcdf([u1 flip(u1)]);

% STATIONARY ANALYSIS
uH = gevinv(fu1,shST,scST,locST);
u2mu = a1 + a2*uH.^(a3);
u2sd = b1 + b2.*exp(b3*uH);
uT = exp(u2mu + (u2.*u2sd));

% NON-STATIONARY ANALYSIS
uHt = gevinv(fu1,shNS,scNS,locNS+loctNS*t);
u2mu = a1 + a2*uHt.^(a3);
u2sd = b1 + b2.*exp(b3*uHt);
uTt = exp(u2mu + (u2.*u2sd));

cplt = [1 0 0;0 0 1];

figure
plot(uH,uT,'-k');hold on
plot(uHt,uTt,'-','Color',cplt(iloc,:))
xlabel('H_s [m]')
ylabel('T_p [s]')
set(gca,'FontSize',12)


