disp('  ');
disp('  plate_corners.m   ver 1.2  October 14, 2009 ');  
disp('  by Tom Irvine  Email: tomirvine@aol.com ');
disp('  ');
disp('  This program calculates the fundamental bending frequency ');
disp('  of a plate supported at each corner. ');
disp(' ');
%
clear all; 
% 
disp(' '); 
fmin=1.0e+90; 
% 
disp(' Select material '); 
disp(' 1=aluminum  2=G10  3=other '); 
imat = input(' '); 
% 
if(imat==1)      % aluminum 
    E=1.0e+07; 
    rho=0.1; 
    mu=0.3;
end 
if(imat==2)      % G10 
    E=2.7e+06; 
    rho=0.065; 
    mu=0.12;  
end 
% 
if(imat==3) 
    disp(' '); 
    disp(' Enter elastic modulus (lbf/in^2)'); 
    E=input(' '); 
    disp(' '); 
    disp(' Enter mass density (lbm/in^3)'); 
    rho=input(' '); 
    disp(' ');
    disp(' Enter Poisson ratio'); 
    mu=input(' '); 
end 
% 
rho=rho/386.; 
%
disp(' '); 
disp(' Enter thickness (in)'); 
h=input(' '); 
% 
disp(' '); 
disp(' Enter length (in)'); 
a=input(' '); 
% 
disp(' '); 
disp(' Enter width (in)'); 
b=input(' '); 
% 
added_mass=0;
disp(' '); 
disp(' Add mass?  1=yes  2=no'); 
am=input(' '); 
if(am==1) 
   disp(' '); 
   disp(' Enter mass (lbm) '); 
   added_mass=input(' '); 
   added_mass=added_mass/386; 
end 
% 
volume = a*b*h; 
rho=rho + (added_mass/volume) 
% 
total_mass=rho*volume*386;
out1=sprintf('\n total mass = %8.4g \n',total_mass);
disp(out1);
D=E*(h^3)/(12*(1-mu^2)); 
DR=D/(rho*h/2); 
out1=sprintf(' Plate Stiffness Factor = %8.4g lbf in\n',D); 
disp(out1); 
% 
num=40; 
nnn=num;
dx=a/num; 
dy=b/num; 
pia=pi/a; 
pib=pi/b; 
% 
knum=2000; 
alpha_r=0; 
beta_r=0;
gamma_r=0;
%
fn1=1.13*sqrt(D/(rho*h))/a^2;
fn2=1.13*sqrt(D/(rho*h))/b^2;
fff=[fn1 fn2];
flow=min(fff);
fup=max(fff);
%
disp(' ');
out1=sprintf(' Expected:  %8.4g < fn < %8.4g Hz ',flow,fup);
disp(out1);
%
disp(' ');
disp(' Calculating natural frequency.... ');
%
c1=2*mu;
c2=2*(1-mu);
%
disp(' ');
for (k=1:knum)
    alpha= 0.8*rand+0.2;
    beta = 0.8*rand+0.2;
%
    if(k>10 & rand > 0.9 )
        alpha=alpha_r*(0.96+0.08*rand);
        beta =beta_r*(0.96+0.08*rand);    
    end
%
    aaa=[alpha beta];
    gamma= min(aaa)*rand;
%    
    if(abs(a-b)<a/100)
        alpha=beta;
    end
%
    sum=alpha+beta+gamma;
    alpha=alpha/sum;
    beta =beta/sum;
    gamma=gamma/sum;
%
    pia_alpha=pia*alpha;
    pib_beta =pib*beta;
    pia2_alpha=(pia^2)*alpha;
    pib2_beta =(pib^2)*beta;    
%
    V=0;
    T=0;
    for(i=0:nnn)
        x=(i*dx-(a/2))*pia;
        cx=cos(x);
        sx=sin(x);
        for(j=0:nnn)
            y=(j*dy-(b/2))*pib;
            sy=sin(y);
            cy=cos(y);
            Z=alpha*cx + beta*cy + gamma*cx*cy;
%    
            dZdx  =  -pia_alpha*(1+gamma*sy)*sx;
            d2Zdx2=  -pia2_alpha*(1+gamma*sy)*cx;   
%
            dZdy  =  -pib_beta*(1+gamma*sx)*sy;
            d2Zdy2=  -pib2_beta*(1+gamma*sx)*cy;  
%    
            d2Zdxdy=pia_alpha*pib_beta*gamma*sy*sx;
%
            V1=(d2Zdx2)^2;
            V2=(d2Zdy2)^2;
            V3=c1*(d2Zdx2)*(d2Zdy2);
            V4=c2*(d2Zdxdy)^2;
%            
            V=V+ V1+V2+V3+V4;
            T=T+Z^2;
        end
    end
    V=V*D/2;
    T=T*rho*h/2;
    omega=sqrt(V/T);
    fn=omega/(2*pi);
    if(fn<fmin && fn>= flow)
        fmin=fn;
        alpha_r=alpha;
        beta_r =beta;
        gamma_r=gamma;
        out1=sprintf('%d %8.4g Hz  alpha=%8.4g  beta=%8.4g  gamma=%8.4g ',k,fn,alpha_r,beta_r,gamma_r);
        disp(out1);
    end
end
%
clear xx;
clear yy;
clear zzr;
%
figure(1);    
zmax=3;
for(i=1:41)
        xx(i)=((i-1)*dx-(a/2));
        cx=cos(xx(i)*pia);
        sx=sin(xx(i)*pia);
        for(j=1:41)
            yy(j)=((j-1)*dy-(b/2));
            sy=sin(yy(j));
            cy=cos(yy(j)*pib);
            zzr(i,j)=alpha_r*cx + beta_r*cy + gamma_r*cx*cy;
            if(zzr>zmax)
                zmax=zzr;
            end
        end
end
clear aaa;
clear abc;
abc=[xx yy];
aaa=max(abc);
xmin=-aaa;
ymin=-aaa;
zmin=0;
xmax=aaa;
ymax=aaa;
surf(xx,yy,zzr);
axis([xmin,xmax,ymin,ymax,zmin,zmax]);