function [ states, waves ] = riemann2x2( A, q_l, q_r )
% 
% Input:
% A: 2 x 2 real matrix
% q_l: 2 element matrix
% q_r: 2 element matrix
% 
% Output:
% Returns the solution:
%   3 states:
%     .left   = 1 x 2 matrix
%     .middle = 1 x 2 matrix
%     .right  = 1 x 2 matrix
%   2 waves:
%     .left  = scalar
%     .right = scalar
% 
% For the P.D.E.:
% 
% u_t + A u_x = 0
% 
  
% 
% Written by:
% -- 
% John L. Weatherwax, Ph.D.         2003-02-11
% MIT Lincoln Laboratory
% 244 Wood Street
% Lexington, MA 02420-9176 USA
% 
% email: weatherwax@ll.mit.edu
% Voice: (781) 981-5370       Fax: (781) 981-0721
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

% Insure the input is a COLUMN vector:
q_l = q_l(:);
q_r = q_r(:);

% Compute the eigen{vectors,values}:
[ R, D ] = eig( A );
% Sort the eigenvalues in increasing order:
eValues = diag( D );
[ s, i ] = sort( eValues );
% Rearraing the eigenvectors/eignvector matrix in the same order:
eValues = eValues( i ); R = R( :, i );

alpha = R \ ( q_r-q_l );

q_m = q_l + alpha(1)*R( :, 1 );

% Pass out into structure:
states.left   = q_l(:).';
states.middle = q_m(:).';
states.right  = q_r(:).';

waves.left  = eValues(1);
waves.right = eValues(2);