disp(' ');
disp('    two_dof_force_frf_alt.m   ver 1.5  June 11, 2010 ');
disp(' ');
disp('    by Tom Irvine   Email: tomirvine@aol.com ');
disp(' ');
disp(' This program finds the frequency response of a ');
disp(' two-degree-of-freedom system subjected to applied forces. ');
disp(' The equation of motion is:   M (d^2x/dt^2) + K x = F(t) ');
disp(' ');
clear k;
clear m;
clear Eigenvalues;
clear ModeShapes;
clear omega;
clear fn;
clear MS;
clear MST;
clear ModalMass;
clear part;
clear QTMQ;
clear r;
clear H1;
clear H2;
%
clear HR1;
clear HR2;
%
tpi=2.*pi;
%
n=2;
%
disp(' Assume symmetric mass and stiffness matrices. ');
disp(' ');
%
disp(' Select units: 1=English  2=metric ');
iunit=input(' ');
%
if(iunit==1)
   disp(' Select mass units: 1=lbm   2=lbf sec^2/in ');
   imass=input(' ');  
end
%
if(iunit==1)
    if(imass==1)
       mass_unit='lbm';
    else
       mass_unit='lbf sec^2/in';        
    end
    stiffness_unit='lbf/in';
    force_unit='lbf';
