%% Raplab - A tool for rip-rap armouring design
clear all
clc
%% Constants 
g=9.81;
%% Log window for input parameters
prompt={'{\bf\fontsize{11}(y)   - Local water depth [m]:}',...
        '{\bf\fontsize{11}(R)   - Centerline radius of curvature of channel bend [m]:}',...
        '{\bf\fontsize{11}(W)   - Width of water surface at upstream end of channel bend [m]:}',...
        '{\bf\fontsize{11}(SF)  - Safety factor (must be >1) [1]:}',...
        '{\bf\fontsize{11}(Cs)  - Stability coefficient [1]} ...                                        (0.30 for angular rock 0.375 for rounded rock):'...
        '{\bf\fontsize{11}(C_T)  - Blanket thickness coefficient given as a function of the uniformity ratio d85/d15 [1]}                                           (Ct=1.0 is recommended):'...
        '{\bf\fontsize{11}(V_[177]) - Channel cross-sectional average velocity [m/s]:}'...
        '{\bf\fontsize{11}(m)   - Bank slope [1:m]:}'...
        '{\bf\fontsize{11}(Sg)  - Relative buoyant density of the protection element [1]:}'...
        '{\bf\fontsize{11}(C_v)  - Velocity distribution coefficient [1]} ...                                  (1) - for straight channels or the inside of bends (1.0); (2) - for the outside of bends (1.283 - 0.2log(Rc/W) and 1.0 for Rc/W > 26); (3) - downstream from concrete channels (1.25); (4) - at the end of dikes  (1.25 )'};
name='US-ARMY rip-rap design parameters';
numlines=1;
defaultanswer={'4','150','15','1.1','0.3','1',' ','1.5','2.65','2'};
options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';
input=inputdlg(prompt,name,numlines,defaultanswer,options);
clear prompt name numlines options defaultanswer
in=str2double(input);
%% Input water depth y
y=in(1,1);
if isnan(y)
    msgbox('Please enter Local water depth "y"', 'Error', 'error')
break
end
%% Input radius and width of the watercourse
rc=in(2,1);
if isnan(rc)
    msgbox('Please enter Centerline radius of curvature of channel bend "rc"', 'Error', 'error')
break
end
W=in(3,1);
if isnan(W)
    msgbox('Please enter Width of water surface at upstream end of channel bend "W"', 'Error', 'error')
break
end
%% Input safety factors, stability and thickness coeff. 
Sf=in(4,1);
if isnan(Sf)
    Sf=1;
end
Cs=in(5,1);
if isnan(Cs)
    Cs=1;
end
Ct=in(6,1);
if isnan(Ct)
    Ct=1;
end
%% Input Average design velocity
Vsr=in(7,1);
if isnan(Vsr)
    msgbox('Please enter Channel cross-sectional average velocity "Vsr"', 'Error', 'error')
break
end
if rc/W<26
    Vproj=Vsr*(1.74-(0.52*log10(rc/W)));
elseif rc/W>26
    Vproj=Vsr;
elseif rc/W==26
    Vproj=Vsr;
end
%% Slope – 1:m shown in degrees
m=1/in(8,1);
if isnan(m)
    msgbox('Please enter Bank slope "m"', 'Error', 'error')
break
end
stupnjeva = atand(m) %tan-1(m)
%% Spec. weight of rip-rap
Sg=in(9,1);
if isnan(Sg)
    Sg=2.65
end
%% Velocity distribution coeficient
krb=in(10,1);
Cv2=1.283-(0.2*log10(rc/W)); %Koeficijent raspodjele brzine [1]
if krb==1
    Cv=1;
elseif krb==2
    Cv=Cv2;
elseif krb==3
    Cv=1.25;
elseif krb==4
    cv=1.25;
elseif isnan(krb)
    msgbox('Please enter Velocity distribution coefficient', 'Error', 'error')
