%
disp(' spring_mass_semidefinite_series.m  ver 1.3   August 31, 2012');
disp(' Natural Frequencies of Spring-Mass Systems in Series ');
disp(' ');
disp(' By Tom Irvine  Email: tomirvine@aol.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;
%
disp(' Enter number of degrees-of-freedom ');
n=input(' ');
%
m=zeros(n,n);
k=zeros(n,n);
ks=zeros(n,1);
%
disp(' Enter mass unit:  1=kg  2=lbm ');
ium=input(' ');
disp(' ');
%
disp(' Enter stiffness unit:  1=N/m  2=lbf/ft  3=lbf/in');
ius=input(' ');
%
disp(' ');
disp(' Select input method ');
disp(' 1=manual  2=mass & stiffness vectors ');
imethod=input(' ');
%
if(imethod==1)
    for(i=1:n)
        out1=sprintf('\n Enter mass %d ',i);
        disp(out1);
        m(i,i)=input(' ');
    end    
%
    for(i=1:n-1)
        out1=sprintf('\n Enter stiffness for spring between masses %d & %d ',i,i+1);
        disp(out1);
        ks(i)=input(' ');
    end
else
    mv=input(' Enter mass vector ');
    ks=input(' Enter stiffness vector ');
    for(i=1:n)
        m(i,i)=mv(i);
    end
end
%
k(1,1)=ks(1);
for(i=2:(n-1))
   k(i-1,i)=-ks(i-1);
   k(i,i)=ks(i-1)+ks(i); 
end
k(n-1,n)=-ks(n-1);
k(n,n)=ks(n-1);
%
for i=1:n
		for j=i:n
            m(i,j)=m(i,j); 
			k(j,i)=k(i,j);
        end
end
%
morig=m;
%
if(ius==2)
   m=m/32.2;
end
if(ius==3)
   m=m/386.;
end
%
disp(' ');
[fn,omega,ModeShapes,MST]=Generalized_Eigen_semidefinite(k,m,2);
%
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);
disp(' ');
disp(' mass matrix ');
m
disp(' ');
disp(' stiffness matrix ');
k
disp(' ');
disp(' mode shapes ');
ModeShapes