function Y = uq_Ex08_BermBreakwater_constraint(X)
%% 1. INPOUT
% 1.1 Optim var
B = X(:,1);
hb = X(:,2);
Ns = X(:,3);
hf = X(:,4);
% 1.2 Enviromental var
Hs = X(:,5);
Tp = X(:,6);
ro_s = X(:,7);
ro_w = X(:,8);
d = X(:,9);
N = X(:,10);
fg = X(:,11);
DD = X(:,12);
CR = X(:,13);
CO1 = X(:,14);
CO2 = X(:,15);

%% 2. LIMIT STATE FUNCTIONS
% 2.1 REC

Hsd = 5.03;
delta=(ro_s./ro_w)-1;
deltad=2.45/1.025-1.;
Dnd = (Hsd./(deltad.*Ns));
Dn = DD.*Dnd; %apply D uncertaintity
%Dn = (Hsd./(deltad.*Ns)); %do not apply D uncertaintity
HSTp = (Hs./(delta.*Dn)).*sqrt(Tp.*sqrt(9.81./Dn));

Y(:,1) = (B) - CR.*(((0.02067.*(HSTp.^1.971).*((d./Dn).^-0.1087).*log((N.*d)./(3000.*Dn)))+(21133.*(HSTp.^2.094).*(hb./Hs).*(d./Dn).*(log(fg).^(7.187.^fg)))-(81.74.*(10.^-5).*(1.378.^(hb./Hs)).*(HSTp.^3).*(N./3000).*((hb./Hs).^2))-1.245).*Dn);

% 2.2 Overtopping
C = 17.8;
hc = (2 .* Dn);
Rc = (C+hc+hf-d);
% 2.2.1 gama_bb
Som = (2.*pi.*Hs./(9.81.*((Tp./1.25).^2)));
gama_bb = (Ns > 2.5) .* (0.7 - 8.2 * Som) + (Ns <= 2.5) .* (0.68 - 4.1 * Som - 0.05 * B ./ Hs);

Y(:,2) = 0.01 - (CO1.*exp(-((CO2.*Rc)./(Hs.*gama_bb.*1)).^1.3)).*(sqrt(9.81.*(Hs.^3)));

end