Added missing file

compute_eigenvalues_of_tensor3d.m

@@ -0,0 +1,56 @@
+function [lambda1 lambda2 lambda3] = compute_eigenvalues_of_tensor3d(t11, t12, t13, t22, t23, t33)
+% COMPUTE_EIGENVALUES_OF_TENSOR3D Estimate the eigenvalues of the real symmetric 3x3 matrix T
+%
+% [lambda1 lambda2 lambda3] = compute_eigenvalues_of_tensor3d(t11, t12, t13, t22, t23, t33)
+%
+
+% Copyright (c) 2012 Daniel Forsberg
+% [email protected]
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+% Taken from http://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices
+% Matrix A is replaced with the elements t11, t12, t13, t22, t23 and t33.
+% Also added part to always make sure that eigenvalues a positive
+
+q = (t11 + t22 + t33)/3;
+p = (t11 - q).^2 + (t22 - q).^2 + (t33 - q).^2 + 2 * (t12.^2 + t13.^2 + t23.^2);
+p = sqrt(p / 6) + eps;
+B11 = (1 ./ p) .* (t11 - q);
+B12 = (1 ./ p) .* (t12);
+B13 = (1 ./ p) .* (t13);
+B22 = (1 ./ p) .* (t22 - q);
+B23 = (1 ./ p) .* (t23);
+B33 = (1 ./ p) .* (t33 - q);
+
+r = (- B33.*B12.^2 + 2*B12.*B13.*B23 - B22.*B13.^2 - B11.*B23.^2 + B11.*B22.*B33) / 2;
+
+% In exact arithmetic for a symmetric matrix  -1 <= r <= 1
+% but computation error can leave it slightly outside this range.
+phi = acos(r) / 3;
+phi(r <= -1) = pi/3;
+phi(r >= 1) = 0;
+
+% the eigenvalues satisfy lambda3 <= lambda2 <= lambda1
+lambda1 = q + 2 * p .* cos(phi);
+lambda3 = q + 2 * p .* cos(phi + pi * (2/3));
+lambda2 = 3 * q - lambda1 - lambda3;     % since trace(A) = lambda1 + lambda2 + lambda3
+
+ind = find(lambda3 < 0);
+if ~isempty(ind)
+    lambda1(ind) = lambda1(ind) - 2*lambda3(ind);
+    lambda2(ind) = lambda2(ind) - 2*lambda3(ind);
+    lambda3(ind) = lambda3(ind) - 2*lambda3(ind);
+end
+ No newline at end of file