function [V,R,eigvecV,eigvalV,eigvecF,eigvalF,eigvecV_] = vr_decomp(F)
%VR decompostion of deformation gradient matrix F
% function [V,R,eigval,eigvec] = vr_decomp(F)  
% This function provides the VR decomposition of a deformation gradient
% matrix F,provides the magnitudes of the principal strain axes, 
% provides the direction cosines for the principal strain axes, plots the
% strain ellipse, and plots the reciprocal strain ellipse
% Input parameters
% F = 2D deformation gradient matrix = [a b;c d]
% Output parameters
% V = symmetric 2D deformation gradient matrix = [A B;C D]
% R = rotation matrix = [A B;C D]
% eigvalV = eigenvalues of V
% eigvecV = eigenvectors of V
% Examples
% [V,R,eigvecV,eigvalV,eigvecF,eigvalF,eigvecV_] = vr_decomp([2,1;0,1]);
% [V,R,eigvecV,eigvalV,eigvecF,eigvalF,eigvecV_] = vr_decomp([2,1;1,1]);

disp('Deformation gradient matrix [F] = [V][R]')
disp(F)

% Perform other calculations
V=sqrtm(F*F');
V_=sqrtm(F'*F);
R =V\F;
[eigvecV,eigvalV] = eig(V);
[eigvecF,eigvalF] = eig(F);
[eigvecV_,eigvalV_] = eig(V_);
[eigvecF_,eigvalF_] = eig(inv(F));

disp('Symmetric component [V]')
disp(V)

disp('Rotation matrix [R]')
disp(R)

disp('Eigenvalues of [V]')
disp(eigvalV)

disp('Eigenvectors of [V] (in columns)')
disp(eigvecV)

disp('Eigenvalues of inverse of [V]')
disp(eigvalV)

disp('Eigenvectors of inverse of [V]')
disp(eigvecV)

disp('Eigenvalues of [F]')
disp(eigvalF)

disp('Eigenvectors of [F] (in columns)')
disp(eigvecF)

% Calculate positions on unit circle
t = 0:pi/180:2*pi;
x = cos(t);
y = sin(t);
X = [x;y];
% Calculate positions on strain ellipse
Xp = F*X;
xp = Xp(1,:);
yp = Xp(2,:);
% Calculate displacement of points on strain ellipse
ux = xp - x;
uy = yp - y;
% Calculate positions on reciprocal strain ellipse
Xpp = inv(F)*X;
xpp = Xpp(1,:);
ypp = Xpp(2,:);
% Calculate displacement of points on reciproal strain ellipse
ux_ = xpp - x;
uy_ = ypp - y;

nt = 6;
dt = (1:360/nt:361-360/nt);

% Plot unit circle, strain ellipse, and selected displacements
subplot(2,2,1)
plot(x,y,xp,yp)
hold on
quiver(x(dt),y(dt),ux(dt),uy(dt),0)
hold off
axis equal
title('Unit circle, strain ellipse, & selected displacements')

% Plot unit circle, strain ellipse, principal axes, and eigenvectors
subplot(2,2,2)
plot(x,y,xp,yp)
hold on
quiver([0,0],[0,0],eigvecV(1,:),eigvecV(2,:),0)
quiver([0,0],[0,0],eigvecF(1,:),eigvecF(2,:),0)
hold off
axis equal
title('Strain ellipse, eigenvectors (blue), principal axes (red)')

% Plot unit circle, reciprocal strain ellipse, and selected displacements
subplot(2,2,3)
plot(x,y,xpp,ypp,'k')
hold on
quiver(x(dt),y(dt),ux_(dt),uy_(dt),0)
hold off
axis equal
title('Unit circle, reciprocal strain ellipse, & selected displacements')

% Plot unit circle, reciprocal strain ellipse, reciprocal principal axes, 
% and reciprocal eigenvectors
subplot(2,2,4)
plot(x,y,xpp,ypp,'k')
hold on
quiver([0,0],[0,0],eigvecF_(1,:),eigvecF_(2,:),0)
quiver([0,0],[0,0],eigvecV_(1,:),eigvecV_(2,:),0)

hold off
axis equal
title('Recip. strain ellipse, eigenvectors (green), principal axes (red)')


end

