% lab_09_ex2c_2008
% lab_09_ex2c_2008
% Prepares plots for Lab 9, exercise 2c, 2008
% Reference frame is x= east, y = north, z = up.
% For this problem a shear zone develops above a fault.
% The fault strikes 90 degrees, is vertical, is very long, and
% is very deep.  The distance from the horizontal top edge of the
% fault to the surface is d.  The slip across the fault is considered
% to be uniform and or magnitude b.  The displacement field is for a screw
% dislocation.
%  
% Find initial positions of points on a grid
x = 0:1:1;
y = 0:1:3;
[Xmesh,Ymesh] = meshgrid(x,y);
% Rearrange x- and y- coordinates into row vectors
xvec = Xmesh(:)';
yvec = Ymesh(:)';

% Define F
% xp = 1.5*x +  1*y
% yp = 0.5*x +  1*y
F = [1.5,1;0.5,1];


% Find positions on deformed grid
Xp_Yp = F*([xvec;yvec]);
% Separate out the transformed x,y coordinates from Xp
% as row vectors
xpvec = Xp_Yp(1,:);
ypvec = Xp_Yp(2,:);
% Rearrange the points one more time to put them on a mesh
Xpmesh = reshape(xpvec,size(Xmesh));
Ypmesh = reshape(ypvec,size(Ymesh));


% Prepare a reference circle of radius 0.5
theta = deg2rad(0:1:360);
xc = 0.5*cos(theta);
yc = 0.5*sin(theta);

% Plot the undeformed grid
figure(1)
plot(Xmesh,Ymesh,'--')
hold on
plot(Xmesh',Ymesh','--')

% Plot deformed grid
plot(Xpmesh,Ypmesh);
plot(Xpmesh',Ypmesh');

% Plot displacements
uxvec = xpvec - xvec;
uyvec = ypvec - yvec;
quiver(xvec,yvec,uxvec,uyvec,0)

% Plot unit circles centered at various point undeformed grid
% At coordinates x,y = 0.5,0.5 
plot(xc+0.5,yc+0.5,'--')
hold on
% At coordinates x,y = 0.5,1.5 
plot(xc+0.5,yc+1.5,'--')
% At coordinates x,y = 0.5,2.5 
plot(xc+0.5,yc+2.5,'--')

% Plot strain ellipses
% For original coordinates x,y = 0.5,0.5
Xp_Yp_a = F*([xc+0.5;yc+0.5]);
plot(Xp_Yp_a(1,:),Xp_Yp_a(2,:))

% For original coordinates x,y = 0.5,1.5
Xp_Yp_a = F*([xc+0.5;yc+1.5]);
plot(Xp_Yp_a(1,:),Xp_Yp_a(2,:))
% For original coordinates x,y = 0.5,2.5
Xp_Yp_a = F*([xc+0.5;yc+2.5]);
plot(Xp_Yp_a(1,:),Xp_Yp_a(2,:))

hold off
axis equal

% Label the plot
title ('Lab 9, Exercise 2c, 2008')
xlabel('x')
ylabel('y')