break
end
%% Parameters for Pilarczyka
prompt2={'{\bf\fontsize{11}(\phi_{SC})  - Stability correction factor [1]:}                                   {\it(a) - exposed edges of gabions/stone matressess (1.0); (b) - exposed edges of rip-rap and armourstone (1.5), (c) - continuous rock protection (0.75), (d) - interlocked blocks and cabled blockmats (0.5)}',...
        '{\bf\fontsize{11}(\psi_C_R) - Mobility parameter [1]:}                                                {\it(a) - rip-rap and armourstone (0.035); (b) - box gabions and gabion matresses (0.070); (c) - rock fill in gabions (<0.100)}',...
        '{\bf\fontsize{11}(k_t^2)   - Turbulence factor [1]:}                                                 {\it(a) - normal turbulence level (kt^2 = 1.0); (b) - non-uniform flow, increased turbulence in outer bends (kt^2 = 1.5); (c) - non-uniform flow, sharp outer bends (kt^2 = 2); (d) - non-uniform flow, special cases (kt^2 > 2)}',...
        '{\bf\fontsize{11}(D_n) - Assumed diameter of rip-rap [m]}',...
        '{\bf\fontsize{11}(x) - Roughness factor}                                                            {\it(k_s from 1 to 3)}',...
        '{\bf\fontsize{11}(\beta) - Longitudinal slope of the watercourse [^o]:}',...
        '{\bf\fontsize{11}(\phi) - Rip-rap angle of response [^o]}'};
name='Pilarczyk rip-rap design parameters';
numlines=1;
defaultanswer={'0.75','0.035','1','0.38','1','2','40'};
options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';
input_Pilarczyk=inputdlg(prompt2,name,numlines,defaultanswer,options);
clear prompt name numlines options defaultanswer
in_Pilarczyk=str2double(input_Pilarczyk);
%% Pilarczyk parameter processing
phiSC=in_Pilarczyk(1,1);
if isnan(phiSC)
    msgbox('Please enter Stability correction factor (\phi_{SC})', 'Error', 'error')
break
end
psi=in_Pilarczyk(2,1);
if isnan(psi)
    msgbox('Please enter critical Mobility parameter of the protection element (\psi)', 'Error', 'error')
break
end
kt2=in_Pilarczyk(3,1);
if isnan(kt2)
    msgbox('Please enter Turbulence factor (k_t^2)', 'Error', 'error')
break
end
Dn=in_Pilarczyk(4,1);
if isnan(Dn)
    msgbox('Please assume rip-rap diameter (D_n)', 'Error', 'error')
break
end
x=in_Pilarczyk(5,1);
if isnan(x)
    msgbox('Please enter roughness factor (x)', 'Error', 'error')
break
end
ks=x*Dn;
beta=in_Pilarczyk(6,1);
if isnan(beta)
    msgbox('Please enter Longitudinal slope of the watercourse (\beta)', 'Error', 'error')
break
end
phi=in_Pilarczyk(7,1);
if isnan(phi)
    msgbox('Please enter rip-rap angle of response ', 'Error', 'error')
break
end
%% Pilarczyk velocity profile
% Profil brzine Pitanje-odgovor
choice = questdlg('Define velocity profile factor', ...
	'Velocity Menu', ...
	'Fully developed logaritmic velocity profile','Not-fully developed logaritmic velocity profile','Fully developed logaritmic velocity profile');
% Handle response
switch choice
    case 'Fully developed logaritmic velocity profile'
        disp([choice ' selected.'])
        velprofile = 1;
    case 'Not-fully developed logaritmic velocity profile'
        disp([choice ' selected.'])
        velprofile = 2;
end
if velprofile==1
    kh=(2/(log10((1+(12*y)/ks))^2));
elseif velprofile==2
    kh=(((1+y)/Dn)^(-0.2));
end
delta=Sg-1;
kd=sqrt(1-(sind(stupnjeva)^2/sind(phi)^2)); %side slope factor
kl=sind(phi-beta)/sind(phi);
ksl=kd*kl;

%% Escarameia & May
prompt3={'{\bf\fontsize{11}(x_b) - Velocity factor [1]:}                                               {\it recommended x_b = 0.74 to 0.90 U}'...
        '{\bf\fontsize{11}(r) - Turbulence level [1]:}                                                {\it(a) - straight river or channel reaches, normal (low) (r=0.12); (b) - Edges of revetments in straight reaches, normal (high) (r=0.20); (c) - Bridge piers, caissons and spur dikes and transitions, medium to high (r=0.35 - 0.50); (d) - Downstream of hydraulic structures, very high (r=0.60)}'...
        '{\bf\fontsize{11}(Sf_{HEC11}) - HEC-11 Stability factor level [1]:}                                                {\it(a) - uniform flow, Rc/W>30 (Sf_{HEC11} = 1.0 to 1.2; (b) - gradualy varying flow, 10<Rc/W<3 (Sf_{HEC11} = 1.3 to 1.6; (c) - rapidly varying flow, Rc/W<10 (Sf_{HEC11} = 1.6 to 2.0)}'};
