disp(' ')
disp(' circular_fixed.m   ver 1.4  February 17, 2012 ')
disp(' ')
disp(' by Tom Irvine   Email: tomirvine@aol.com ')
disp(' ');
disp(' This program calculates the natural frequencies ');
disp(' of the bending modes of a fixed circular plate ');
disp(' The solution is derived from Bessel functions. ');
%
close all hidden;
%
clear n;
clear fn;
clear a;
clear b;
clear q;
clear xo;
%
clear x;
clear y;
clear root;
clear rootn;
clear rootk;
%
clear PF;
%
clear Dc;
clear Cc;
%
tpi=2*pi;
%
disp(' ');
diam=input(' Enter diameter (inch) ');
radius=diam/2;
disp(' ');
h=input(' Enter thickness (inch) ');
%
disp(' ');
imat=input(' Enter material: \n 1=aluminum  2=steel  3=G10  4=other ');
%
if(imat==1) % Al
    rho=0.1;
	E=1.0e+07;
    mu=0.3;
end
if(imat==2) % steel
    rho=0.29;
    E=3.0e+007;
    mu=0.3;
end
if(imat==3) % G10
    rho=0.065 ;
    E=2.7e+006;
    mu=0.12;
end
if(imat ~=1 && imat ~=2 && imat ~=3)
    disp(' ');   
    E = input(' Enter the elastic modulus (lbf/in^2) ');
    disp(' ');
    rho = input(' Enter the mass density (lbm/in^3) '); 
    disp(' ');
    mu=input(' Enter Poisson ratio ');    
end
%
D=E*h^3/(12*(1.-mu^2));
%
smass=rho*(diam^2)*h*pi/4;
disp(' ');
out1=sprintf(' Structural mass = %8.4g lbm',smass);
disp(out1);
disp(' ');
disp(' Add uniformly distributed non-structural mass ?');
disp(' 1=yes  2=no ')
imass=input(' ');
%
NSM=0;
if(imass==1)
  disp(' ');
  disp(' Enter added mass(lbm) ');
  NSM=input(' ');
end
%
mass=NSM+smass;
M=mass;
%
disp(' ');
out1=sprintf(' Total mass = %8.4g lbm ',M);
disp(out1);
%
mass=mass/386;
M=M/386;
rho=rho/386;   % mass per volume
%
disp(' ');
out1=sprintf(' Plate Stiffness Factor D = %8.4g lbf in',D);
disp(out1);
%
DD=sqrt(D/(rho*h));
%
delta=0.1;
%
num=6;
num2=num^2;
nm1=num-1;
%
for(ijk=1:num2)
    root(ijk)=0;
    rootn(ijk)=0;
    rootk(ijk)=0;
    Cc(ijk)=0;
    Dc(ijk)=0;
    PF(ijk)=0;
end
%
ijk=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
num_roots=(nm1+1)^2;
%
progressbar;
%
for(n=0:nm1)
%
    clear a;
    clear b;
    clear x;
    clear y;
%
    j=1;
    delta=0.001;
    for(i=1:40000)
        xx=i*delta;
        x(i)=xx;
%        
        ya=besselj(n,xx)*besseli(n+1,xx)+besselj(n+1,xx)*besseli(n,xx);
        y(i)=ya;   
%
        if(i>2)
            if(ya*yb<=0)
                b(j)=xx;
                j=j+1;
            end
        end
        yb=ya;    
        if(j>4)
            break;
        end
%        
    end
    a=b-delta;
%
%   secant
%
    k=0;
    for i=1:max(size(a))
%
        for j=1:150
%            
            xx=a(i);
            ya=besselj(n,xx)*besseli(n+1,xx)+besselj(n+1,xx)*besseli(n,xx);
%
            xx=b(i);
            yb=besselj(n,xx)*besseli(n+1,xx)+besselj(n+1,xx)*besseli(n,xx);
%
            if(abs(b(i)-a(i))<0.00001)
                break;
            end
%
            m=(yb-ya)/(b(i)-a(i));
            bb=ya-m*a(i);
            ab=-bb/m;
%
            xx=ab;
            yab=besselj(n,xx)*besseli(n+1,xx)+besselj(n+1,xx)*besseli(n,xx);   
%
            if(ya*yab<=0)
                b(i)=ab;    
            else
                a(i)=ab;
            end
