disp('  ');
disp('  fixed_free_fixed_free_plate.m   ver 1.1  August 3, 2011 ');  
disp('  by Tom Irvine  Email: tomirvine@aol.com ');
disp('  ');
disp('  This program calculates the fundamental bending frequency ');
disp('  of a fixed-free-fixed_free plate via the Rayleigh method. ');
disp(' ');
%
clear fmin;
clear P;
clear W;
clear V;
clear T;
clear x;
clear y;
clear xx;
clear yy;
clear arg;
clear E;
clear rho;
clear mu;
clear h;
clear a;
clear b;
clear dx;
clear dy;
% 
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 of fixed side (in)'); 
a=input(' '); 
% 
disp(' '); 
disp(' Enter length of free side (in)'); 
b=input(' '); 
% 
added_mass=0;
volume = a*b*h;
%
disp(' '); 
disp(' Add nonstructural mass?  1=yes  2=no'); 
am=input(' '); 
if(am==1) 
   disp(' '); 
   disp(' Enter total mass (lbm) of plate + nonstructural '); 
   total_mass=input(' ');
   rho=(total_mass/volume)/386;
else
   total_mass=rho*volume*386;
end
%
out1=sprintf('\n total mass = %8.4g lbm \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=200; 
nnn=num;
dx=b/num;   % needed switch
dy=a/num; 
pia=pi/a; 
pib=pi/b; 
% 
knum=800; 
alpha_r=0; 
beta_r=0;
%
fn1=(22.373/(2*pi))*sqrt(D/(rho*h))/b^2;
%
disp(' ');
out1=sprintf(' Expected:  fn ~ %8.4g Hz ',fn1);
disp(out1);
%
disp(' ');
disp(' Calculating natural frequency.... ');
%
L=b;
%
c1=2*mu;
c2=2*(1-mu);
%
root=4.73004;
bL=root;
beta=root/L;
beta2=beta^2;
C=(sinh(bL)+sin(bL))/(cosh(bL)-cos(bL));  
%
for(i=1:nnn+1)
   ii=i-1;
   x(i)=ii*dx;
   y(i)=(ii*dy-(a/2));
end 
%    
disp(' ');
for (k=1:knum)
    if(k<knum/10 || rand<0.5)
        alpha=  -1.0 +2*rand;
        theta=   5*rand/a;
    else
        alpha=alpha_r*(0.98+0.04*rand);
        theta=theta_r*(0.98+0.04*rand);    
    end
    theta2=theta^2;
%
    if(k==1)
        alpha=0;
    end
%
    V=0;
    T=0;   
%    
    argx=beta*x;
    argy=theta*y;
%    
    for(i=1:nnn+1)
        cosh_arg=cosh(argx(i));
        cos_arg=cos(argx(i));
        sinh_arg=sinh(argx(i));
        sin_arg=sin(argx(i));
%        
                    P=(cosh_arg-cos_arg)-C*(sinh_arg-sin_arg);     
            dPdx=beta*((sinh_arg+sin_arg)-C*(cosh_arg-cos_arg));
         d2Pdx2=beta2*((cosh_arg+cos_arg)-C*(sinh_arg+sin_arg));
%        
        for(j=1:nnn+1)
            cos_arg=cos(argy(j));
            sin_arg=sin(argy(j));
                 W=1+alpha*cos_arg;
              dWdy=-alpha*theta*sin_arg;
            d2Wdy2=-alpha*theta2*cos_arg;            
%
            Z=P*W;
%    
            dZdx  =   dPdx*W; 
            d2Zdx2= d2Pdx2*W;     
%
            dZdy  = P*dWdy; 
            d2Zdy2= P*d2Wdy2;  
%    
            d2Zdxdy=dPdx*dWdy;
%
            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)
        fmin=fn;
        alpha_r=alpha;
        theta_r =theta;
        root_r=root;
        out1=sprintf('%d %8.4g Hz  alpha=%8.4g  theta=%8.4g %8.4g',k,fn,alpha_r,theta_r);
        disp(out1);
    end
end
%
clear xx;
clear yy;
clear zzr;
%
figure(1); 
zmax=0;
zmin=0;
%
num=40; 
nnn=num;
dx=b/num;   % needed switch
dy=a/num; 
%
clear arg;
clear x;
clear y;
%
for(i=1:(num+1))
          x=((i-1)*dx-(b/2));
          yy(i)=x;
          arg=beta*(x+b/2);
          P=(cosh(arg)-cos(arg))-C*(sinh(arg)-sin(arg)); 
%          
        for(j=1:(num+1))
            y=((j-1)*dy-(a/2));
            xx(j)=y;
            arg=theta_r*y;
              W=1+alpha_r*cos(arg);
%
             zzr(i,j)=P*W;
%
            if(zzr(i,j)>zmax)
                zmax=zzr(i,j);
            end
%
            if(zzr(i,j)<zmin)
                zmin=zzr(i,j);
            end            
        end
end
%
zzr=zzr/zmax;
%
clear aaa;
clear abc;
abc=[xx yy];
aaa=max(abc);
xmin=-aaa;
ymin=-aaa;
xmax=aaa;
ymax=aaa;
surf(xx,yy,zzr);
zmax=xmax;
axis([xmin,xmax,ymin,ymax,zmin,zmax]);
out1=sprintf(' fn=%8.4g Hz',fmin);
title(out1)
%
disp(' ');
out1=sprintf(' Expected:  fn ~ %8.4g Hz   Rayleigh: fn=%8.4g Hz',fn1,fmin);
disp(out1);
%