disp(' ');
disp(' mdof_reduced_CB.m   ver 1.9  March 13, 2010 ');
disp(' ');
disp(' by Tom Irvine   Email: tomirvine@aol.com ');
disp(' ');
disp(' This program solves the following equation of motion: ');
disp('    M (d^2x/dt^2) + K x = 0 ');
disp('  ');
disp(' where M and K are reduced according to the users specification. ');
disp('  ');
disp(' The discarded dof may be brought back in using the Craig-Bampton method. ');
disp(' ');
%
clear acc;
clear dis;
clear vel;
clear k;
clear m;
clear damp;
clear Eigenvalues;
clear ModeShapes;
clear omega;
clear omegan;
clear fn;
clear MST;
clear ModalMass;
clear part;
clear QTMQ;
clear r;
clear t;
clear f1;
clear f2;
clear x;
clear d;
clear v;
clear na;
clear nv;
clear nd;
clear Q;
clear mass;
clear stiffness;
clear max;
clear disp;
%
tpi=2.*pi;
%
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
%
disp(' ');
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
stiffness=k;
%
size(m);
size(k);
%
num=max(size(m));
%
disp(' ');
disp(' The mass matrix is');
m
disp(' ');
disp(' The stiffness matrix is');
k
%
disp(' ');
iret=input(' Enter the number of degrees-of-freedom to retain: ');
disp(' ');
clear ret;
for(i=1:iret)
   ret(i)=input(' Enter retained dof number: ');
end
%
clear mr;
clear kr;
kr=zeros(iret,iret);
mr=zeros(iret,iret);
%
clear M11;
clear M12;
clear M21;
clear M22;
%
clear K11;
clear K12;
clear K21;
clear K22;
%
clear ngw;
ijk=iret+1;
ngw=zeros(num,1);
ngw(1:iret)=ret;
for(i=1:num)
    iflag=0;
    for(nv=1:iret)
      if(i==ret(nv))
          iflag=1;
      end
    end
    if(iflag==0)
        ngw(ijk)=i;
        ijk=ijk+1;
    end
end
%
%%%%%%%%%%%%%%%%%%%%% M11 K11 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
ic=1;
for(i=1:num)
   iflag=0;          
   for(nv=1:iret)
      if(i==ret(nv))
         iflag=1;
         break;
      end
   end
   if(iflag==1)
      jc=1;      
      for(j=1:num)
          jflag=0;
          for(nv=1:iret)
             if(j==ret(nv))
                jflag=1;
                break;
             end   
          end
          if(jflag==1)
 %            out1=sprintf(' i=%d j=%d ic=%d jc=%d \n',i,j,ic,jc);
 %            disp(out1);             
             M11(ic,jc)=m(i,j);
             K11(ic,jc)=k(i,j);
             jc=jc+1;
          end
%         out1=sprintf(' i=%d j=%d iflag=%d jflag=%d \n',i,j,iflag,jflag);
%         disp(out1);
      end
   end
   if(iflag==1)
       ic=ic+1;
   end
end
%
%%%%%%%%%%%%%%%%%%%%% M12 K12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
ic=1;
for(i=1:num)
   iflag=0;          
   for(nv=1:iret)
      if(i==ret(nv))
         iflag=1;
         break;
      end
   end
   if(iflag==1)
      jc=1;      
      for(j=1:num)
          jflag=0;
          for(nv=1:iret)
             if(j==ret(nv))
                jflag=1;
                break;
             end   
          end
          if(jflag==0)
 %            out1=sprintf(' i=%d j=%d ic=%d jc=%d \n',i,j,ic,jc);
 %            disp(out1);             
             M12(ic,jc)=m(i,j);
             K12(ic,jc)=k(i,j);
             jc=jc+1;
          end
%         out1=sprintf(' i=%d j=%d iflag=%d jflag=%d \n',i,j,iflag,jflag);
%         disp(out1);
      end
   end
   if(iflag==1)
       ic=ic+1;
   end
end
%
%%%%%%%%%%%%%%%%%%%%% M21 K21 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
ic=1;
for(i=1:num)
   iflag=0;          
   for(nv=1:iret)
      if(i==ret(nv))
         iflag=1;
         break;
      end
   end
   if(iflag==0)
      jc=1;      
      for(j=1:num)
          jflag=0;
          for(nv=1:iret)
             if(j==ret(nv))
                jflag=1;
                break;
             end   
          end
          if(jflag==1)
 %            out1=sprintf(' i=%d j=%d ic=%d jc=%d \n',i,j,ic,jc);
 %            disp(out1);             
             M21(ic,jc)=m(i,j);
             K21(ic,jc)=k(i,j);
             jc=jc+1;
          end
