disp(' ');
disp(' ss_plate_rotary_shear.m  ver 1.1  September 21, 2011 ');
disp(' by Tom Irvine ');
disp(' ')
disp(' This script calculates the natural frequencies of a ');
disp(' rectangular plate simply-supported on all four sides. ');
disp(' The script includes the effects of rotary inertia '); 
disp(' and shear deformation. '); 
disp(' ');
disp(' Reference:  Herbert, Reismann & Yu-Chung Lee ');
disp(' ');
%
disp(' Enter units:  1=English  2=metric ');
iu=input(' ');
%
if(iu==1)
    disp(' Enter length (inch) ');
else
    disp(' Enter length (m) ');   
end
L=input(' ');
%
if(iu==1)
    disp(' Enter width (inch) ');
else
    disp(' Enter width (m) ');   
end
W=input(' ');
%
if(iu==1)
    disp(' Enter thickness (inch) ');
else
    disp(' Enter thickness (m) ');   
end
h=input(' ');
%
disp(' ');
imat=input(' Enter material: \n 1=aluminum  2=steel  3=G10  4=other ');
%
    if(imat==1) % Al
        if(iu==1)
            rho=0.1;
            E=1.0e+07;
        else
            rho=2768;
            E=6.891e+010;
        end
        mu=0.3;
    end
    if(imat==2) % steel
        if(iu==1)
            rho=0.29;
            E=3.0e+007;
        else
            rho=7887;
            E=2.067e+011;
        end
        mu=0.3;
    end
    if(imat==3) % G10
        if(iu==1)
            rho=0.065;
            E=2.7e+006;
        else
            rho=1799;
            E=1.861e+010;        
        end
        mu=0.12;
    end
    if(imat ~=1 && imat ~=2 && imat ~=3)
        if(iu==1)
            disp(' ');   
            E = input(' Enter the elastic modulus (lbf/in^2) ');
            disp(' ');
            rho = input(' Enter the mass density (lbm/in^3) '); 
        else
            disp(' ');   
            E = input(' Enter the elastic modulus (N/m^2) ');
            disp(' ');
            rho = input(' Enter the mass density (kg/m^3) ');        
        end
        disp(' ');
        mu=input(' Enter Poisson ratio ');    
    end
%
    if(iu==1)
        rho=rho/386;
    end
%
disp(' ');
disp(' Enter shear coefficient (typically 0.833) ');
kappa2=input(' ');
%
c=sqrt(E/rho);
D=E*h^3/(12*(1.-mu^2));
%
a=W;
alpha=(h/a)/sqrt(12);
beta=(L/W);
%
epsilon=1/sqrt(2*(1+mu));
eta=1/sqrt(2*(1-mu));
gamma=1/sqrt(1-mu^2);
k=sqrt( kappa2/(2*(1+mu)) );
%
disp(' ');
disp(' Enter the total mode number index per axis ');
num_modes=input(' ');
%
m_max=num_modes;
n_max=num_modes;
%
A_alpha=zeros(m_max,n_max,3);
B=zeros(m_max,n_max,3);
C=zeros(m_max,n_max,3);
fn=zeros(m_max,n_max,3);
%
for(m=1:m_max)
    for(n=1:n_max)
%
        clear stiff;
        clear mass;
        clear Eigenvalues;
        clear ModeShapes;
        clear DL_omega;
%        
        out1=sprintf('\n\n m=%d n=%d \n',m,n);
        disp(out1);
%
        mass=eye(3,3);
%
        stiff(1,1)=k^2*pi^2*(m^2+(n^2/beta^2));
        stiff(1,2)=m*pi*k^2/alpha;
        stiff(1,3)=(n*pi/beta)*(k^2/alpha);
%
        stiff(2,2)=(gamma*pi*m)^2+(epsilon*pi*n/beta)^2+(k/alpha)^2;
        stiff(2,3)=eta^2*pi^2*m*n/beta;
%
        stiff(3,3)=(epsilon*pi*m)^2+(gamma*pi*n/beta)^2+(k/alpha)^2;
%
        for(i=1:3)
            for(j=i+1:3)
                stiff(j,i)=stiff(i,j);
            end
        end
%
        [ModeShapes,Eigenvalues]=eig(stiff,mass);
%
        DL_omega=zeros(3,1); 
%
        for(i=1:3)
            DL_omega(i)=sqrt(Eigenvalues(i,i));
        end
        omega=DL_omega*c/a;
        ffn=omega/(2*pi);
        ffn=sort(ffn);
        fn(m,n,:)=ffn(:);
%
        disp(' r     fn(Hz)       A           B           C    ');

        for(r=1:3)
%
%           Find normalized mode shape coefficients
%
            A_alpha(m,n,r)=1;
%            
            aaa=stiff;
            
            aaa(1,1)=aaa(1,1)-Eigenvalues(r,r);
            aaa(2,2)=aaa(2,2)-Eigenvalues(r,r);
            aaa(3,3)=aaa(3,3)-Eigenvalues(r,r);
%
            aaa(1,1)=1;
            aaa(1,2)=0;
            aaa(1,3)=0;
%
            bbb=zeros(3,1);
            bbb(1)=1;
%            
            xxx=aaa\bbb;
%
            term=(alpha^2*beta/4)*(1+(xxx(2))^2+(xxx(3))^2);
            A(m,n,r)=1/sqrt(term);
            B(m,n,r)=xxx(2)*A(m,n,r)/alpha;
            C(m,n,r)=xxx(3)*A(m,n,r)/alpha;
