disp(' ')
disp(' arbitrary_isolators.m   ver 1.6  October 13, 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 an arbitrary number ');
disp(' of isolators.');
disp(' ');
disp(' Refer to isolators_arbitrary.pdf for a diagram. ');
disp(' ');
disp(' The equation of motion is:   M (d^2x/dt^2) + K x = 0 ');
%
clear a;
clear b;
clear c;
clear k;
clear m;
clear Eigenvalues;
clear ModeShapes;
clear omega;
clear fn;
clear lambda;
clear MST;
clear ModalMass;
clear part;
clear QTMQ;
clear r;
clear x;
clear y;
clear z;
clear m6;
clear k6;
%
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(' Enter the number of isolators ');
nis=input(' ');
%
disp(' ');
disp(' Assume uniform stiffness within a given axis. ');
disp(' The stiffness values are for individual isolator springs within the axis. ');
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 the distance of each spring mounting point to the C.G.');
disp(' The dimensions are in inches. ');
disp(' ');
disp(' Select option: ');
disp('  1=preloaded matrix ');
disp('  2=keyboard entry ');
op=input(' ');
%
if(op==1)
    clear THM;
    disp(' Enter matrix name ');
    THM=input(' ');
    x=THM(:,1);
    y=THM(:,2);
    z=THM(:,3);
else
%
    for(i=1:nis)
%
        out1=sprintf('Enter x(%d)',i);
        disp(out1)
        x(i)=input(' ');    
%
        out1=sprintf('Enter y(%d)',i);
        disp(out1)
        y(i)=input(' '); 
%
        out1=sprintf('Enter z(%d)',i);
        disp(out1)
        z(i)=input(' '); 
    end
end


ox=0.;
oy=0.;
figure(1);
plot(x,y,'s');
hold on;
plot(ox,oy,'ro');
hold off;
title(' Isolator Locations XY Plane ');
xlabel('X');
ylabel('Y');
grid on;
%
figure(2);
plot(x,z,'s');
hold on;
plot(ox,oy,'ro');
hold off;
title(' Isolator Locations XZ Plane ');
xlabel('X');
ylabel('Z');
grid on;
%
figure(3);
plot(y,z,'s');
hold on;
plot(ox,oy,'ro');
hold off;
title(' Isolator Locations YZ Plane ');
xlabel('Y');
ylabel('Z');
grid on;
%
a=x;
b=y;
c=z;
%
disp(' ');
disp(' Isolator Spring Location Summary (inch) ');
disp(' ');
disp(' isolator    X        Y        Z ');
for(i=1:nis)
    out1=sprintf('\t%d\t%7.4g\t%7.4g\t%8.4g',i,a(i),b(i),c(i));
    disp(out1)
end
%
clear ab;
clear ac;
clear bc;
%
a2=a.*a;
b2=b.*b;
c2=c.*c;
ab=a.*b;
ac=a.*c;
bc=b.*c;
%
disp(' ');
%
n=nis;
k(1,1)=n*kx;
%
k(1,2)=0;
%
k(1,3)=0;
%
k(1,4)=0;
%
k(1,5)=kx*sum(c);
%
k(1,6)=kx*sum(b);
%
k(2,2)=n*ky;
%
k(2,3)=0;
%
k(2,4)=-ky*sum(c);
%
k(2,5)=0;
%
k(2,6)=ky*sum(a);
%	
k(3,3)=n*kz;
%
k(3,4)=-kz*sum(b);
%
k(3,5)=kz*sum(a);
%
k(3,6)=0;
%
k(4,4)=ky*sum(c2)+kz*sum(b2);
%
k(4,5)=kz*sum(ab);
%
k(4,6)=-ky*sum(ac);
%
k(5,5)=kx*sum(c2) + kz*sum(a2);
%
k(5,6)=kx*sum(bc);
%
k(6,6)=ky*sum(a2) + kx*sum(b2);
%
% Symmetry
%
for(i=1:6)
    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:6)
    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;
%
MST=ModeShapes';
temp=zeros(6,6);
temp=m*ModeShapes;
QTMQ=MST*temp;
%   
for(i=1:6)
    scale(i)=1./sqrt(QTMQ(i,i));
end
%
for(i=1:6)
    ModeShapes(:,i) = ModeShapes(:,i)*scale(i);  
end
%
MST=ModeShapes';
%
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
%
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
%
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
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
clear big;
clear ff
%
ff(:,1)=fn;
%
clear fn;
%
big=[ff,MSTT,PF,EMM];
big=sortrows(big,1);
%
fn    =big(:,1);
MSTT  =big(:,2:7);
PF    =big(:,8:13);
EMM   =big(:,14:19);
%
disp(' ');
disp('  Natural Frequencies = ');
%
for(i=1:6)
    out1=sprintf(' %d.   %8.4g Hz',i,fn(i));
    disp(out1);
end
%
disp(' ');
disp('  Modes Shapes (rows represent modes) ');
disp(' ');
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(' ');
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(' ');
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;    
%
    [rd_trans_x,rd_trans_y,rd_trans_z,trans_x,trans_y,trans_z,mass_x,stiff_x,stiff_y,stiff_z]=six_dof_n_iso_frf(mass,kx,ky,kz,m6,k6,nis);
%
    power_trans_x(:,1)=trans_x(:,1);
%
    clear length;
    nnn=length(trans_x(:,1));
%
    for(i=1:nnn)
        for(j=2:7)
           power_trans_x(i,j)=(trans_x(i,j))^2; 
           power_trans_y(i,j)=(trans_y(i,j))^2;    
           power_trans_z(i,j)=(trans_z(i,j))^2;      
        end
    end
%
    disp(' ');
    disp(' *** Transmissibility Output Files *** ');
    disp(' ');
    disp('   Output File Columns:  ');
    disp('     Freq    X     Y     Z     alpha       beta        theta      ');
    disp('     (Hz)  (G/G) (G/G) (G/G)  (rad/s^2/G) (rad/s^2/G) (rad/s^2/G) ');
%
    disp(' ');
    disp('    Output Filenames:  ');
    disp('      trans_x  -  X-axis input ')
    disp('      trans_y  -  Y-axis input ')    
    disp('      trans_z  -  Z-axis input ')
%
    disp(' ');
    disp(' *** Power Transmissibility Output Files *** ');
    disp(' ');
    disp('   Output File Columns:  ');
    disp('     Freq      X         Y         Z       alpha            beta  theta ');
    disp('     (Hz)  (G^2/G^2) (G^2/G^2) (G^2/G^2) ((rad/s^2)^2/G^2)    <-   <-   ');
%
    disp(' ');
    disp('    Output Filenames: ');    
    disp('      power_trans_x  -  X-axis input ')
    disp('      power_trans_y  -  Y-axis input ')    
    disp('      power_trans_z  -  Z-axis input ')   
%
    [q]=iso_apply_sine(rd_trans_x,rd_trans_y,rd_trans_z,trans_x,trans_y,trans_z);
end