disp(' ')
disp(' arbitrary_isolators.m   ver 2.6  December 19, 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(' ');
%
close all;
%
clear a;
clear b;
clear c;
clear k;
clear m;
clear Eigenvalues;
clear ModeShapes;
clear omega;
clear omegan;
clear fn;
clear lambda;
clear MST;
clear ModalMass;
clear part;
clear QTMQ;
clear r;
clear x;
clear y;
clear z;
clear m6;
clear k6;
clear sum;
%
fig_num=1;
%
tpi=2.*pi;
%
n=6;
%
disp(' ');
disp(' Enter the number of isolators ');
niso=input(' ');
nis=niso;
%
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(' 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(fig_num);
fig_num=fig_num+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(fig_num);
fig_num=fig_num+1;
plot(x,z,'s');
hold on;
plot(ox,oy,'ro');
hold off;
title(' Isolator Locations XZ Plane ');
xlabel('X');
ylabel('Z');
grid on;
%
figure(fig_num);
fig_num=fig_num+1;
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
omegan=2*pi*fn;
%
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;
    clear x;
%
    [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);