%         out1=sprintf(' i=%d j=%d iflag=%d jflag=%d \n',i,j,iflag,jflag);
%         disp(out1);
      end
   end
   if(iflag==0)
       ic=ic+1;
   end
end
%
%%%%%%%%%%%%%%%%%%%%%%% M22 K22 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
ic=1;
for(i=1:num)
   iflag=0;          
   for(nv=1:iret)
      if(i==ret(nv))
         iflag=1;
         break;
      end
   end
   if(iflag==0)
      jc=1;      
      for(j=1:num)
          jflag=0;
          for(nv=1:iret)
             if(j==ret(nv))
                jflag=1;
                break;
             end   
          end
          if(jflag==0)
 %            out1=sprintf(' i=%d j=%d ic=%d jc=%d \n',i,j,ic,jc);
 %            disp(out1);             
             M22(ic,jc)=m(i,j);
             K22(ic,jc)=k(i,j);
             jc=jc+1;
          end
%         out1=sprintf(' i=%d j=%d iflag=%d jflag=%d \n',i,j,iflag,jflag);
%         disp(out1);
      end
   end
   if(iflag==0)
       ic=ic+1;
   end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% M11
%% if(iret<num)
%%    M12
%%    M21
%%    M22
%%end
%
%% K11
%% if(iret<num)
%%    K12
%%    K21
%%    K22
%% end
%
clear C;
clear CC;
%
d2=iret;
nmds=num-d2;
if(iret<num)
 %%
    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;
end
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
         if(iret<num)
            C(i+d2,j)=CC(i,j);
         else
            C(i+d2,j)=0;    
         end
     end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Calculate eigenvalues and eigenvectors
%
clear mq;
clear kq;
clear m_partition;
clear k_partition;
%
m_partition(1:iret,1:iret)        =M11;
m_partition(1:iret,iret+1:num)    =M12;
m_partition(iret+1:num,1:iret)    =M21;
m_partition(iret+1:num,iret+1:num)=M22;
%
k_partition(1:iret,1:iret)        =K11;
k_partition(1:iret,iret+1:num)    =K12;
k_partition(iret+1:num,1:iret)    =K21;
k_partition(iret+1:num,iret+1:num)=K22;
%
disp(' Partitioned Matrices ');
%
m_partition
k_partition
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  
inm=0;
disp(' ');
disp(' Restore any modes via the CB method? ');
disp(' 1=yes  2=no ');
icb=input(' ');
%
if(icb==1)
    disp(' ');
    disp(' How many modes to restore? ');
    imm=input(' ');
    if(imm>nmds)
        imm=mnds;
    end
    if(imm>0)
        clear idofm;
        for(i=1:imm)
            disp(' Enter restored dof number ');
            idof = input(' ');
            idofm(i)=idof;
        end
%
        d2p1=(d2+1);
        nnn=d2+imm;
        C(:,d2p1:nnn)=0;
%
        [fn,omega,ModeShapes,MST]=Generalized_Eigen(K22,M22,1);
%
        ijk=d2+1;
        for(i=d2p1:num)
            for(j=1:imm)
                if(idofm(j)==ngw(i))
                    idel=i-d2;
                    out1=sprintf(' i=%d j=%d idofm=%d  ngw=%d idel=%d',i,j,idofm(j),ngw(i),idel);
                    disp(out1);
                    C(d2p1:num,ijk)=ModeShapes(:,idel);
                    ijk=ijk+1;
                    break;
                end
            end
        end
%
    end
%
end
disp(' ');
disp(' Transformation Matrix ');
C
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' Transformed matrices ');
%
mq=C'*m_partition*C
kq=C'*k_partition*C
%
clear ModeShapes;
clear omega;
clear fn;
%
[fn,omega,ModeShapes,MST]=Generalized_Eigen(kq,mq,1);
%
omegan=omega;
%
clear MS_switched
MS_switched=C*ModeShapes;
%
disp(' ');
disp(' node switch vector ');
ngw
%
clear MS_ordered;
MS_ordered=zeros(num,iret+imm);
%
for(i=1:num)
    for(j=1:iret+imm)
        MS_ordered(ngw(i),j)=MS_switched(i,j);
    end    
end
disp(' ');
disp(' Modeshapes with original dof order ');
MS_ordered
%
disp(' Natural Frequencies ');
disp(' No.      f(Hz)');
%
for(i=1:dof)
   out1=sprintf('%d.  %12.5g ',i,fn(i));
   disp(out1);
end