function [eigr, lambr, eigrinv] = realeig(am)
% Function to reformulate an eigenvalue decomposition of a real-valued
% matrix provide by 'am' in the form of a near-diagonalization as 
%    an =  eigr * lambr * eigrinv
% You may freely download, copy, use, or apply this function as you wish.
%
[nr,nc] = size(am)
eigr=zeros(nr,nr);
eigrinv=zeros(nr,nr);
lambr=zeros(nr,nr);
if  nr != nc
  disp('matrix is not square')
  return
end
[eigvc,eigvl] = eig(am);

jj=0;
for i=1:nc
  jj=jj+1;
  if imag(eigvl(jj,jj)) < 1.0e-11
    eigr(:,jj) = eigvc(:,jj);
    lambr(jj,jj) = eigvl(jj,jj);
  else
    eigr(:,jj) = real(eigvc(:,jj));
    eigr(:,jj+1) = imag(eigvc(:,jj));
    lambr(jj,jj) = real(eigvl(jj,jj));
    lambr(jj+1,jj+1) = real(eigvl(jj,jj));
    lambr(jj,jj+1) = imag(eigvl(jj,jj));
    lambr(jj+1,jj) = -imag(eigvl(jj,jj));
    jj = jj+1;
  end
  if jj>=nc
    break;
  end
end
eigrinv = eigr^-1;
return