name='HEC-11, Escarameia and May rip-rap design parameters';
numlines=1;
defaultanswer={'0.74','0.15','1','2'};
options.Resize='on';
options.WindowStyle='normal';
options.Interpreter='tex';
input_EM=inputdlg(prompt3,name,numlines,defaultanswer,options);
clear prompt name numlines options defaultanswer
in_EM=str2double(input_EM);

%% Escarameia & May parameters processing
xb=in_EM(1,1);
if isnan(xb)
    msgbox('Please enter Velocity factor (x_{b})', 'Error', 'error')
break
end
r=in_EM(2,1);
if isnan(r)
    msgbox('Please enter Turbulence level (r)', 'Error', 'error')
break
end
Sf_hec11=in_EM(3,1);
if isnan(Sf_hec11)
    msgbox('Please enter HEC-11 Stability factor level Sf_{HEC11})', 'Error', 'error')
break
end
ub=xb*Vsr;
ct_EM=12.3*r-0.2;

%% Calculation for US_ARMY_EM_1601
tic
K1=sqrt(1-((sind(stupnjeva-14)/sind(32))^1.6));
d30=y*Sf*Cs*Ct*Cv*((Vproj/sqrt(K1*(Sg-1)*g*y))^2.5);
US_ARMY_EM_1601=1.2*d30*1000;
%% Pilarczyk calculation
d_Pilarczyk=phiSC/delta*0.035/psi*kh/ksl*kt2*((Vsr)^2)/2/g;
d50_Pilarczyk=d_Pilarczyk/0.84*1000;
%% Calculation for EM
d50_EM=ct_EM*(ub^2)/2/g/delta*1000;
%% Calculation for HEC11
K1_hec11=kd;
Csg_hec11=2.12/(delta^1.5);
Csf_hec11=(Sf_hec11/1.2)^1.5;
d50_HEC11=0.00594*Csg_hec11*Csf_hec11*(Vsr^3)/sqrt(y)/(K1_hec11^1.5)*1000;
%% Record input data
myfolder = uigetdir;
f1 = fullfile(myfolder,'inputparameters.xls');

inpa={'(y)   - Local water depth [m]:';
        '(R)   - Centerline radius of curvature of channel bend [m]:';
        '(W)   - Width of water surface at upstream end of channel bend [m]:';
        '(SF)  - Safety factor (must be >1) [1]:';
        '(Cs)  - Stability coefficient [1]:';
        '(C_T)  - Blanket thickness coefficient given as a function of the uniformity ratio d85/d15 [1]:';
        '(Vsr) - Channel cross-sectional average velocity [m/s]:';
        '(m)   - Bank slope [1:m]:';
        '(Sg)  - Relative buoyant density of the protection element [1]:';
        '(C_v)  - Velocity distribution coefficient [1]:';
        '(phiSC)  - Stability correction factor [1]:';
        '(psiCR) - Mobility parameter [1]:';
        '(kt^2)   - Turbulence factor [1]:';
        '(Dn) - Assumed diameter of rip-rap [m]:';
        '(x) - Roughness factor [1]:';
        '(beta) - Longitudinal slope of the watercourse [o]:';
        '(phi) - Rip-rap angle of response [o]:';
        '(x_b) - Velocity factor [1]:';
        '(r) - Turbulence level [1]:'
        '(Sf_HEC11) - HEC-11 Stability factor level [1]:'};

inpall=[input(1,1);input(2,1);input(3,1);input(4,1);input(5,1);input(6,1);input(7,1);input(8,1);input(9,1);input(10,1); 
        input_Pilarczyk(1,1);input_Pilarczyk(2,1);input_Pilarczyk(3,1);input_Pilarczyk(4,1);input_Pilarczyk(5,1);input_Pilarczyk(6,1);input_Pilarczyk(7,1);
        input_EM(1,1);input_EM(2,1);input_EM(3,1)];

    inpc=[inpa, inpall];
    
xlswrite(f1,inpc);

%% Graf
% Graphics plot Question-answer
choice = questdlg('Chose graph', ...
    'Graph Menu', ...
	'Normal scale on "x" and "y" axis','Logaritmic scale on "x" and "y" axes','Logaritmic scale on "y" axis','Logaritmic scale on "x" axis');
