function [ux,uy,uxi,uyi,e1,e2,e1i,e2i] = lycra(x,y,d,xp,yp,a,b,t)
% function lycra(x,y,xp,yp,a,b,t,txy)
% Plots deformation fields for the stretched lycra sheet with a hole.
% Calculates and plots the displacement field, contours the maximum 
% elongation and the minimum elongation, and draws trajectories 
% perpendicular to the direction of the maximum extension.  
% All the input matrices must have the same number of rows and columns.
% Input parameters
% x = matrix of the x-coordinates of the undeformed nodes
% y = matrix of the y-coordinates of the undeformed nodes 
% r = matrix of the radii of the undeformed ellipses 
% xp = matrix of the x-coordinates of the deformed nodes
% yp= matrix of the y-coordinates of the deformed nodes
% a = matrix of the major axes of the deformed ellipses
% b = matrix of the minor axes of the deformed ellipses
% t = matrix of the trends of the minor axes of the def. Ellipses (deg)
% txy= matrix of the angles between the deformed x and y axes (deg)
% txy is measured from the deformed x-axis to the deformed y-axis.
% Output parameters
% u = matrix of displacement in the x-direction
% v = matrix of displacement in the y-direction
% e1 = matrix of the maximum elongation
% e2 = matrix of the minimum elongation
% tau_xy = matrix of the shear strain tau_xy
% Example:
% [x,y] = meshgrid([0:2:8],[0:2:6]]);
% r =  [NaN NaN NaN 1.1 1.1 1.1; NaN NaN 1.1 1.1 1.1;  1.1 1.1 1.1 1.1 1.1; 1.1 1.1 1.1 1.1 1.1];
% xp = [NaN NaN 53.3 72.9; 94.5; NaN NaN 50.9 71.4 94.0; 00.0 24.5 47.7 68.7 92.1; 00.0; 00.0 22.9 46.2 67.1 90.6]
% yp = [0 y.12 y.13 y.14; y.21 y.22 y.23 y.24; y.31 y.32 y.33 y.34; y.41 y.42 y.43 y.44];
% a = [NaN a.12 a.13 a.14; a.21 a.22 a.23 a.24; a.31 a.32 a.33 a.34; a.41 a.42 a.43 a.44];
% b = [NaN b.12 b.13 b.14; b.21 b.22 b.23 b.24; b.31 b.32 b.33 b.34; b.41 b.42 b.43 b.44];
% t = [0 t.12 t.13 t.14; t.21 t.22 t.23 t.24; t.31 t.32 t.33 t.34; t.41 t.42 t.43 t.44]; 
% txy = [90 t.12 t.13 t.14; t.21 t.22 t.23 t.24; t.31 t.32 t.33 t.34; t.41 t.42 t.43 t.44];

% Check to make sure angles are in the range -90 degrees to 90 degrees to
% avoid poor interpolation results (Change made 11/1/06, SJM)
t = mod(t+90,180)-90;

% Calculate the displacement in the x-direction (ux)
ux = xp - x;
% Calculate the displacement in the y-direction (uy)
uy = yp - y;
% Calculate the maximum principal elongations
e1 = (a-d)./d;
% Calculate the minimum principal elongations
e2 = (b-d)./d;
% Calculate the shear strain tau_xy based on the change in right angles
% Make sure any conversions to radians are done correctly.
%tau_xy = tan(deg2rad(90-txy));

% Prepare finer mesh for interpolation
xi = [0:4:80];
yi = [0:4:60];
[xxi,yyi] = meshgrid(xi,yi);
% Interpolate the values for ux, uy, e1, e2, and t
uxi = interp2(x,y,ux,xxi,yyi,'cubic');
uyi = interp2(x,y,uy,xxi,yyi,'cubic');
e1i = interp2(x,y,e1,xxi,yyi,'cubic');
e2i = interp2(x,y,e2,xxi,yyi,'cubic');
ti = interp2(x,y,t,xxi,yyi,'cubic');

% Set values of displacements and strains inside hole to 0
ri = sqrt(xxi.*xxi + yyi.*yyi);
%index = find(ri<32.5);
index = find(ri<36);
uxi(index) = 0*uxi(index);
uyi(index) = 0*uyi(index);
e1i(index) = NaN*e1i(index);
e2i(index) = NaN*e2i(index);
ti(index) = 0*ti(index);

% Prepare plots
% Plot the displacement field
figure(1)
clf     % Clear pre-existing figure
quiver(xxi,yyi,uxi,uyi)
xlabel('x (cm)')
ylabel('y (cm')
axis equal
axis([min(xi) max(xi) min(yi) max(yi)]);
title('Displacement field')

% Contour the field of maximum principal elongation
figure(2)
clf     % Clear pre-existing figure
c = contour(xxi,yyi,e1i);
clabel(c)
xlabel('x (cm)')
ylabel('y (cm')
axis equal
axis([min(xi) max(xi) min(yi) max(yi)]);
title('Maximum principal elongation field')

% Contour the field of maximum principal elongations
figure(3)
clf     % Clear pre-existing figure
d = contour(xxi,yyi,e2i);
clabel(d)
xlabel('x (cm)')
ylabel('y (cm')
axis equal
axis([min(xi) max(xi) min(yi) max(yi)]);
title('Minimum principal elongation field')

% Plot trajectories of the field of minimum principal elongations
figure(4)
clf     % Clear pre-existing figure
traj2(xxi,yyi,sin(ti*pi/180),cos(ti*pi/180))
xlabel('x (cm)')
ylabel('y (cm')
axis equal
axis([min(xi) max(xi) min(yi) max(yi)]);
title('Trajectories of the minimum principal elongations')

% Plot the tau_xy field
%figure(45)
%contour(x,y,tau_xy)
%xlabel('x (cm)')
%ylabel('y (cm')
%axis equal
%title('Shear strain /tau_{xy}')