%
disp(' ')
disp(' plate_4.m   ver 1.0   January 14, 2009 ');
disp(' ');
disp(' by Tom Irvine   Email: tomirvine@aol.com ');
disp(' ');
disp(' This program calculates the natural frequencies ');
disp(' of a rigid plate or rectangular frame mounted via springs ');
disp(' at each corner. ');
disp(' ');
%
clear fn;
clear omega;
clear MJ;
clear MM;
clear MG;
clear KG;
%
%*****************************************************************
%
%  Input variables
%
%  There are four degrees-of-freedom:
%  Translation in the X-axis at each corner.
%
%  The X-axis is normal to the plane of the plate.
%
L1=5;      % length (in) of the plate along the Y-axis
Ly=2.5;    % CG location (in) along the Y-axis
%
L2=5;      % length (in) of the plate along the Z-axis
Lz=2.5;    % CG location (in) along the Z-axis
%
thick=0.2  % thickness (inch)
%
m=1;       % mass (lbm)
k=250;     % stiffness (lbf/in )
           % assume constant stiffness at each of the four springs
%
%*****************************************************************
%
Jy=m*(L2^2+thick^2)/12;
Jz=m*(L1^2+thick^2)/12;
%
alpha=Ly/L1;
lambda=1-alpha;
beta=Lz/L2;
sigma=1-beta;
%
alpha2 = alpha^2;
lambda2 = lambda^2;
beta2  = beta^2;
sigma2 = sigma^2;
%
alphap=1/L1;
betap=1/L2;
%
tpi=2*pi;
%
n=4;
%
MG=zeros(n,n);
KG=zeros(n,n);
MJ=zeros(n,n);
MM=zeros(n,n);
%
MM(1,1)= sigma2*lambda2;
MM(1,2)= sigma*beta*lambda2;
MM(1,3)= alpha*beta*sigma*lambda;
MM(1,4)= alpha*lambda*sigma2;
MM(2,2)= beta2*lambda2;
MM(2,3)= MM(1,2);
MM(2,4)= MM(1,3);
MM(3,3)= alpha2*beta2;
MM(3,4)= MM(1,2);
MM(4,4)= alpha2*sigma2;
%
MJy(1,1)=  lambda2;
MJy(1,2)= -lambda2;
MJy(1,3)= -alpha*lambda;
MJy(1,4)=  alpha*lambda;
MJy(2,2)=  lambda2;
MJy(2,3)=  alpha*lambda;
MJy(2,4)= -alpha*lambda;
MJy(3,3)=  alpha2;
MJy(3,4)= -alpha2;
MJy(4,4)=  alpha2;
%
MJz(1,1)=  sigma2;
MJz(1,2)=  sigma*beta;
MJz(1,3)= -sigma*beta;
MJz(1,4)= -sigma2;
MJz(2,2)=  beta2;
MJz(2,3)= -beta2;
MJz(2,4)= -sigma*beta;
MJz(3,3)=  beta2;
MJz(3,4)=  sigma*beta;
MJz(4,4)=  sigma2;
%
MM
for(i=2:n)
    for(j=1:i-1)
        MM(i,j)=MM(j,i);
        MJy(i,j)=MJy(j,i);    
        MJz(i,j)=MJz(j,i);       
    end
end
MM
%
for(i=1:n)
    KG(i,i)=1;
end
%
MM=MM*m;
MJy=MJy*Jy*betap^2;
MJz=MJz*Jz*alphap^2;
%
MG=(MM+MJy+MJz)*(1./386.);
KG=KG*k;
%
MG
KG
%
[ModeShapes,Eigenvalues]=eig(KG,MG);
%
%  98.9 Hz
%
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;
%
clear ttt;
%
for(i=1:n-1)
    for(j=i+1:n)
        if(fn(j)<fn(i))
          ttt=fn(i);
          fn(i)=fn(j);
          fn(j)=ttt;
%
          for(k=1:n)
            ttt=ModeShapes(k,i);
            ModeShapes(k,i)=ModeShapes(k,j);     
            ModeShapes(k,j)=ttt;                          
          end
%
        end
    end
end
%
disp(' ');
disp('  Natural Frequencies = ');
%
for(i=1:n)
    out1=sprintf('   %8.4g Hz',fn(i));
    disp(out1);
end
%
MST=ModeShapes';
temp=zeros(n,n);
temp=MG*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
%
ModeShapes