function strain1(a,b,c,d)
% function strain1(a,b,c,d)
% Plots an undeformed square of unit half length,
% the deformed parallelogram the square transforms into,
% a unit circle inscribed within the square
% and the strain ellipse for homogeneous deformation described
% by the coefficients a,b,c,d of the deformation gradient matrix F:
%	x'	= ax + by
%	y'	= cx + dy
%
%  [x'] = [a b][x]
%  [y'] = [c d][y]
%
%  [X']  = [F]*[X]
%
% For the figures that are produced, the locations of the points are:
% A: lower left corner
% B: upper left corner
% C: upper right corner
% D: lower right corner
%
% Input parameters (separated by commas in the argument list)
% a = d(x')/dx  (partial derivative of x' with respect to x
% b = d(x')/dy  (partial derivative of x' with respect to y
% c = d(y')/dx  (partial derivative of y' with respect to x
% d = d(y')/dy  (partial derivative of y' with respect to y
% Examples
% strain1(pi/2,pi/2,-pi/2,pi/2)
% strain1(1,1,0,1)
% strain1(1,2,0,1)

% Set the coordinates of the square
A = [-1;-1];
B = [-1;1];
C = [1;1];
D = [1;-1];
X = [A B C D A];	% So x is a 2,5 matrix that represents an initial state

% Plot the square with a dashed black line
clf
plot(X(1,:),X(2,:),'k--')
hold on
axis('equal')
% Label the corners of the square
text(X(1,1),X(2,1),'A')
text(X(1,2),X(2,2),'B')
text(X(1,3),X(2,3),'C')
text(X(1,4),X(2,4),'D')

% Find points around the unit circle
a0 = 1;
thetac = 0:pi/360:2*pi;
xc = a0*cos(thetac);
yc = a0*sin(thetac);
XC = [xc;yc];
% Plot the unit circle with a black dashed line
plot(XC(1,:),XC(2,:),'k--')

% Now perform calculations for the deformation
% Transformation matrix F
F = [a b; c d]
I = [1 0; 0 1];
Ju = F - I
E = 0.5*(F'*F - I)

% Perform the transformation
Xprime = F*X;
XCprime = F*XC;
% Calculate the displacements
U = Ju*X;
UC = Ju*XC;

% Plot the transformed square
plot(Xprime(1,:),Xprime(2,:))
% Label the corners of the transformed square
text(Xprime(1,1),Xprime(2,1),'A*')
text(Xprime(1,2),Xprime(2,2),'B*')
text(Xprime(1,3),Xprime(2,3),'C*')
text(Xprime(1,4),Xprime(2,4),'D*')
% Draw vectors connecting points A,B,C,D to A',B',C',D'
quiveralt(X(1,:),X(2,:),U(1,:),U(2,:),'r');

% Plot the strain ellipse
plot(XCprime(1,:),XCprime(2,:))

% Label the axes and add a title
xlabel('x')
ylabel('y')
title('Initial (dashed) and Final (solid) Objects Under Homogeneous Deformation')