disp(' ')
disp(' six_isolators.m   ver 1.9  October 27, 2011 ')
disp(' ')
disp(' by Tom Irvine   Email: tomirvine@aol.com ')
disp(' ');
disp(' This program finds the eigenvalues and eigenvectors for a ');
disp(' six-degree-of-freedom system mounted on six isolators.');
disp(' ');
disp(' Refer to six_isolators.pdf for a diagram. ');
disp(' ');
disp(' The equation of motion is:   M (d^2x/dt^2) + K x = 0 ');
%
close all;
%
clear k;
clear m;
clear k6;
clear m6;
clear Eigenvalues;
clear ModeShapes;
clear omega;
clear fn;
clear lambda;
clear MST;
clear ModalMass;
clear part;
clear QTMQ;
clear r;
clear sum;
%
fig_num=1;
%
tpi=2.*pi;
%
n=6;
%
disp(' ');
disp(' Enter m (lbm)');
mass=input(' ');
%
disp(' Enter Jx (lbm in^2)');
jx=input(' ');
%
disp(' Enter Jy (lbm in^2)');
jy=input(' ');
%
disp(' Enter Jz (lbm in^2)');
jz=input(' ');
%
mass=mass/386.;
jx=jx/386.;
jy=jy/386.;
jz=jz/386.;
%
m=zeros(n,n);
k=zeros(n,n);
%
m(1,1)=mass;
%
m(2,2)=mass;
%	
m(3,3)=mass;
%	
m(4,4)=jx;
%
m(5,5)=jy;
%
m(6,6)=jz;
%
disp(' ');
disp(' Note that the stiffness values are for individual springs ');
disp(' ');
disp(' Enter kx (lbf/in)');
kx=input(' ');
%
disp(' ');
disp(' Enter ky (lbf/in)');
ky=input(' ');
%
disp(' ');
disp(' Enter kz (lbf/in)');
kz=input(' ');
%
disp(' ');
disp(' Enter a1 (in)');
a1=input(' ');
%
disp(' ');
disp(' Enter a2 (in)');
a2=input(' ');
%
disp(' ');
disp(' Enter b (in)');
b=input(' ');
%
disp(' ');
disp(' Enter c1 (in)');
c1=input(' ');
%
disp(' ');
disp(' Enter c2 (in)');
c2=input(' ');
%
disp(' ');
disp(' Enter c3 (in)');
c3=input(' ');
%
disp(' ');
%
k(1,1)=6*kx;
%
k(1,2)=0;
%
k(1,3)=0;
%
k(1,4)=0;
%
k(1,5)=2*kx*(-c1+c2+c3);
%
k(1,6)=6*kx*b;
%
k(2,2)=6*ky;
%
k(2,3)=0;
%
k(2,4)=2*ky*(c1-c2-c3);
%
k(2,5)=0;
%
k(2,6)=3*ky*(-a1+a2);
%	
k(3,3)=6*kz;
%
k(3,4)=-6*kz*b;
%
k(3,5)=3*kz*(a1-a2);
%
k(3,6)=0;
%
k(4,4)= 6*kz*b^2 + 2*ky*(c1^2 + c2^2 + c3^2);
%
k(4,5)=3*kz*(-a1+a2)*b;
%
k(4,6)=ky*(-a1+a2)*(c1-c2-c3);
%
k(5,5)=2*kx*(c1^2 + c2^2 + c3^2) + 3*kz*(a1^2+a2^2);
%
k(5,6)=2*kx*(-c1+c2+c3)*b;
%
k(6,6)=6*kx*b^2 + 3*ky*(a1^2+a2^2);
%
% Symmetry
%
for(i=1:n)
    for(j=1:i-1)
        k(i,j)=k(j,i);    
    end
end
%
disp(' ');
disp(' The mass matrix is');
m
disp(' ');
disp(' The stiffness matrix is');
k
%
%  Calculate eigenvalues and eigenvectors
%
[ModeShapes,Eigenvalues]=eig(k,m);
%
k6=k;
m6=m;
%
for(i=1:n)
    lambda(i)=Eigenvalues(i,i);
    if(lambda(i)<0.)
        disp(' ');
        disp(' Warning: negative eigenvalue ');
        disp(' ');
    end
end
%
disp(' ');
disp(' Eigenvalues ');
lambda
%
omega = sqrt(lambda);
fn=omega/tpi;
%
disp(' ');
disp('  Natural Frequencies = ');
%
for(i=1:n)
    out1=sprintf(' %d.   %8.4g Hz',i,fn(i));
    disp(out1);