%
            out1=sprintf(' %d   %7.3e  %7.3e  %7.3e  %7.3e ',r,fn(m,n,r),A(m,n,r),B(m,n,r),C(m,n,r));
            disp(out1);
        end
%
    end
end
%
disp(' ');
disp(' Perform impulse force analysis?  ');
disp(' 1=yes  2=no ');
nfa=input(' ');
%
if(nfa==1)
%
    clear tt;
    clear Q;
    clear q;
    clear omega;
    clear X;
    clear Y;
    clear w;
    clear psi_x;
    clear psi_y;
%
    omega=2*pi*fn; 
%
    Q=zeros(m_max,n_max,3);
%
    disp(' ');
    disp(' Select impact location ');
    disp(' 1=center 2=elsewhere');
    iil=input(' ');
%
    if(iil==1)
        xa=L/2;
        ya=L/2;
    else
        if(iu==1)
            disp(' Enter impact location X-axis (inch) ');
        else
            disp(' Enter impact location X-axis (m) ');   
        end
        xa=input(' ');
%    
        if(iu==1)
            disp(' Enter impact location Y-axis (inch) ');
        else
            disp(' Enter impact location Y-axis (m) ');   
        end
        ya=input(' ');
    end
%
%
    disp(' ');
    disp(' Select response location ');
    disp(' 1=center 2=elsewhere');
    iir=input(' ');
%
    if(iir==1)
        xr=L/2;
        yr=W/2;
    else
        if(iu==1)
            disp(' Enter response location X-axis (inch) ');
        else
            disp(' Enter response location X-axis (m) ');   
        end
        xr=input(' ');
%    
        if(iu==1)
            disp(' Enter response location Y-axis (inch) ');
        else
            disp(' Enter response location Y-axis (m) ');   
        end
        yr=input(' ');
    end
%
    xa=xa/a;
    ya=ya/a;
    xr=xr/a;
    yr=yr/a;    
%
    if(iu==1)
        disp(' Enter total impulse (lbf sec) ');
    else
        disp(' Enter total impulse (N sec) ');   
    end
    Imp=input(' ');
%
    for(m=1:m_max)
        for(n=1:n_max)
            for(r=1:3)
                Q(m,n,r)=A(m,n,r)*sin(m*pi*xa/beta)*sin(n*pi*ya);  % switch beta from paper
            end
        end
    end
    Q=Q*Imp;
%
    dur=2/fn(1,1,1);
    dt=(1/fn(m_max,n_max,3))/10;
    num_init=round(dur/dt);
%
    nnn=round(log(num_init)/log(2));
    num=2^nnn;
%
    tt=linspace(0,dt*(num-1),num);
%
    q=zeros(m_max,n_max,3,num);
%
    for(m=1:m_max)
        for(n=1:n_max)
            for(r=1:3)
%
               s=sin(omega(m,n,r)*tt);
               X=Q(m,n,r,:);
               Y=s;
               q(m,n,r,:) = ifft((fft(X)/num).*(fft(Y)/num));              
%                
            end
        end
    end  
%     
%
    w=zeros(num,1);
    psi_x=zeros(num,1);
    psi_y=zeros(num,1);   
%
    for(m=1:m_max)
        for(n=1:n_max)
%
            arg1=m*pi*xr;
            arg2=n*pi*yr/beta;
            c1=cos(arg1);
            s1=sin(arg1);
            c2=cos(arg2);
            s2=sin(arg2);
%                
            for(r=1:3)
%
                 WW=A(m,n,r)*s1*s2;
                PPX=B(m,n,r)*c1*s2;
                PPY=C(m,n,r)*s1*c2;
%
                for(i=1:num)
%
                        w(i)=     w(i) +  WW*q(m,n,r,i);
                    psi_x(i)= psi_x(i) + PPX*q(m,n,r,i);
                    psi_y(i)= psi_y(i) + PPY*q(m,n,r,i); 
%                
                end
            end
        end
    end  
%
    wmax=max(w);
    wmin=min(w);
%
    if(iu==1)
        out1=sprintf('\n max disp = %8.3e inch ',wmax);
        out2=sprintf(' min disp = %8.3e inch ',wmin);
    else
        out1=sprintf('\n max disp = %8.3e m ',wmax);
        out2=sprintf(' min disp = %8.3e m ',wmin);        
    end
%    
    disp(out1);
    disp(out2);
%
    sz=size(tt);
    if(sz(2)>sz(1))
        tt=tt';
    end
    sz=size(w);
    if(sz(2)>sz(1))
        w=w';
    end
    sz=size(psi_x);
    if(sz(2)>sz(1))
        psi_x=psi_x';
    end    
    sz=size(psi_y);
    if(sz(2)>sz(1))
        psi_y=psi_y';
    end        
%
    figure(1);
    plot(tt,w);
    grid on;
    xlabel('Time (sec) ');
    if(iu==1)
        ylabel('Disp (in)');
    else
        ylabel('Disp (m)');   
    end    
%    
    figure(2);
    plot(tt,psi_x);
    grid on;
    title('Rotation about X-axis ');
    xlabel('Time (sec) ');
    ylabel('Rotation (rad)');
%    
    figure(3);
    plot(tt,psi_y);
    grid on;
    title('Rotation about Y-axis ');
    xlabel('Time (sec) ');
    ylabel('Rotation (rad)');
%
end