disp(' ');
disp(' mass_condensation.m  ver 2.1  March 13, 2010 ');
disp(' by Tom Irvine ');
%
disp(' ');
disp(' This program performs mass condensation (Guyan Reduction)');
disp(' on a set of mass and stiffness matrices. ');
%
clear C;
clear CC;
clear dof;
clear sz;
clear k;
clear m;
clear mass;
clear stiff;
clear fn;
clear S2;
clear M2;
clear ModeShapes;
clear Eigenvalues;
%
clear M11;
clear M12;
clear M21;
clear M22;
%
clear K11;
clear K12;
clear K21;
clear K22;
%
disp(' ');
disp(' Enter the units system ');
disp(' 1=English  2=metric ');
iu=input(' ');
%
disp(' Assume symmetric mass and stiffness matrices. ');
%
mass_scale=1;
%
if(iu==1)
     disp(' Select input mass unit ');
     disp('  1=lbm  2=lbf sec^2/in  ');
     imu=input(' ');
     if(imu==1)
         mass_scale=386;
     end
else
    disp(' mass unit = kg ');
end
%
if(iu==1)
    disp(' stiffness unit = lbf/in ');
else
    disp(' stiffness unit = N/m ');
end
%
%
disp(' ');
disp(' Select file input method ');
disp('   1=file preloaded into Matlab ');
disp('   2=Excel file ');
file_choice = input('');
%
disp(' ');
disp(' Mass Matrix ');
%
if(file_choice==1)
        m = input(' Enter the matrix name:  ');
end
if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%
        m = xlsread(xfile);
%
end
%
m=m/mass_scale;
%
mass=m;
%
disp(' ');
disp(' Stiffness Matrix ');
%
if(file_choice==1)
        k = input(' Enter the matrix name:  ');
end
if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%
        k = xlsread(xfile);
%
end
stiff=k;
%
%
disp(' Input Matrices ');
mass
stiff
%
sz=size(mass);
dof=(sz(1));
n=dof;
%
clear ngw;
ngw=1:n;
%
% Diagonal Pivot
%
for i=1:n-1
    max=abs(mass(i,i)/stiff(i,i));
%
    for k=i+1:n
        if( abs(mass(k,k)/stiff(k,k))>max)
            max=abs(mass(k,k)/stiff(k,k));
%
            temp=ngw(i);
            ngw(i)=ngw(k);
            ngw(k)=temp;
%
% Row switch
            for j=1:n
                temp=mass(i,j);
                mass(i,j)=mass(k,j);
                mass(k,j)=temp;
%
                temp=stiff(i,j);
                stiff(i,j)=stiff(k,j);
                stiff(k,j)=temp;
            end
% Column switch
            for j=1:n
                temp=mass(j,i);
                mass(j,i)=mass(j,k);
                mass(j,k)=temp;
%
                temp=stiff(j,i);
                stiff(j,i)=stiff(j,k);
                stiff(j,k)=temp;
            end
        end
   end
end
%
disp(' Matrices with rows and columns switched');
mass
stiff
%
disp(' ');
disp(' Enter number of degrees-of-freedom ');
disp(' to retain in reduced model ');
nr=input(' ');
if(nr>n)
    nr=n;
end
d2=nr;
nmds=n-d2;
%
for i=1:d2
    for j=1:d2
        M11(i,j)=mass(i,j);
        K11(i,j)=stiff(i,j);
    end
end
%
for i=1:nmds
    for j=1:d2
        M21(i,j)=mass(i+d2,j);
        K21(i,j)=stiff(i+d2,j);
    end
end
%
for i=1:nmds
    for j=1:nmds
        M22(i,j)=mass(d2+i,d2+j);
        K22(i,j)=stiff(d2+i,d2+j);
    end
end
%
for i=1:n
    for j=1:nmds
        M12(i,j)=mass(i,j+d2);
        K12(i,j)=stiff(i,j+d2);
    end
end
%
clear invK22;
clear QQQ;
disp(' ');
disp(' Select submatrix inverse method ');
disp(' 1=inv   2=pinv  3=backslash ');
im=input(' ');
if(im==1)
    invK22=inv(K22);
end
if(im==2)
    invK22=pinv(K22);
end
if(im==3)
    QQQ=eye(nmds,nmds);
    invK22=K22\QQQ;
end
%
CC=-invK22*K21;
for i=1:d2
    for j=1:d2
         C(i,j)=0.;
    end
    C(i,i)=1.;
end
for i=1:nmds
    for j=1:d2
         C(i+d2,j)=CC(i,j);
     end
end
%
M11
M12
M21
M22
%
K11
K12
K21
K22
%
disp(' ');
disp(' inverse of K22 ');
inv(K22)
%
C
%
disp(' ');
disp(' Reduced Matrices ');
%
S2=C'*stiff*C
M2=C'*mass*C
dof=d2;
%
%  Calculate eigenvalues and eigenvectors
%
[fn,omega,ModeShapes,MST]=Generalized_Eigen(S2,M2,1);
%
disp(' ');
disp(' node switch vector ');
ngw'
%
clear MS_switched
MS_switched=C*ModeShapes
%
clear MS_ordered;
num=n;
MS_ordered=zeros(num,dof);
%
for(i=1:num)
    for(j=1:dof)
        MS_ordered(ngw(i),j)=MS_switched(i,j);
    end    
end
disp(' ');
disp(' Modeshapes with original dof order ');
MS_ordered
%
disp(' Mass Normalized Modeshapes with original dof order ');
%
clear MSTT;
clear temp;
clear QTMQ;
clear MN_MS_ordered;
%
MSTT=MS_ordered';
temp=zeros(dof,dof);
temp=m*MS_ordered;
QTMQ=MSTT*temp;
%   
clear scale;
for(i=1:dof)
    scale(i)=1./sqrt(QTMQ(i,i));
    MN_MS_ordered(:,i) = MS_ordered(:,i)*scale(i);    
end
%
MN_MS_ordered
%
disp(' Natural Frequencies ');
disp(' No.      f(Hz)');
%
for(i=1:dof)
   out1=sprintf('%d.  %12.5g ',i,fn(i));
   disp(out1);
end