function nPos= OrthoEig(Q)
% Q is an n-by-n real orthogonal matrix
% nPos is the number of Q's eigenvalues that are strictly to the right of the
% imaginary axis.

[n,n] = size(Q);
% Tridiagnalize Q + Q' ...
S = Q + Q';
T = hess(S);
% Use Sturm sequence property...

nPos = n - Sturm(T,0);
npos = sum(eig(Q+Q')>0);

function s = Sturm(T,mu)
% T is an n-by-n symmetric tridiagonal matrix with nonzero subdiagonal
% elements.
% mu is a real number.
% s is the number of T's eigenvalues <= mu

[n,n] = size(T);
s = 0; pkm2 = 1; pkm1 = T(1,1)-mu;
if pkm1<=0
s = s+1;
end
for k=2:n
pk = (T(k,k)-mu)*pkm1 - T(k,k-1)^2*pkm2;
if sign(pk) ~= sign(pkm1)
s = s+1;
end
pkm2 = pkm1; pkm1 = pk;
end