%
disp(' ');
disp(' spring_mass_series.m  ver 1.6   January 8, 2014');
disp(' Natural Frequencies of Spring-Mass Systems in Series ');
disp(' ');
disp(' By Tom Irvine  Email: tom@vibrationdata.com');
disp(' ');
%
clear k;
clear m;
clear v;
clear pf;
clear pff;
clear emm;
clear mmm;
clear mv;
clear ks;
clear fn;
clear omega;
clear ModeShapes;
clear MST;
clear nn;
clear size;
%
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(' Enter number of masses');
n=input(' ');
%
m=zeros(n,n);
%
disp(' Enter the boundary conditions ');
disp('  1=fixed-fixed  ');
disp('  2=fixed-free   ');
disp('  3=free-free    ');
ibc=input(' ');
%
disp(' ');
disp(' Select input method ');
disp(' 1=manual  2=mass & stiffness vectors ');
imethod=input(' ');
%
if(ibc==1) % fixed-fixed
        out1=sprintf('\n The boundary nodes are 0 & %d ',n+1);
        disp(out1);
end        
%
inum=1;
%
if(ibc==2) % fixed-free
        if(n>1)
            disp(' Select numbering ');
            disp('  1=  node 1 is next to fixed boundary (typical) ');
            disp('  2=  node 1 is at free end ');
            inum=input(' ');
        end
end        
%
if(ibc==1) % fixed-fixed
        nn=n+2;
end
%
if(ibc==2) % fixed-free
        nn=n+1;
end
%
if(ibc==3) % free-free
        nn=n;    
end    
%
if(imethod==1)
    for i=1:n
        out1=sprintf('\n Enter mass at node %d ',i);
        disp(out1);
        m(i,i)=input(' ');    
    end    
%
    if(ibc==1) % fixed-fixed
        for i=1:n+1
            disp(' ');        
            disp(' Enter stiffness for spring ');
            out1=sprintf(' between nodes %d & %d ',i-1,i);
            disp(out1);
            ks(i)=input(' ');
        end
        NS=n+1;
    end
%
    if(ibc==2) % fixed-free      
        for i=1:n
            disp(' ');
            disp(' Enter stiffness for spring ');
      
            if(inum==1)
                if(i==1)
                    out1=sprintf(' between fixed boundary and node %d ',i);
                else
                    out1=sprintf(' between nodes %d & %d ',i-1,i);         
                end
            else
                if(i==n)
                    out1=sprintf(' between node %d and fixed boundary ',i);
                else
                    out1=sprintf(' between nodes %d & %d ',1,i+1); 
                end             
            end    
            disp(out1);
            ks(i)=input(' ');                
        end
        NS=n;
    end
%
    if(ibc==3) % free-free
        for i=1:n-1
            disp(' ');
            disp(' Enter stiffness for spring ');
            out1=sprintf(' between nodes %d & %d ',i,i+1);
            disp(out1);
            ks(i)=input(' ');
        end
        NS=n-1;
    end    
    k=zeros(nn,nn);     
else
    mv=input(' Enter mass vector ');
    ks=input(' Enter stiffness vector ');
    for i=1:n
        m(i,i)=mv(i);
    end
    NS=max(size(ks));   
    k=zeros(nn,nn);    
end   
%
for i=1:NS
    T=zeros(2,nn);
%   
    T(1,i)=1;
    T(2,i+1)=1; 
%
    KL=[ks(i) -ks(i);-ks(i) ks(i)];
    k=k+T'*KL*T;
end
%
temp=zeros(NS,NS);
%
% Apply BCs
%
if(ibc==1) % fixed-fixed
   temp=k;
   clear k;
   k=temp(2:nn-1,2:nn-1);
end
if(ibc==2) % fixed-free
   temp=k;
   clear k;
   if(inum==1)
      k=temp(2:nn,2:nn);
   else
      k=temp(1:nn-1,1:nn-1);       
   end
end
%
m=m/mass_scale;
m
k
disp(' ');
[fn,omega,ModeShapes,MST]=Generalized_Eigen(k,m,1);
%
disp(' ');
dof=n;
%
disp(' ');
%
for i=1:dof
  v(i)=1.;
end
%
disp('        Natural    Participation    Effective  ');
disp('Mode   Frequency      Factor        Modal Mass ');
%
    LM=MST*m*v';
    pf=LM;
    sum=0;
%    
    mmm=MST*m*ModeShapes;   
%
for i=1:dof
        pff(i)=pf(i)/mmm(i,i);
        emm(i)=pf(i)^2/mmm(i,i);
        out1 = sprintf('%d  %10.4g Hz   %10.4g   %10.4g',i,fn(i),pff(i),emm(i) );
        disp(out1)
        sum=sum+emm(i);
end
out1=sprintf('\n modal mass sum = %8.4g \n',sum);
disp(out1);
%
stiffness=k;
mass=m;
%
disp(' ');
disp(' Output arrays:   mass, stiffness, fn, ModeShapes    ');