freq_damp_FEM(K,M,a,b)

function [wn,z] = freq_damp_FEM(K,M,a,b)

K = full(K); M = full(M);
n     = max( size(K));
C     = a*M + b*K;   		         

A = [ zeros(n) eye(n); -M\K -M\C];

[V,D] = eig(A);            
ReD   = (D+D')/2;
ImD   = (D'-D)*i/2;

wni   = (ReD^2+ImD^2).^.5;         
zi    = -ReD/wni;                  

for k = 1:n
   z(k)  = zi(k*2,k*2);            
   wns(k) = wni(k*2,k*2);
end;
[wn,idx] = sort(abs(wns));      	 
wn    = abs(wns(idx));
zf    = z(idx);

disp('--------------------------------------------');
disp('Damped natural frequencies and modal damping');
disp('    V, Hz       T, sec     Damping');
disp('--------------------------------------------');
for j = 1:n
     BB=sprintf(' %10.5f %10.5f %10.5f \n',wn(j)/(2*pi),(2*pi)/wn(j),zf(j));
     disp(BB);
end;