function FS = factor_of_Safety(c_prime, h, theta, v_r, rho_r, m, phi, v)
%FACTOR_OF_SAFETY
%
% Usage:
% FS = factor_of_safety(c_prime, h, theta, vr, rho_r, m, phi, v)
%
% This function calculates the factor of safety for an infinite slope
% stability problem.
%
% See Arrowsmith notes: http://arrowsmith362.asu.edu/Lectures/Lecture15/
%
% Author:
% JRA, 11/8/2006
%
% Variables:
%        External
% c_prime is the (effective) cohesion (as reduced by loss of surface tension)
% h is the thickness of the potential slide
% theta is the dip angle of the potential failure plane
% rho_b = v_r rho_r + m(1-v_r)rho_w is the wet bulk density
% v_r is the volume fraction of solid material; (1-v_r) is the porosity
% or about 0.3 to 0.4, so v_r ~ 0.6-0.7
% rho_r is rock density (2670 kg/m^3)
% m is the portion of saturated thickness of the slide; m = 1 for a
% fully saturated slide and m = 0 for a completely dry slide
% phi is the angle of internal friction
% v stands for verbose. If v ==0, then supress output. 
%
%       Internal
% g is acceleration of gravity
% rho_w is water density (1000 kg/m^3)
% theta_rads is theta in radians
% phi_rads is phi in radians
% denominator = driving stresses
% numerator = resisting stresses

%write the constants back to make sure that it looks right
if v ~=0
    fprintf(1,'Constants\n')
    fprintf(1,'c = %.3f, h = %.3f, theta = %.3f, rho_r = %.3f, m = %.3f, phi = %.3f\n', c_prime,h,theta,rho_r,m,phi);
end

%Specify constants
g = 9.81;
rho_w = 1000;
theta_rads = theta.*2.*pi./360;
phi_rads = phi.*2.*pi./360;

%Compute wet bulk desnity of debris load
rho_b = v_r.*rho_r + m.*(1-v_r).*rho_w;

%Resisting forces
numerator = c_prime + h.*g.*(cos(theta_rads).^2).*(rho_b-rho_w.*m).*tan(phi_rads);

%Driving forces
denominator = rho_b.*h.*g.*sin(theta_rads).*cos(theta_rads);

%Factor of safety
FS = numerator./denominator;

if v ~=0
    fprintf(1,'Results\n')
    fprintf(1,'rho_b = %.3f, FS = %.3f, theta = %.3f, rho_r = %.3f, m = %.3f, phi = %.3f\n', rho_b,FS);
    if FS>1
            fprintf(1,'SAFE!\n');
        elseif FS<1
            fprintf(1,'UNSAFE!\n');
        elseif FS==1
            fprintf(1,'Critical!\n');
        else
            fprintf(1,'Something is wrong\n');
    end
end


