% displacement factor is alpha * t where alpha is a factor and t is time
alphat = 5e-1;
F = [1 alphat 0; 0 1 0; 0 0 1];
H = F - I
I = eye(3)
C = transpose(F) * F
B = F * transpose(F)
fprintf(1, 'polar decompositions\n')
U = sqrtm(C)
R = F / U
V = F * R'

% various strains
fprintf(1, 'strain tensors\n')
fprintf(1, 'extact strain Um1\n')
Um1 = U - I
fprintf(1, 'Green St Venant strain G\n')
G = (C - I) / 2
fprintf(1, 'Infinitesimal strain E\n')
E = 0.5 * (H + H')

% rotations
fprintf(1, 'rotation exact R\n')
R
fprintf(1, 'infinitesimal rotaiton W, approximating R - I:\n')
W = 0.5 * (H - H')

%%% volumetric strain
fprintf(1, 'J and volumetric strains\n')
J = det(F)
fprintf(1, 'exact vol strain J - 1\t')
epsVexact = J - 1
fprintf(1, 'infinitesimal vol strain J - 1\t')
epsVinfinitesimal = trace(E)

%%%% surface transfer matrix
sx2sy = J * transpose(inv(F))

% sample vectors of directions used to calculate change of segments, areas,
% etc
cntr = 1;
vecs{cntr} = [1 0 0];
cntr = cntr + 1;
vecs{cntr} = [0 1 0];
cntr = cntr + 1;
vecs{cntr} = [0 0 1];
cntr = cntr + 1;
vecs{cntr} = [1 1 1];
%
sz = length(vecs);

% make vectors of size 1 so it's easier to interpret results
for i = 1:sz
    vecs{cntr} = vecs{cntr} / norm(vecs{cntr});
end


for i = 1:sz
    vec = vecs{i}';
    fprintf(1, 'outputs for vector # %d: [%f, %f, %f]\n', i, vec(1), vec(2), vec(3));
    fprintf(1, 'dx (= vec) to dy, dy = F dx:\n')
    dy = F * vec
    fprintf(1, 'dSx (= vec) to dSy, dSy = det F F^-T dSx:\n')
    dSy = F * vec
    
    fprintf(1, 'normal strains along n = vec\n')
    fprintf(1, 'Exact strain Un\n')
    strainExact = Um1 * vec
    fprintf(1, 'Green St-Venant strain Gn\n')
    strainG = G * vec
    fprintf(1, 'Infinitesimal strain Gn\n')
    strainE = E * vec
end