clear

% Kristina Løfman
% Project, due Nov 27 2012
% Calculating factor of safety (FS), shear stress versus shear strength
% Models are based on equations from Sidle&Swanston

% A is image

% constants
Ca = 2.00 % apparent cohesion (kPa), based on Sidle&Swanston
Gammaf = 19.32 % unit weight of soil and field moisture (kN/m3)
z = 0.914 % vertical soil depth. 36 inches
Phi = 43 % internal angle of friction (degrees)
Gammaw = 9.807 % unit weight of water

% Variables
% Alpha = %slope angle (degrees)
h = 0.65 % vertical height of maximum water level in piezometer
% run h for 0.05, 0.1, 0.13, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65,
% 0.7

A = load('clip_output.txt');
A(A == -9999) = NaN; % convert no data values from arc to NaN, to avoid scaling from -9999 to 100

[M,N] = size(A);
FS = zeros(M,N); % matrix

% setting noise elevations (<<<< 1%) to NaN
Aoriginal = A;
A(A < 10) = NaN;

% convert degrees to radians
Ar = A*(pi/180);
Phir = Phi*(pi/180);

% use 2 for-loops, one for rows (M), one for columns (N)
for i = 1:M
    for j = 1:N
        Alpha = Ar(i,j);
        u = Gammaw*h*(cos(Alpha))^2;  %pore water pressure at the failure surface (kPa), depends on precipitation
        S = Ca + (Gammaf*z*(cos(Alpha))^2 - u)*tan(Phir); %shear strength
        T = Gammaf*z*(cos(Alpha))*sin(Alpha);% shear stress
        FS(i,j) = S/T;
    end
end

%plotting
figure(1)
imagesc(Aoriginal)
axis image
title('Slopes (degree)')
colorbar
saveas(gcf,'Slopes','fig');

figure(2)
clims = [ 0 6 ];
imagesc(FS,clims)
axis image
title('FS, h = 0.65')
colormap('hot')
colorbar
saveas(gcf,'h065','fig')

%colorscheme is:
%black = this area was set to NaN (no data) indicating minimal slope
%yellow = safe
%red = highest risk