% Handle response
switch choice
    case 'Normal scale on "x" and "y" axis'
        disp([choice ' coming right up.'])
        mjerilo = 1;
    case 'Logaritmic scale on "x" and "y" axes'
        disp([choice ' coming right up.'])
        mjerilo = 2;
    case 'Logaritmic scale on "y" axis'
        disp([choice ' coming right up.'])
        mjerilo = 3;
    case 'Logaritmic scale on "x" axis'
        disp([choice ' coming right up.'])
        mjerilo = 4;
end
xpl = 0.25:0.25:7; %korak
%US ARMY
if rc/W<26
    ypl_USARMY = y*Sf*Cs*Ct*Cv*(((xpl*(1.74-(0.52*log10(rc/W))))/sqrt(K1*(Sg-1)*g*y)).^2.5)*1.2*1000;
elseif rc/W>26
    ypl_USARMY = y*Sf*Cs*Ct*Cv*((xpl/sqrt(K1*(Sg-1)*g*y)).^2.5)*1.2*1000;
elseif rc/W==26
    ypl_USARMY = y*Sf*Cs*Ct*Cv*((xpl/sqrt(K1*(Sg-1)*g*y)).^2.5)*1.2*1000;
end
%PILARCZYK
ypl_Pilarczyk = (phiSC/delta*0.035/psi*kh/ksl*kt2*((xpl).^2)/2/g)/0.84*1000;
%EM
ypl_EM = ct_EM*((xb*xpl).^2)/2/g/delta*1000;
%HEC11
ypl_hec11=0.00594*Csg_hec11*Csf_hec11*(xpl.^3)/sqrt(y)/(K1_hec11^1.5)*1000;
%uvjeti max(ypl_USARMY,ypl_Pilarczyk,ypl_EM,ypl_hec11)
max1=max(ypl_USARMY);
max2=max(ypl_Pilarczyk);
max3=max(ypl_EM);
max4=max(ypl_hec11);

maxy=[max1 max2 max3 max4];

maxyy=max(maxy)+100;

if mjerilo==1
    plot(xpl,ypl_USARMY,':k^', xpl,ypl_Pilarczyk, '--ko',xpl,ypl_EM,'-.kd',xpl,ypl_hec11,'-ks')
    set(gca,'XTick',0:0.5:7,'YTick',0:100:max(maxyy)) %x os od nula : s korakom : do 7
elseif mjerilo==2
    loglog(xpl,ypl_USARMY,':k^', xpl,ypl_Pilarczyk, '--ko',xpl,ypl_EM,'-.kd',xpl,ypl_hec11,'-ks')
elseif mjerilo==3
    semilogy(xpl,ypl_USARMY,':k^', xpl,ypl_Pilarczyk, '--ko',xpl,ypl_EM,'-.kd',xpl,ypl_hec11,'-ks')
elseif mjerilo==4
    semilogx(xpl,ypl_USARMY ,':k^', xpl,ypl_Pilarczyk, '--ko',xpl,ypl_EM,'-.kd',xpl,ypl_hec11,'-ks')
end
hlegend=legend('US-ARMY-EM-1601','Pilarczyk','Escarameia and May','HEC-11');
title('Rip-rap median diameter relative to flow velocity');
xlabel('Flow velocity [m/s]');
ylabel('Median rip-rap diameter d_{50} [mm]');
set(hlegend,'FontAngle','italic','Location','Best')
%% Plot results
disp('rezultati:')
str1=['US_ARMY_EM_1601: ', 'd50= ', num2str(US_ARMY_EM_1601), '[mm]'];
disp(str1)
str2=['Pilarczyk: ', 'd50= ', num2str(d50_Pilarczyk), '[mm]'];
disp(str2)
str3=['Escarameia and May: ', 'd50= ', num2str(d50_EM), '[mm]'];
disp(str3)
str4=['HEC-11: ', 'd50= ', num2str(d50_HEC11), '[mm]'];
disp(str4)
msgbox({num2str(str1);num2str(str2);num2str(str3);num2str(str4)}, 'Rezultati:')

%% Record output results
f2 = fullfile(myfolder,'DesignDiameters.xls');
inpa2={'US ARMY - EM 1601 d50 =';
        'Pilarczyk d50 =';
        'Escarameia and May d50 =';
        'HEC-11 d50 ='};
inpall2={US_ARMY_EM_1601; d50_Pilarczyk; d50_EM; d50_HEC11};
unit = {'mm';'mm';'mm';'mm'};
inpc2=[inpa2, inpall2, unit];
xlswrite(f2,inpc2)
%% End
toc
disp('***Normal run complete***')
clock