disp(' ');
disp(' matrix_deflation_beam.m  ver 1.0  August 31, 2012 ');
disp(' ');
disp(' by Tom Irvine  Email: tom@vibrationdata.com ');
disp(' ');
%
clear u;
%
clear mass;
clear stiff;
%
clear Kct
clear Mct
%
clear Kctr;
clear Mctr;
%
clear CT;
clear CR;
%
clear Kct_temp;
%
disp(' ');
mass=input(' Enter the mass matrix name:  ');
disp(' ');
stiff=input(' Enter the stiffness matrix name:  ');
%
n=max(size(mass));
%
disp(' ');
disp(' Natural Frequencies prior to Deflation ');
%
[fn,omega,ModeShapes,MST]=Generalized_Eigen_semidefinite(stiff,mass,2);
%
disp(' ');
disp(' Eigenvalues ');
disp(' ');
for i=1:n
    out1=sprintf(' %8.4f ',omega(i)^2);
    disp(out1);
end
%
disp(' ');
disp(' Form the translational constraint matrix ');
%
clear ut;
ut=zeros(n,1);
for i=1:2:n
    ut(i)=1;
end
%
CT=zeros(n,n-1);
%
for i=1:n-1
    CT(i,i)=1;
end
%   
den=0;
for i=1:n
    den=den+ut(i)*mass(i,n);  
end
%   
for i=1:n-1
%
    num=0;
    disp(' ');
    for j=1:n
        num=num+ut(j)*mass(j,i);
    end    
%
    CT(n,i)=-num/den;
%    
end
%
CT
%
Mct=CT'*mass*CT;
Kct=CT'*stiff*CT;
%
[fn,omega,ModeShapes,MST]=Generalized_Eigen_semidefinite(Kct,Mct,2);
%
disp(' ');
disp(' Eigenvalues ');
disp(' ');
for i=1:n-1
    out1=sprintf(' %8.4f ',omega(i)^2);
    disp(out1);
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
nct=n-1;
%
Kct_temp=Kct;
%
zm=zeros(n-1,1);
%
zm(1)=1;
Kct_temp(1,:)=0;
Kct_temp(1,1)=1;
%
clear ur;
ur=Kct_temp\zm;
%
disp(' ');
disp(' Rotational rigid-body mode shape ');
ur
%
CR=zeros(nct,nct-1);
%
for i=1:nct-1
    CR(i,i)=1;
end
%    
den=0;
%
for i=1:nct
    den=den+ur(i)*Mct(i,nct);
end
%   
for i=1:nct-1
%
    num=0;
    for j=1:nct
        num=num+ur(j)*Mct(j,i);
    end    
%
    CR(nct,i)=-num/den;
%    
end
%
CR
%
Mctr=CR'*Mct*CR;
Kctr=CR'*Kct*CR;
%
[fn,omega,ModeShapes,MST]=Generalized_Eigen_semidefinite(Kctr,Mctr,2);
%
disp(' ');
disp(' Eigenvalues ');
disp(' ');
for i=1:n-2
    out1=sprintf(' %8.4f ',omega(i)^2);
    disp(out1);
end
%
ModeShapes
%
disp(' ');
disp(' Restored Eigenvalues ');
disp(' ');
CT*CR*ModeShapes