function [val1,val2,vec1x,vec1y,vec2x,vec2y] = eig2d(a,b,c,d)
% function [val1,val2,vec1x,vec1y,vec2x,vec2y] = eig2d(a,b,c,d)
% Calculates the eigenvalues and eigenvectors of 2x2 matrices
% using the method of Newton (Newton, T.A., 1990, A simple method
% for finding eigenvalues and eigenvectors for 2x2 matrices: 
% The American Mathematical Monthly, v. 97, p. 57-60.
% The input and output are arranged as vectors rather than
% matrices to allow the code to be used without the need for
% looping over a grid.
% I AM HAVING TROUBLE WITH [a b; c d] = [1 0 ; 2 2] !!!!
% Input parameters (arguments)
% a = vector of coefficients (n x 1)
% b = vector of coefficients (n x 1)
% c = vector of coefficients (n x 1)
% d = vector of coefficients (n x 1)
% Output parameters (arguments)
% val1 = vector of greatest eigenvalues (n x 1)
% val2 = vector of least eigenvalues (n x 1)
% vec1x = vector of x-components for eigenvalues for val1 (n x 1)
% vec1y = vector of y-components for eigenvalues for val1 (n x 1)
% vec2x = vector of x-components for eigenvalues for val2 (n x 1)
% vec2y = vector of y-components for eigenvalues for val2 (n x 1)
% Examples
% [val1,val2,vec1x,vec1y,vec2x,vec2y] = eig2d(1,2,3,4)   % This checks!
% [val1,val2,vec1x,vec1y,vec2x,vec2y] = eig2d(1,-2,-1,2) % This checks!
% [val1,val2,vec1x,vec1y,vec2x,vec2y] = eig2d(1,1,-2,3)  % Checks.  Odd.
% [val1,val2,vec1x,vec1y,vec2x,vec2y] = eig2d(3,4,4,-3)  % This checks!
% [val1,val2,vec1x,vec1y,vec2x,vec2y] = eig2d(1,0,2,2)   % This checks!
% [val1,val2,vec1x,vec1y,vec2x,vec2y] = eig2d(2,0,0,2)   % This checks!
% [val1,val2,vec1x,vec1y,vec2x,vec2y] = eig2d(2,0,6,1)   % This checks!

% Dimension output vectors
%val1 = zeros(length(a),1);
val1 = zeros(size(a));
val2 = val1;
vec1x = val1;
vec1y = val1;
vec2x = val1;
vec2y = val1;
m1 = val1;
m2 = val2;

% Three cases
% Case a of Newton
i = find(b~=0);  
C = ( a(i)+d(i) )/2;
R = sqrt( (d(i)-a(i)).^2 + 4*b(i).*c(i) )/2;
val1(i) = C + R;
val2(i) = C - R;
m1(i) = ( val1(i) - a(i) )./b(i);
m2(i) = ( val2(i) - a(i) )./b(i);
% Find the corresponding (normalized) eigenvectors
vec1x(i) = 1./sqrt(1+m1(i).^2);
vec1y(i) = m1(i).*vec1x(i);
vec2x(i) = 1./sqrt(1+m2(i).^2);
vec2y(i) = m2(i).*vec2x(i);

% Case b of Newton (two identical eigenvalues)
j = find(a==d & b==0);
val1(j) = d(j);
val2(j) = val1(j);


%vec1x(j) = zeros(j,1); % Possible vector dimension problem here
vec1x(j) = zeros(size(val1(j))); % Possible vector dimension problem here
%vec1y(j) = ones(j,1);   % Possible vector dimension problem here
vec1y(j) = ones(size(val1(j)));   % Possible vector dimension problem here
vec2x(j) = vec1x(j);  % Possible vector dimension problem here
vec2y(j) = vec1y(j);   % Possible vector dimension problem here

    % Subcase of isotropic tensor.  Eigenvectors set to be orthogonal
    j1 = find (a==d & b==0 & c==0);
    vec2x(j1) = vec1y(j1);  % Possible vector dimension problem here
    vec2y(j1) = vec1x(j1);  % Possible vector dimension problem here


% Mix of Case a and Case b of Newton (two different eigenvalues)
k1 = find(a~=d & b==0 & a>d);
k2 = find(a~=d & b==0 & d>a);

val1(k1) = a(k1);   % Maximum eigenvalue = a
val2(k1) = d(k1);   % Minimum eigenvalue = d
m = c(k1)./(a(k1)-d(k1));
vec1x(k1) = 1./sqrt(1+m.^2);
vec1y(k1) = m.*vec1x(k1);
%vec2x(k1) = 0; 
%vec2y(k1) = 1;  
vec2x(k1) = zeros(size(a(k1))); 
vec2y(k1) = ones(size(a(k1)));   

val1(k2) = d(k2);   % Maximum eigenvalue = d
val2(k2) = a(k2);   % Minimum eigenvalue = a
m = c(k2)./(a(k2)-d(k2));
%vec1x(k2) = 0;
%vec1y(k2) = 1;
vec1x(k2) = zeros(size(d(k2)));
vec1y(k2) = ones(size(d(k2)));
vec2x(k2) = 1./sqrt(1+m.^2); 
vec2y(k2) = m.*vec2x(k2);