%
        end
        root(ijk)=0.25*(a(i)+b(i))^2;
        rootn(ijk)=n;
        rootk(ijk)=k;
        progressbar(ijk/num_roots);
%
%  Eigenvectors & Participation Factors
%
        Dc(ijk)=0;
        Cc(ijk)=0;
        if(root(ijk)>0.01)
            L=sqrt(root(ijk));
            Dc(ijk)=besselj(n,L)/besseli(n,L);
            A=Dc(ijk);
            beta=L/radius;
            x=beta*radius;
%
            Q=0;
            Q = quad(@circular_ss_bessel_function,0,x,[],[],n,A);  % also valid for fixed
            Cc(ijk)=L/sqrt(Q);
%
            Q=0;
            if(n==0)
                Q = quad(@circular_ss_bessel_function_PF,0,x,[],[],n,A);  % also valid for fixed
            end
            PF(ijk)=2*Cc(ijk)*Q/L^2;
%
        end
%
        if(k==nm1)
            break;
        end
        k=k+1;
        ijk=ijk+1;
%        
    end 
end
%
progressbar(1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% plot(x,y);
% grid on;
% axis([0,max(x),-10,10]);
disp(' ');
root';
%
clear RM;
clear n;
clear k;
%
sz=size(Cc);
if(sz(2)>sz(1))
    Cc=Cc';
end
sz=size(Dc);
if(sz(2)>sz(1))
    Dc=Dc';
end
sz=size(PF);
if(sz(2)>sz(1))
    PF=PF';
end
%
RM=[root' rootn' rootk' Cc Dc PF];
www=RM;
%
RM = sortrows(RM,1);
%
sz=size(root);
%
for i=max(sz):-1:1
    if(RM(i,1)<=0.001)
       RM(i,:)=[];
    end
end    
%
clear Cc;
clear Dc;
%
root=RM(:,1);
n=RM(:,2);
k=RM(:,3);
Cc=RM(:,4);
Dc=RM(:,5);
PF=RM(:,6);
%
disp(' ');
disp('  Natural ');
disp(' Frequency    ');
disp('   (Hz)    n  k       C        D        root      PF    EMM');
% 
DTerm=sqrt(D/(rho*h))/(tpi*radius^2);
clear f;
clear nn;
clear kk;
j=1;
%
sz=size(root);
NT=max(sz);
if(NT>18)
    NT=18;
end
%
for i=1:NT
        fn(i)=root(i)*DTerm;
        EMM=(PF(i))^2;
        out1=sprintf('%8.2f   %d  %d   %8.5f %9.6f %8.4f %8.4f %8.4f',fn(i),n(i),k(i),Cc(i),Dc(i),sqrt(root(i)),PF(i),EMM);
        disp(out1);
        f(j)=fn(i);
        nn(j)=n(i);
        kk(j)=k(i);
        j=j+1;
end
%
disp(' ')
disp('   n = nodal diameters ');
disp('   k = nodal circles ');
disp('  PF = participation factor ');
disp(' EMM = effective modal mass '); 
%
%*********************************************************
%
fig_num=1;
%
disp(' ');
disp(' plot mode shapes? ');
disp(' 1=yes 2=no ');
ims=input(' ');
if(ims==1)
    clear n;
    mm=90;
    mmp=mm+1;
    for i=1:4
        clear Z;
        clear r;
        clear theta;
        n=nn(i);
        lambda(i)=sqrt(root(i));
%
        L=lambda(i);
%
        for ii=1:mmp
            theta(ii)=(ii-1)*2*pi/mm;
            arg=n*theta(ii);
            for(jj=1:mmp)
                r(jj)=(jj-1)/mm;  % note that the radius is non-dimensional
                kr=L*r(jj);
                Z(ii,jj)=Cc(i)*( besselj(n,kr)- Dc(i)*besseli(n,kr))*cos(arg);
            end
        end
%        
        figure(fig_num);
        fig_num=fig_num+1;
%  
        polarplot3d(Z','PlotType','surfn','meshscale',0.5,'PolarGrid',[20 72]);
%
        out1=sprintf('Mode %d  fn=%8.4g Hz  n=%d  k=%d',i,f(i),nn(i),kk(i));
        title(out1);
%
    end
end
%
