%
disp(' ')
disp(' plate_k2_base.m   ver 1.2   January 23, 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;
%
clear THM;
clear y;
clear Y;
clear H;
clear tt;
%
clear n;
clear nv;
%
clear lambda;
clear sigma;
clear alpha;
clear beta;
clear lambda2;
clear sigma2;
clear beta2;
clear alpha2;
%
%*****************************************************************
%
%  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
disp(' ');
k=input(' Enter the individual spring stiffness (lbf/in) ')
%
%*****************************************************************
%
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;
%
mp=3;
%
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
%
n=max(size(MG));
%
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
%
MST=ModeShapes';
temp=zeros(n,n);
temp=MG*ModeShapes;
QTMQ=MST*temp
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
for(i=1:4)
    r(i)=1;
end
%
part = MST*MG*r';
%
disp(' ')
disp(' Particpation Factors = ');
disp(' ');
out1=sprintf('  %8.4g   ', part(1));
out2=sprintf('  %8.4g   ', part(2));
out3=sprintf('  %8.4g   ', part(3));
out4=sprintf('  %8.4g   ', part(4));
disp(out1);
disp(out2);
disp(out3);
disp(out4);
%
ModalMass = (part.*part)/1.;
%
disp(' ');
disp(' Effective Modal Mass =');
disp(' ');
out1=sprintf('  %8.4g   ', ModalMass(1));
out2=sprintf('  %8.4g   ', ModalMass(2));
out3=sprintf('  %8.4g   ', ModalMass(3));
out4=sprintf('  %8.4g   \n', ModalMass(4));
disp(out1);
disp(out2);
disp(out3);
disp(out4);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
for(i=1:4)
    damp(i)=0.05;
    omegad(i)=omega(i)*sqrt(1-(damp(i))^2);
    omega2(i)=(omega(i))^2;
end
%
clear n;
clear nv;
clear na;
clear THM;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
for(i=1:4)
%
%  Time history with 2^n, with n as an integer is preferred to eliminate
%  error from disjointed-FFTs.
%
    out1=sprintf(' Enter base input displacement time history  %d',i);
    disp(out1)
%
    [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename);
    fid = fopen(filename,'r');
    THM = fscanf(fid,'%g %g ',[2 inf]);
%
       for(j=1:20)
           if(max(size(THM(1,:)))< (2^j))
               nnn=(2^(j-1));
               tt(1:nnn)=THM(1,1:nnn);   
                y(i,1:nnn)=THM(2,1:nnn);
               break;
           end
       end
%
end
%
num=max(size(tt));
%
clear n;
clear nv;
clear na;
%
n=zeros(4,num);
nv=zeros(4,num);
na=zeros(4,num);
%
dt=(max(tt)-min(tt))/( num - 1 );
%
disp(' ');
disp(' Fade in ');
%
UM=0.01
UL=0.04;
%
for(j=1:4)
    for(i=1:num)
        if(tt(i)<=UM)
            y(j,i)=0.; 
        end
        if(tt(i)>UM && tt(i)<UL)  
            y(j,i)=y(j,i)*(tt(i)-UM)/(UL-UM);
        end
        if(tt(i)>UL)
            break;
        end
    end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
out1=sprintf('\n dt = %8.4g sec \n',dt);
disp(out1);
%
for(j=1:mp)   % row  mp = modes participating
%
    clear h;
%
    out1=sprintf('  j = %d  ',j);
    disp(out1);
%
    dwn=-damp(j)*omega(j);
    ee=exp(dwn*tt); 
        
    ot=tt*omegad(j);    
%
    sss=sin(ot);
    ccc=cos(ot);
%
    h=ee.*sss;   
    H=fft(h);     
%
    hv=ee.*( dwn*sss +omegad(j)*ccc );
    Hv=fft(hv);
%
%%%    for(i=1:1)  % column     
%
        clear w;
%                
%       Convolution via Fourier and inverse Fourier transforms
%
        arg=MST(j,1)*y(1,:) + MST(j,2)*y(2,:) + MST(j,3)*y(3,:) + MST(j,4)*y(4,:);
        Y=fft(arg);
%
        w=real(ifft(Y.*H));   
%        
        n(j,:) = (dt/omegad(j))*(k*w);
%
        v=real(ifft(Y.*Hv));  
%
        nv(j,:)= (dt/omegad(j))*(k*v);      
%
%%%  end
    nsd(j)=n(j,1);
    nsv(j)=nv(j,1);
%
end
%%%%%%%%%%%%%%%%%%%%
%
%  Starting correction
%
for(j=1:mp)
%    
    dwn(j)=damp(j)*omega(j);
    ee=exp(-dwn(j)*tt); 
    ot=tt*omegad(j);  
%
    SC=( nsv(j) + damp(j)*omega(j)*nsd(j) )/omegad(j);
%       
    for(i=1:num)
        cc=cos(ot(i));
        ss=sin(ot(i));
%
        dc=ee(i)*( nsd(j)*cc + SC*ss );
        dee=-dwn(j)*dc;
%
        n(j,i)=n(j,i) - dc;
%        
        nv(j,i)=nv(j,i) - dee -omegad(j)*ee(i)*( -nsd(j)*ss + SC*cc );
    end
%
    arg=MST(j,1)*y(1,:) + MST(j,2)*y(2,:) + MST(j,3)*y(3,:) + MST(j,4)*y(4,:);
    na(j,:)= -2*dwn(j)*nv(j,:)-omega2(j)*n(j,:)+k*arg; 
%
end    
%%%%%%%%%%%%%%%%%%%%
%
disp(' Physical Coordinates ');
%
MST*(MG*ModeShapes);
%
clear X;
clear V;
clear A;
%
X=zeros(4,num);
V=zeros(4,num);
A=zeros(4,num);
%
for(j=1:4)
%    
    for(i=1:num)
%
        for(km=1:mp)
            X(j,i)= X(j,i) + ModeShapes(j,km)*n(km,i);
            V(j,i)= V(j,i) + ModeShapes(j,km)*nv(km,i);      
            A(j,i)= A(j,i) + ModeShapes(j,km)*na(km,i);       
        end
%
    end
end
%
clear Z;
%
for(i=1:4)
    Z(i,:)=X(i,:)-y(i,:);
end
%
Z=k*Z;
%
for(i=1:4)
    zz(i)=std(Z(i,:));
end
%
disp(' ');
disp(' Spring Transmitted Forces (lbf)');
out1=sprintf(' %8.4g  %8.4g  %8.4g  %8.4g  \n',zz(1),zz(2),zz(3),zz(4));
disp(out1);
%
%%% for(j=1:4)
%%%    V(j,:)=diff(X(j,:))/dt;
%%% end
%
% for(j=1:4)
%     A(j,:)=diff(V(j,:))/dt;
% end
%
%
% tta=tt(1:tta_size);
%
A=A/386;
%
% figure(1);
% plot(tt,X(1,:));
%
% figure(2);
% plot(tt,X(2,:));
%
% figure(3);
% plot(tt,X(3,:));
%
% figure(4);
% plot(tt,X(4,:));
%
tta=tt;
A1=[tta',A(1,:)'];
A2=[tta',A(2,:)'];
A3=[tta',A(3,:)'];
A4=[tta',A(4,:)'];
%
Z1=[tt' Z(1,:)'];
Z2=[tt' Z(2,:)'];
Z3=[tt' Z(3,:)'];
Z4=[tt' Z(4,:)'];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   CG
%
alpha=Ly/L1;
lambda=1-alpha;
beta=Lz/L2;
sigma=1-beta;
%
CGamp =  A(1,:)*(lambda*sigma) + A(2,:)*( beta*lambda) + A(3,:)*(beta*alpha) + A(4,:)*(sigma*alpha );
%
plot(tta',CGamp');
title(' CG Accel ' );
xlabel(' Time (sec)' );
ylabel(' Accel (G) ' );
%
%%%% figure(1);
%%%% hold on
%%%% plot(tt,V(1,:));
%%%% plot(tt,V(2,:));
%%%% plot(tt,V(3,:));
%%%% plot(tt,V(4,:));
%
CGA=[tta',CGamp'];
%
%
sa=STD(CGamp);
%
out1=sprintf(' CG overall level = %8.4g grms \n',sa);
disp(out1);
out1=sprintf(' Max level = %8.4g G   Min level = %8.4g G  \n',max(CGamp),min(CGamp));
disp(out1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' Write the response PSD data to a file? ')
choice = input(' 1=yes 2=no ');
%
if(choice == 1)
   disp(' ');   
   disp(' Enter the output filename ');
   filename1 = input(' ','s');
   fid1 = fopen(filename1,'w');  
   for(i=1:max(size(tta)))
      fprintf(fid1,'%11.4e %11.4e \n',tta(i), CGamp(i));    
   end
   fclose(fid1);
end    