%
%  partition_matrices_accel.m  ver 1.5  October 8, 2012
%
%  need clear statements to minimized storage memory for large systems
%
function[nff,ngw,Mwd,Mww,fn,omegan,ModeShapes,MST,part,T1,T2]=...
    partition_matrices_accel(~,nem,ea,mass,stiff)
%
num=max(size(mass));
%
nff=num-nem;
%%
%
%  Track changes
%
ijk=nem+1;
ngw=zeros(num,1);
ngw(1:nem)=ea;
for i=1:num
    iflag=0;
    for nv=1:nem
      if(i==ea(nv))
          iflag=1;
      end
    end
    if(iflag==0)
        ngw(ijk)=i;
        ijk=ijk+1;
    end
end
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Partition Mass & Stiffness Matrices
%
%%%%%%%%%%%%%%%%%%%%% Mdd Kdd %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Mdd Kdd ');
ic=1;
%
Mdd=zeros(nem,nem);
Kdd=zeros(nem,nem);
%
for i=1:num
   iflag=0;          
   for nv=1:nem
      if(i==ea(nv))
         iflag=1;
         break;
      end
   end
   if(iflag==1)
      jc=1;      
      for j=1:num
          jflag=0;
          for nv=1:nem
             if(j==ea(nv))
                jflag=1;
                break;
             end   
          end
          if(jflag==1)             
             Mdd(ic,jc)=mass(i,j);
             Kdd(ic,jc)=stiff(i,j);
             jc=jc+1;
          end
      end
   end
   if(iflag==1)
       ic=ic+1;
   end
end
%
%%%%%%%%%%%%%%%%%%%%% Mdf Kdf %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' Mdf Kdf ');
ic=1;
%
Mdf=zeros(nem,nff);
Kdf=zeros(nem,nff);
%
progressbar;
%
for i=1:num
   progressbar(i/num);
   iflag=0;          
   for nv=1:nem
      if(i==ea(nv))
         iflag=1;
         break;
      end
   end
   if(iflag==1)
      jc=1;      
      for j=1:num
          jflag=0;
          for nv=1:nem
             if(j==ea(nv))
                jflag=1;
                break;
             end   
          end
          if(jflag==0)            
             Mdf(ic,jc)=mass(i,j);
             Kdf(ic,jc)=stiff(i,j);
             jc=jc+1;
          end
      end
   end
   if(iflag==1)
       ic=ic+1;
   end
end
progressbar;
%
%%%%%%%%%%%%%%%%%%%%% Mfd Kfd %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' Mfd Kfd ');
ic=1;
%
Mfd=zeros(nff,nem);
Kfd=zeros(nff,nem);
%
progressbar;
for i=1:num
   progressbar(i/num);
   iflag=0;          
   for nv=1:nem
      if(i==ea(nv))
         iflag=1;
         break;
      end
   end
   if(iflag==0)
      jc=1;      
      for j=1:num
          jflag=0;
          for nv=1:nem
             if(j==ea(nv))
                jflag=1;
                break;
             end   
          end
          if(jflag==1)             
             Mfd(ic,jc)=mass(i,j);
             Kfd(ic,jc)=stiff(i,j);
             jc=jc+1;
          end
      end
   end
   if(iflag==0)
       ic=ic+1;
   end
end
progressbar(1)
%
%%%%%%%%%%%%%%%%%%%%%%% Mff Kff %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' Mff Kff ');
ic=1;
%
Mff=zeros(nff,nff);
Kff=zeros(nff,nff);
%
progressbar;
for i=1:num
   progressbar(i/num);
   iflag=0;          
   for nv=1:nem
      if(i==ea(nv))
         iflag=1;
         break;
      end
   end
   if(iflag==0)
      jc=1;      
      for j=1:num
          jflag=0;
          for nv=1:nem
             if(j==ea(nv))
                jflag=1;
                break;
             end   
          end
          if(jflag==0)     
             Mff(ic,jc)=mass(i,j);
             Kff(ic,jc)=stiff(i,j);
             jc=jc+1;
          end
      end
   end
   if(iflag==0)
       ic=ic+1;
   end
end
pause(1);
progressbar(1)
%
clear mass;
clear stiff;
%
I=eye(nem,nem);
T2=eye(nff,nff);
nKff=max(size(Kff));
%
disp(' invKff');
invKff=pinv(Kff);
%
T1=-invKff*Kfd;
clear invKff;
%
TT=zeros(num,num);
%
TT(1:nem,1:nem)=I;
TT(nem+1:num,1:nem)=T1;
TT(nem+1:num,nem+1:num)=T2;
%
disp(' MP');
MP=zeros(num,num);
MP(1:nem,1:nem)=Mdd;
clear Mdd;
MP(1:nem,nem+1:num)=Mdf;
clear Mdf;
MP(nem+1:num,1:nem)=Mfd;
clear Mfd;
MP(nem+1:num,nem+1:num)=Mff;
clear Mff;
%
disp(' KP');
KP=zeros(num,num);
KP(1:nem,1:nem)=Kdd;
clear Kdd;
KP(1:nem,nem+1:num)=Kdf;
clear Kdf;
KP(nem+1:num,1:nem)=Kfd;
clear Kfd;
KP(nem+1:num,nem+1:num)=Kff;
clear Kff;
%
disp(' TTt*MP*TT ');
%
MT=TT'*MP*TT;
clear MP;
%
disp(' TTt*KP*TT ');
%
KT=TT'*KP*TT;
clear KP;
clear TT;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% [fn,omega,ModeShapes,MST]=Generalized_Eigen(KT,MT,2);
%
disp(' ');
clear Mwd;
clear Mww;
clear Kww;
%
disp(' Mwd');
Mwd=MT(nem+1:num,1:nem);
disp(' Mww');
Mww=MT(nem+1:num,nem+1:num);
clear MT;
disp(' Kww');
Kww=KT(nem+1:num,nem+1:num);
clear KT;
%
disp(' Eigenvalues ');
%
disp(' ');
disp(' Does the system have rigid-body modes?  1=yes 2=no ');
rigid=input(' ');
[fn,omega,ModeShapes,MST,nv]=Generalized_Eigen_fea(Kww,Mww,rigid);
fn=abs(fn);
omega=abs(omega);
ModeShapes=abs(ModeShapes);
MST=real(MST);
%
omegan=omega;
%
clear r;
clear part;
clear ModalMass;
%
sz=size(Mww);
r=ones(sz(2),1);
%
part = MST*Mww*r;
%
% disp(' Participation Factors  ');
% part
%
n=length(fn);
%
if(n>30)
    n=30;
end
%
disp(' n  fn(Hz) ');
%
for i=1:n
    out1=sprintf(' %d  %8.4g',i,fn(i));
    disp(out1);
end