function lab14(b,a)
% function lab14(b,a).  Draws a profile of predicted displacement at
% the ground surface as a function of distance from a long vertical
% strike-slip fault with constant slip using a screw dislocation model.
% Parameter "b" is the slip across the fault (in meters).
% The slip is TWICE the fault-parallel component displacement 
% along one of the fault walls!
% Parameter "a" is the depth of the lower edge of the dislocation (in km). 
% Both parameters "b" and "a" must be placed between parentheses.
% For example, to start and just see the data type
% lab12(0,0)
% To get model curves you need to provide non-zero values for "b" and "a".
% If your curve is below the data, the slip and/or fault depth is too low.
% If your curve is above the data, the slip and/or fault depth is too high.
% Plots will be superposed.  To clear the screen to start over type
% clf
% The surface displacements are elastic displacements calculated
% using a screw dislocation solution (see lecture 23).  
% The displacements are calculated along a horizontal plane
% that bisects a vertical screw dislocation in an infinite body.
% This dislocation extends from a depth of "a" km below the surface
% to "a" km above the surface.
% The horizontal plane represents the surface of a half-space,
% and here that is the ground surface.
% Slip across the dislocation results in no tractions on this
% plane (i.e., no normal and shear stresses act ON this plane),
% so the displacements on or below this plane are appropriate
% for those in the Earth around the central portion of a long vertical 
% strike slip fault with a constant slip.
% Data for fault-parallel displacements (with error bars) are from the
% 1906 San Francisco earthquake as reported by Pollard and Segall (1987).
% The reference frame has the x-axis vertical and in the plane of the fault.
% The y-axis is normal to the fault and at the ground surface.
% The z-axis is horizontal and parallels fault strike.
% Estimate the slip to +/- 1 meter and the depth of faulting to +/- 5 km.

% Set the grid to calculate displacements on
y = 0:0.1:14;
x = zeros(size(y));


% Calculate displacement w parallel to the fault
w = (b/(2*pi)) * ( atan2(y,(x-a)) - atan2(y,(x+a)) );

% 1906 Displacement data
y6 = [0.18, 0.18, 0.18]; 	w6 = [2.05, 2.45, 2.87];
y5 = [0.50, 0.50, 0.50];	w5 = [2.11, 2.50, 2.91];
y7 = [1.48, 1.48, 1.48];	w7 = [1.69, 2.09, 2.50];
y4 = [3.65, 3.65, 3.65];	w4 = [1.43, 1.83, 2.23];
y3 = [3.92, 3.92, 3.92];	w3 = [1.38, 1.79, 2.19];
y8 = [5.72, 5.72, 5.72];	w8 = [1.15, 1.55, 1.95];
y9 = [6.40, 6.40, 6.40];	w9 = [0.97, 1.36, 1.79];
y10= [6.71, 6.71, 6.71];	w10 = [1.08, 1.48, 1.89];
y11= [6.82, 6.82, 6.82];	w11 = [1.28, 1.70, 2.10];
y12= [7.66, 7.66, 7.66];	w12 = [1.05, 1.45, 1.85];
y2= [11.26, 11.26, 11.26];	w2 = [0.60, 1.00, 1.41];
y1= [13.56, 13.56, 13.56];	w1 = [0.60, 1.00, 1.41];

% Plot 1906 data
figure(1)
plot  ( y6,w6,'-',y5,w5,'-',y7,w7,'-',y4,w4,'-',y3,w3,'-',y8,w8,'-',...
		y9,w9,'-',y10,w10,'-',y11,w11,'-',y12,w12,'-',y2,w2,'-',y1,w1,'-')
hold on
plot  ( y6(2),w6(2),'o',y5(2),w5(2),'o',y7(2),w7(2),'o',y4(2),w4(2),'o',...
		y3(2),w3(2),'o',y8(2),w8(2),'o',y9(2),w9(2),'o',y10(2),w10(2),'o',...
		y11(2),w11(2),'o',y12(2),w12(2),'o',y2(2),w2(2),'o',y1(2),w1(2),'o')

if b~=0	
	% Plot model curve
	plot (y,w)
	aa = num2str(a);
	bb = num2str(b);
	text(y(100),w(100)+0.05,['a=',aa,' km, b=',bb,' m'])
end

xlabel('Distance from fault (km)')
ylabel('Displacement parallel to fault (m)')
title('1906 Displacements - Point Arena')