function geldike(Fz,num)
% Calculates the trajectories perpendicular to the direction
% of the least compressive stress and contours of the
% maximum shear stress for the gelatin dike experiment for
% GG303, assuming no displacements in the direction of the
% short dimension of the tank.
% Fz = the force (weight) of the surface load in Newtons (kg m / sec^2)
% If Fz = 0, one obtains the stresses in a tank with no surface load
% num = number of contours on contour plot

% Set material parameters and grid dimensions
rho = 1.02e03;
g = 9.8;
nu = 0.5;
E = 1332;
G = 444;
%X = -0.245:0.01:0.245;
%X = -0.2475:0.005:0.2475;
X = -0.0975:0.005:0.0975;
%Z = 0:0.01:0.2;
Z = 0:0.005:0.2;
Pair = 1e05;

[x,z] = meshgrid(X,Z);
r = sqrt(x.*x + z.*z);
theta = atan2(x,z);
zip = zeros(size(x));


% Find direction cosines between x,z and r,t axes
axr = x./r;
azr = z./r;
axt = azr;
azt = -axr;

% Find stresses due to gravity
szzg = -rho*g*z - Pair;
sxxg = (nu/(1-nu)) * szzg;
sxzg = zip;
szxg = zip;
syyg = (nu/(1-nu)) * szzg;

% Find stresses due to surface force Fz
%sint = sin(theta);
sint = -z./r;
srrF = (2*Fz/pi)*(sint./r);
srtF = zip;
strF = zip;
sttF = zip;

% Convert stresses due to surface forces to x,y frame
sxxF = axr.*axr.*srrF...
     + axr.*axt.*srtF...
	 + axt.*axr.*strF...
	 + axt.*axt.*sttF;
sxzF = axr.*azr.*srrF...
     + axr.*azt.*srtF...
	 + axt.*azr.*strF...
	 + axt.*azt.*sttF;
szxF = azr.*axr.*srrF...
     + azr.*axt.*srtF...
	 + azt.*axr.*strF...
	 + azt.*axt.*sttF;
szzF = azr.*azr.*srrF...
     + azr.*azt.*srtF...
	 + azt.*azr.*strF...
	 + azt.*azt.*sttF;

% Superpose the gravitational stresses and those due to
% the surface force
sxx = sxxg + sxxF;
sxz = sxzg + sxzF;
szx = szxg + szxF;
szz = szzg + szzF;

% Contour the maximum shear stress
% The formula below can be derived from a Mohr circle
figure(1)
clf
maxshear = (1/2)*sqrt( (sxx-szz).^2 + (sxz-szx).^2 );
colormap(jet)
brighten(jet,0.8)
contourf(x,-z,maxshear,num);
hold on
line([min(X),max(X)],[min(Z),min(Z)]);

xlabel('x')
ylabel('z')
axis('equal');
title('Contours of maximum shear stress')

% Contour the maximum shear stress
% The formula below can be derived from a Mohr circle
figure(2)
clf
maxshear = (1/2)*sqrt( (sxx-szz).^2 + (sxz-szx).^2 );
cc = contour(x,-z,maxshear,num);
clabel(cc);
hold on
line([min(X),max(X)],[min(Z),min(Z)]);

xlabel('x')
ylabel('z')
axis('equal');
title('Contours of maximum shear stress')

% Find the trajectories perpendicular to the least
% compressive stress
figure(3)
clf
ang = 0.5*atan2(sxz,(szz-sxx)/2);
axis('equal');
index = find(x<0);
ang(index) = ang(index) + pi;
traj2(x,-z,cos(ang),sin(ang))
xlabel('x')
ylabel('z')
axis('equal');
title('Direction of most compressive stress')