Calculation of Vertical Wall Breakwaters using Goda formulation (1984)

Hi there!
a few months ago I got an assignment from my boss to calculate the wave force infront of a vertical wall breakwater using Goda formulation. i found it quite 'a work' if you have to calculate it one by one using a calculator so that i'm (again) doing it using matlab as always, the case that i use is taken from problem example in the book itself (Example 4.1 Goda, 1984). here it is:

Calculate the wave pressure, uplift pressure, and their moment produced by wave of the following characteristicsincident on the upright section of the vertical wall breakwater shown in Fig below:

Wave : H’0 = 6.3 [m] , T1/3= 11.4  [s] , b=15 [deg]

Tide level : W.L = + 0.6 [m]

Sea bottom slope tanθ = 1/100

And this is my matlab script to calculate the : wave pressure, uplift pressure and their moments

clc; clear all; close all; 

%  Wave force calculation according to GODA, 1985
%  mei.ramdhani@gmail.com
%% UPDATE
% 27/10/15 : Iteration for calculating alfa 1 and 1/cosh(2*pi()*h/L) added
% tobe Updated
% 25/10/15 : tobe added --> Newton rapshon iteration for calculation L0

%% INPUT
% depth definition
h01= 9.5; %from SWL to the seabed
h02= 6.5; %from SWL to the end of the vertical wall  or h'
h03= 5 ;  %from SWL to the bagian kaki dari vertical wall atau d
h04= 4;   %from crest of BW to the SWL or hc

% Waves
H_0     =6.3;  %[m] Significant wave height in offshore (H'0)
Tsig    =11.4; %[s] Significant wave period in offshore (T1/3)
beta    =15;   %[deg] angle between the direction of wave approach and a line normal to BW

% Tide Level
LWL=0 ;  % [m] Lowest water level
HWL=0.6; % [m] Highest water level
MSL=0.3; % [m] Mean sea level

% Sea Bottom Slope
tan_teta=1/100;  

%Width of The bottom of the upright section 
B=15 ; %[m]
H1_3    =5.8  ;  % H1/3 taken from mike SW
Hmax    =8    ;  % eq. 3.26 (Goda, 1984) or from wave model simulation

%% LOADING GODA CURVE
g10=load('goda10.dat'); % calculation diagram for parameter alpha1 (Goda, 1984)
g11=load('goda11.dat');
g20=load('goda20.dat');
g21=load('goda21.dat');

%% WATER DEPTH AND CREST ELEVATION
h = h01+HWL; %[m] h
h1= h02+HWL; %[m] h'
d = h03+HWL; %[m] d
hc= h04-HWL; %[m] hcrest

%% CALCULATION OF WAVE LENGTH AND WAVE HEIGHT
L0      =202.7 % [m] calculated using newton rapshon (later)
H1L1    =H_0/L0 ; %H'0/L0
hL0     =h/L0   ; %h/L0
hb      =h+(h03*H1_3/100);

%% COEFICIENT FOR WAVE PRESSURE
the coefficient of alpha 1 is taken from this curve:

Calculation diagram for alpha 1

Digitzed the diagram save it to g10 and g11 file..

%alfa1
temp1=h/L0;
if temp1 <= 0.1, alfa1 = polyval(g10,temp1),...
else,alfa1 = polyval(g11,temp1);...  
end;

%alfa2
tema1=((hb-d)/(3*hb)*(Hmax/d)^2);
tema2=((2*d)/Hmax);
alfa2=min([tema1,tema2]);

and so does for the par1 parameter,  i get it from this curve

Calculation diagram for the factor of 1/cosh(2*pi*h/L)

if temp1 <= 0.15, par1= polyval(g20,temp1),...
else,par1 = polyval(g21,temp1);...  
end;

alfa3=(1)-(7.1/10.1)*(1-par1)

%% MAXIMUM ELEVATION OF THE WAVE PRESSURE

cosb=cos(beta*pi/180);
etas=0.75*(1+cosb)*Hmax;

%% PRESSURE COMPONENT
w0=1.03 %[ton/m3] spesific wight of sea water
p1=(1/2)*(1+cosb)*(alfa1+(alfa2*(cosb)^2))*(w0)*(Hmax);  %[ton/m2]
p2=p1/(1/par1)
p3=alfa3*p1; %[ton/m2]

if etas > hc,...
    p4=p1*(1-(hc/etas)),... %[ton/m2]
else,...
    p4=0  %[ton/m2]
end

% toe pressure
pu=(1/2)*(1+cosb)*alfa1*alfa3*w0*Hmax %[ton/m2]

%% TOTAL PRESSURE AND UPLIFT
hcs=min(etas,hc) %[m] this is hc*
% Total Wave Pressure
P=(0.5*(p1+p3)*(h1))+(0.5*(p1+p4)*hcs) %[ton/m]
% Total Uplift Pressure
U=0.5*(pu)*B %[ton/m3]

%% MOMENT OF WAVE PRESSURE
MP=((1/6)*(2*p1+p3)*(h1^2))+(0.5*(p1+p4)*h1*hcs)+((1/6)*(p1+2*p4)*hcs^2)  %[ton-m/m]

%% MOMENT OF UPLIFT PRESSURE
MU=(2/3)*U*B


Reference:
Goda (1984), Random Seas and Design of Maritime Structures, University of Tokyo Press