%  Hydrostatic and Deviatoric Stresses

%  Example of a three dimensional state of stress

sigma = [ 2 -5 4
         -5  3  -3
         4  -3  -1]
     
%  Calculate average of the normal stresses

sig_h = mean(diag(sigma))

%  Here is the hydrostatic stress tensor

sigma_hyd = sig_h * eye(3)

%  Here's the deviatoric stress tensor

sigma_dev = sigma - sigma_hyd

%  Now let's get the principal stresses

[V,D] = eigs(sigma);
sigps = sort(diag(D));
s1 = sigps(3);
s2 = sigps(2);
s3 = sigps(1);
taumax = (s1-s3)/2

%  We use this to find the stresses on the octahedral
%  planes.  Coordinate system is referred to principal
%  directions.

noct = 1/sqrt(3) * [ 1 1 1]';
poct = D * noct;
sig_oct = noct' * poct
tau_oct = sqrt( norm(poct)^2 - sig_oct^2)


%  See what happens when we transform the deviatoric stress
% tensor.

sigma_dev_prime = V'*sigma_dev*V;
poct1 = sigma_dev_prime * noct;
sig_oct1 = noct' * poct1
tau_oct1 = sqrt( norm(poct1)^2 - sig_oct1^2)

%  Finally, compare this with the formula on page 34

tau_oct2 = 1/3*sqrt( (s1-s2)^2+(s2-s3)^2+(s3-s1)^2)

