function [V, D] =stensor_decomp(S)
[V0, D0] = eig(S);
[Y, I] = sort(diag(D0));
D = diag( [D0(I(3),I(3)), D0(I(2),I(2)), D0(I(1),I(1))] );
V = [ V0(:,I(3)), V0(:,I(2)), V0(:,I(1))];
format short
mystring = sprintf('\nThe eigenvalues are:\n');
disp(mystring);
mystring = sprintf('      %6.2f %6.2f %6.2f\n',diag(D));
disp(mystring);
mystring = sprintf('\nThe eigenvectors are:\n');
disp(mystring);
mystring = sprintf('     V1     V2     V3\n');
disp(mystring);
mystring = sprintf('      N: %6.3f %6.3f %6.3f\n',V(1,:));
disp(mystring);
mystring = sprintf('      E: %6.3f %6.3f %6.3f\n',V(2,:));
disp(mystring);
mystring = sprintf('      D: %6.3f %6.3f %6.3f\n',V(3,:));
disp(mystring);