else
    mass_unit='kg';
    stiffness_unit='N/m';
    force_unit='N';
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    
out1=sprintf(' Enter m11 (%s)',mass_unit);
disp(out1);
m(1,1)=input(' ');
%
out1=sprintf(' Enter m22 (%s)',mass_unit);
disp(out1);	
m(2,2)=input(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
m(1,2)=0;
m(2,1)=m(1,2);
%
if(iunit==1 && imass==1)
    m=m/386;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
out1=sprintf(' Enter k11 (%s)',stiffness_unit);
disp(out1);
k(1,1)=input(' ');
%
out1=sprintf(' Enter k12 (%s)',stiffness_unit);
disp(out1);
k(1,2)=input(' ');
%	
out1=sprintf(' Enter k22 (%s)',stiffness_unit);
disp(out1);
k(2,2)=input(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
k(2,1)=k(1,2);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
out1=sprintf(' Enter the force applied to dof 1 (%s)',force_unit);
disp(out1);
B(1)=input(' ');
%
out1=sprintf(' Enter the frequency (Hz) of dof 1 force');
disp(out1);
alpha(1)=input(' ');
%
out1=sprintf(' Enter the phase (deg) of dof 1 force');
disp(out1);
phase(1)=input(' ');
%
out1=sprintf(' Enter the force applied to dof 2 (%s)',force_unit);
disp(out1);
B(2)=input(' ');
%
out1=sprintf(' Enter the frequency (Hz) of dof 2 force');
disp(out1);
alpha(2)=input(' ');
%
out1=sprintf(' Enter the phase (deg) of dof 2 force');
disp(out1);
phase(2)=input(' ');
%
phase=phase*pi/180;
alpha=alpha*2*pi;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Enter modal damping ratio 1 ');
damp(1)=input(' ');
%
disp(' Enter modal damping ratio 2 ');
damp(2)=input(' ');
%
disp(' ');
disp(' The mass matrix is');
m
disp(' ');
disp(' The stiffness matrix is');
k
%
%  Calculate eigenvalues and eigenvectors
%
[fn,omegan,ModeShapes,MST]=Generalized_Eigen(k,m,1);
%
for(i=1:2)
    domegan(i)=damp(i)*omegan(i);
end
%
clear C;
C=zeros(2,1);
for(i=1:2)
    for(j=1:2)
        p(1,1)=   omegan(i)^2-alpha(j)^2;
        p(1,2)=-2*domegan(i)*alpha(j);
        p(2,1)=-p(1,2);
        p(2,2)= p(1,1);
%
        MM=ModeShapes(j,i)*B(j);
%        
        C(1)=MM*cos(phase(j));
        C(2)=MM*sin(phase(j));
%
        D=p\C;
%
        N(i,j)=D(1);
        P(i,j)=D(2);     
    end
end
%
MN=zeros(2,2);
MP=zeros(2,2);    
%
for(i=1:2)
    for(j=1:2)
       for(n=1:2)
          MN(i,j)=MN(i,j)+ModeShapes(i,n)*N(n,j);
          MP(i,j)=MP(i,j)+ModeShapes(i,n)*P(n,j);
       end   
    end    
end
%
MT=zeros(2,1);
RMT=zeros(2,1);
for(i=1:2)
    for(j=1:2)    
       M(i,j)=sqrt(MN(i,j)^2+MP(i,j)^2);
       RMT(i)=RMT(i)+M(i,j)^2;
       MT(i)=MT(i)+M(i,j);
       phi(i,j)=atan2(MN(i,j),MP(i,j));
%       
       while(phi(i,j)<0)
          phi(i,j)=phi(i,j)+2*pi;
       end
%       
    end
    RMT(i)=sqrt(RMT(i)/2.);
end
%
disp(' ');
disp(' Phase angle convention:  arctan( sine amplitude/cosine amplitude)');
disp(' ');
disp('Displacement Results');
disp(' ');
disp(' Response  Excitation    Freq    Magnitude    Phase       ');
disp('    dof       dof        (Hz)     (inch)      (deg)       ');
%
for(i=1:2)
   for(j=1:2)  
      out1=sprintf('     %d         %d      %6.2f     %7.3e   %6.1f  ',i,j,alpha(j)/(2*pi),M(i,j),phi(i,j)*180/pi);
      disp(out1);
   end  
end     
disp(' ');
disp(' dof   RMS(inch)  Peak Magnitude(inch) ');  
for(i=1:2) 
    out1=sprintf('   %d   %8.4g    %8.4g  ',i,RMT(i),MT(i));
    disp(out1);    
end    
%
%  velocity
%
MN=zeros(2,2);
MP=zeros(2,2);  
%
for(i=1:2)
    for(j=1:2)
       for(n=1:2)
          MN(i,j)=MN(i,j)+ModeShapes(i,n)*N(n,j)*alpha(j);
          MP(i,j)=MP(i,j)+ModeShapes(i,n)*P(n,j)*alpha(j);
       end   
    end    
end
%
MT=zeros(2,1);
RMT=zeros(2,1);
for(i=1:2)
    for(j=1:2)    
       M(i,j)=sqrt(MN(i,j)^2+MP(i,j)^2);
       RMT(i)=RMT(i)+M(i,j)^2;
       MT(i)=MT(i)+M(i,j);
       phi(i,j)=atan2(-MP(i,j),MN(i,j));
%       
       while(phi(i,j)<0)
          phi(i,j)=phi(i,j)+2*pi;
       end
%       
    end
    RMT(i)=sqrt(RMT(i)/2.);
end
%
%
disp(' ');
disp('Velocity Results');
disp(' ');
disp(' Response  Excitation    Freq    Magnitude    Phase       ');
disp('    dof       dof        (Hz)    (in/sec)     (deg)       ');
%
for(i=1:2)
   for(j=1:2)  
      out1=sprintf('     %d         %d      %6.2f     %7.3e   %6.1f  ',i,j,alpha(j)/(2*pi),M(i,j),phi(i,j)*180/pi);
      disp(out1);
   end  
end     
disp(' ');
disp(' dof   RMS(in/sec)  Peak Magnitude(in/sec) ');  
for(i=1:2) 
    out1=sprintf('   %d   %8.4g    %8.4g  ',i,RMT(i),MT(i));
    disp(out1);    
end    
%
%
%  acceleration
%
MN=zeros(2,2);
MP=zeros(2,2);  
%
for(i=1:2)
    for(j=1:2)
       for(n=1:2)
          MN(i,j)=MN(i,j)+ModeShapes(i,n)*N(n,j)*(alpha(j))^2;
          MP(i,j)=MP(i,j)+ModeShapes(i,n)*P(n,j)*(alpha(j))^2;
       end   
    end    
end
%
MT=zeros(2,1);
RMT=zeros(2,1);
for(i=1:2)
    for(j=1:2)    
       M(i,j)=sqrt(MN(i,j)^2+MP(i,j)^2);
       RMT(i)=RMT(i)+M(i,j)^2;
       MT(i)=MT(i)+M(i,j);
       phi(i,j)=atan2(-MN(i,j),-MP(i,j));
%       
       while(phi(i,j)<0)
          phi(i,j)=phi(i,j)+2*pi;
       end
%       
    end
    RMT(i)=sqrt(RMT(i)/2.);
end
%
%
disp(' ');
disp('Acceleration Results');
disp(' ');
disp(' Response  Excitation    Freq    Magnitude    Phase       ');
disp('    dof       dof        (Hz)   (in/sec^2)     (deg)       ');
%
for(i=1:2)
   for(j=1:2)  
      out1=sprintf('     %d         %d      %6.2f     %7.3e   %6.1f  ',i,j,alpha(j)/(2*pi),M(i,j),phi(i,j)*180/pi);
      disp(out1);
   end  
end     
disp(' ');
disp(' dof   RMS(in/sec^2)  Peak Magnitude(in/sec^2) ');  
for(i=1:2) 
    out1=sprintf('   %d   %8.4g    %8.4g  ',i,RMT(i),MT(i));
    disp(out1);    
end    
%