end
%
MST=ModeShapes';
temp=zeros(n,n);
temp=m*ModeShapes;
QTMQ=MST*temp;
%   
for(i=1:n)
    scale(i)=1./sqrt(QTMQ(i,i));
end
%
for(i=1:n)
    ModeShapes(:,i) = ModeShapes(:,i)*scale(i);  
end
%
MST=ModeShapes';
%
disp(' ');
disp('  Modes Shapes (rows represent modes) ');
disp(' ');
clear MSTT;
MSTT=MST;
maximum = max(max(abs(MST)));  
for(i=1:6)
   for(j=1:6)
       if(abs(MSTT(i,j))< (maximum/1.0e+05) )
           MSTT(i,j)=0.;
       end
   end
end
out1 = sprintf('            x           y           z         alpha       beta       theta   ');
disp(out1);
for(i=1:6)
    out1 = sprintf('%d.\t%9.3g\t%9.3g\t%9.3g\t%9.3g\t%9.3g\t%9.3g ',i,MSTT(i,1),MSTT(i,2),MSTT(i,3),MSTT(i,4),MSTT(i,5),MSTT(i,6));
    disp(out1);
end
%
disp(' ');
disp('  Participation Factors (rows represent modes) ');
disp(' ');
L=MST*m;
PF=L;
PF=PF*386.;
maximum = max(max(abs(PF)));  
for(i=1:6)
   for(j=1:6)
       if(abs(PF(i,j))< (maximum/1.0e+04) )
           PF(i,j)=0.;
       end
   end
end
out1 = sprintf('            x           y           z         alpha       beta       theta   ');
disp(out1);
for(i=1:6)
    out1 = sprintf('%d.\t%9.3g\t%9.3g\t%9.3g\t%9.3g\t%9.3g\t%9.3g ',i,PF(i,1),PF(i,2),PF(i,3),PF(i,4),PF(i,5),PF(i,6));
    disp(out1);
end
%
disp(' ');
disp('  Effective Modal Mass (rows represent modes) ');
disp(' ');
EMM=L.*L;
EMM=EMM*386.;
maximum = max(max(abs(EMM)));  
for(i=1:6)
   for(j=1:6)
       if(abs(EMM(i,j))< (maximum/1.0e+04) )
           EMM(i,j)=0.;
       end
   end
end
out1 = sprintf('            x           y           z         alpha       beta       theta   ');
disp(out1);
for(i=1:6)
    out1 = sprintf('%d.\t%9.3g\t%9.3g\t%9.3g\t%9.3g\t%9.3g\t%9.3g ',i,EMM(i,1),EMM(i,2),EMM(i,3),EMM(i,4),EMM(i,5),EMM(i,6));
    disp(out1);
end
disp(' ')
disp(' Total Modal Mass ')
disp(' ')
%
out1 = sprintf(' \t%9.3g\t%9.3g\t%9.3g\t%9.3g\t%9.3g\t%9.3g ',sum(EMM(:,1)),sum(EMM(:,2)),sum(EMM(:,3)),sum(EMM(:,4)),sum(EMM(:,5)),sum(EMM(:,6)));
disp(out1);
%
%
disp(' ');
disp(' Calculate base excitation frequency response functions? ')
disp(' 1=yes  2=no ');
ifrf=input(' ');
if(ifrf==1)
%
    clear rd_trans_x;
    clear rd_trans_y;
    clear rd_trans_z;
    clear trans_x;
    clear trans_y;
    clear trans_z;
    clear power_trans_x;
    clear power_trans_y;
    clear power_trans_z;    
%
    [damp,rd_trans_x,rd_trans_y,rd_trans_z,trans_x,trans_y,trans_z,...
        power_trans_x,power_trans_y,power_trans_z,rd_power_trans_x,rd_power_trans_y,rd_power_trans_z,fig_num]...
            =iso_frf(mass,ModeShapes,omegan,fig_num);
%
    [q]=iso_apply_sine(rd_trans_x,rd_trans_y,rd_trans_z,trans_x,trans_y,trans_z);
%    
    [q,fig_num]=...
        iso_apply_random(rd_power_trans_x,rd_power_trans_y,rd_power_trans_z,power_trans_x,power_trans_y,power_trans_z,fig_num);
%
    [fig_num]=iso_apply_hs(m,fn,damp,ModeShapes,PF,fig_num);
%
    [fig_num]=iso_apply_arbit(m,fn,damp,ModeShapes,PF,fig_num);
%
end
[x]=iso_apply_rb_accel(m6,k6);