% The third step is to use the function fmincon embedded in MATLAB constraint
%to search for the minimum factor of safety. 

%The function will take the given information of soil and slope as 
%input data and return the center of the critical circle, the 
%corresponding minimum factor of safety.

%The function considers two kinds of different critical mechanism for flat
%slopes and steep slopes.

function [x,FS]=slope_u(D,H,beta,c,gamma)

D=1.5;
H=10.0;
beta=52.1;
c=50.0;
gamma=20.0;

B = beta/180*pi;

if beta==90
    xt = 0;
    yt = -H;
else    
    xt = H/tan(B);
    yt = -H;
end
% seperate three parts to set up an initial start-searching point, get
% three optimization values respectively, and then compare the reslus,
% take the minor.

if beta>=45
    for i=1:3 
%initial start-searching point on the left of toe
        if i==1
            x0 = [0,H];  
%initial start-searching point on the line above toe
        elseif i==2
%if beta is less than 9 degree, need to increase y0
            y0=H+30*H*abs(tan(B-pi/3))
            x0 = [xt,y0]; 
%initial start-searching point on the right of toe
        else
            x0 = [xt+H,y0]; 
        end

%call the function of constraint and FoS as constraint and objective 
%function to use in fmincon. @(x) is function handle in MATLAB
        nonlcons = @(x)constraint_u(x,D,H,beta);
        obj = @(x)FoS_u(x,D,H,beta,c,gamma);
        options = optimoptions('fmincon','MaxFunEvals',10000,'TolCon',1e-16,...
           'TolX',1e-32,'MaxIter',2000);
        [x,FS(i)]=fmincon(obj,x0,[],[],[],[],[],[],nonlcons,options)
        y(i,:)=x
    end
    if FS(1)<FS(2) && FS(1)<FS(3)
        x=y(1,:)
    elseif FS(2)<FS(1) && FS(2)<FS(3)
        x=y(2,:)
    else
        x=y(3,:)
    end   
    FS=min(FS)
else
    x0 = [0,H]; 
 %call the function of constraint and FoS as constraint and objective 
%function to use in fmincon. @(x) is function handle in MATLAB
    nonlcons = @(x)constraint_u(x,D,H,beta);
    obj = @(x)FoS_u(x,D,H,beta,c,gamma);
    options = optimoptions('fmincon','MaxFunEvals',10000,'TolCon',1e-16,...
       'TolX',1e-32,'MaxIter',2000);
    [x,FS]=fmincon(obj,x0,[],[],[],[],[],[],nonlcons,options)
end
    
%get (x1,y1) and (x2,y2) for drawing the geometric properties of slope and 
%the critical circle
if x(1)<xt
    R = x(2)+D*H;
       
    x1 = x(1)-sqrt(R^2-x(2)^2);

    if (R^2-(yt-x(2))^2) > 0
     xh = x(1)+sqrt(R^2-(yt-x(2))^2);
    else
     xh = 0;
    end

    if xh > xt
        x2 = xh;
        y2 = -H;
    else if ((x(1)-x(2)*tan(B))^2-(1+(tan(B))^2)*(x(1)^2+x(2)^2-R^2))>0
        x2 = (x(1)-x(2)*tan(B)+sqrt((x(1)-x(2)*tan(B))^2 ...     
        -(1+(tan(B))^2)*(x(1)^2+x(2)^2-R^2)))/(1+(tan(B))^2);
        y2 = -x2*tan(B);
        else
            x2=xt;
            y2=yt;
        end
    end
else
    R = sqrt((x(2)-yt)^2+(xt-x(1))^2);
    
    x1 = x(1)-sqrt(R^2-x(2)^2);
    y1 = 0;

    x2=xt;
    y2=yt;
end

%drawing geometric properties of slope 
set(gcf,'color','w');

if D==1.0
    dl=-2.0; dr=0.0;
elseif D==2.0
    dl=-3.0; dr=2.0;
elseif D==4.0
    dl=-5.0; dr=4.0;
end


dl=-20; dr=20;
xp1=[dl 0 H/tan(B) H/tan(B)+dr];
axis off;
yp1=[   0 0      -H    -H];
line(xp1,yp1,'color','black');
line([dl H/tan(B)+dr],[-D*H -D*H],'color','black');
xp2=linspace(x1,x2);
yp2=x(2)-sqrt(R^2-(xp2-x(1)).^2);
line(xp2,yp2,'color','red','LineWidth',2);
%xlim([dl,H/tan(B)+dr]);
axis equal
hold on
end

