%
disp(' ');
disp(' circular_plate  ver 1.7  February 28, 2012 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' This program calculates the natural frequencies ');
disp(' of the bending modes of a circular plate ');
%
close all;
%
clear f;
clear kk;
clear nn;
clear lambda;
clear lambda2;
%
[iu,D,rho,mass,radius,h,mu]=circular_plate_materials();
%
disp(' ');
disp(' Enter boundary condition:')
disp(' 1=fixed  2=simply supported  3=free');
bc=input(' ');
%
DD=sqrt(D/(rho*h));
den=2*pi*radius^2;
DD=DD/den;
%
clear f;
%
if(bc==1) % fixed
    lambda2(1)=10.22;
    lambda2(2)=21.26;
    lambda2(3)=34.88;
    lambda2(4)=39.77;
    nn(1)=0;
    kk(1)=0;
    nn(2)=1;
    kk(2)=0;
    nn(3)=2;
    kk(3)=0; 
    nn(4)=0;
    kk(4)=1;     
    alpha=0.379;
%
end
if(bc==2) % SS
    lambda2(1)=4.977;
    lambda2(2)=13.94;
    lambda2(3)=25.65;
    lambda2(4)=29.76;
    nn(1)=0;
    kk(1)=0;
    nn(2)=1;
    kk(2)=0;
    nn(3)=2;
    kk(3)=0; 
    nn(4)=0;
    kk(4)=1;       
    alpha=0.441;
%
end
if(bc==3) % free
    lambda2(1)=5.253;
    lambda2(2)=9.084;
    lambda2(3)=12.23;
    lambda2(4)=20.52;
    nn(1)=2;
    kk(1)=0;
    nn(2)=0;
    kk(2)=1;
    nn(3)=3;
    kk(3)=0; 
    nn(4)=1;
    kk(4)=1;
    alpha=0.680;
%
end
f=lambda2*DD;
%
disp(' ');
disp('  Natural ');
disp(' Frequency ');
disp('   (Hz)       n    k');
for(i=1:max(size(f)))
    out1=sprintf('%9.4g \t  %d \t%d',f(i),nn(i),kk(i));
    disp(out1);
end
%
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(lambda2(i));
%
        L=lambda(i);
%               
        if(bc==3)
%
                Jn=besselj(n,L);
               dJn=0.5*(besselj(n-1,L)-besselj(n+1,L));
             dJnp1=0.5*(besselj(n,L)  -besselj(n+2,L));
              ddJn=-dJnp1+(n/L)*dJn;
%
                In=besseli(n,L);
               dIn=0.5*(besseli(n-1,L)+besseli(n+1,L));
             dInp1=0.5*(besseli(n,L)  +besseli(n+2,L));
              ddIn= dInp1+(n/L)*dIn;
%
            aa=(mu/L);
            bb=-mu*(n/L)^2;
 %           
            num=   ddJn +aa*dJn +bb*Jn;
            den=   ddIn +aa*dIn +bb*In;
            c=num/den;
        else
            c=besselj(n,L)/besseli(n,L);
        end
%
        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)=( besselj(n,kr)- c*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]);
        if(i==1)
            out1=sprintf('Mode 1  fn=%8.4g Hz  n=%d  k=%d',f(1),nn(1),kk(1));
            title(out1);
        end
        if(i==2)
            out1=sprintf('Mode 2  fn=%8.4g Hz  n=%d  k=%d',f(2),nn(2),kk(2));
            title(out1);        
        end
        if(i==3)
            out1=sprintf('Mode 3  fn=%8.4g Hz  n=%d  k=%d',f(3),nn(3),kk(3));
            title(out1);   
        end     
        if(i==4)
            out1=sprintf('Mode 4  fn=%8.4g Hz  n=%d  k=%d',f(4),nn(4),kk(4));
            title('Mode 4');
            title(out1);  
        end 
    end
end
%
disp(' ');
disp(' plot phase diagrams? ');
disp(' 1=yes 2=no ');
iph=input(' ');
if(iph==1)
%
    clear Xout;
    clear Yout;
    clear X1;
    clear Y1;
    NOP=360;
    rc=1;
    delta = rc*(1-alpha);
    THETA=linspace(0,2*pi,NOP);
    RHO=ones(1,NOP)*rc;
    [Xout,Yout] = pol2cart(THETA,RHO);
    RHO=RHO*alpha;
    [X1,Y1] = pol2cart(THETA,RHO);
%
    D1X(1)=0;
    D1X(2)=0;
    D1Y(1)=-1;
    D1Y(2)=1;
%
    D2X(1)=-1;
    D2X(2)=1;
    D2Y(1)=0;
    D2Y(2)=0;
%
    theta=60*pi/180;
    D3X(1)=cos(theta);
    D3X(2)=-cos(theta);
    D3Y(1)=sin(theta);
    D3Y(2)=-sin(theta);
%
    D4X(1)=D3X(2);
    D4X(2)=D3X(1);
    D4Y(1)=D3Y(1);
    D4Y(2)=D3Y(2);
%
    for i=1:max(size(f))
%
        figure(fig_num);
        fig_num=fig_num+1;      
%
            hold on;
            plot(Xout,Yout);
            out1=sprintf(' Mode %d  fn=%8.4g Hz  n=%d k=%d',i,f(i),nn(i),kk(i));
            title(out1);
            if(nn(i)==1)
                plot(D1X,D1Y);
            end
            if(nn(i)==2)
                plot(D1X,D1Y);
                plot(D2X,D2Y);           
            end    
            if(nn(i)==3)
                plot(D2X,D2Y); 
                plot(D3X,D3Y);
                plot(D4X,D4Y);           
            end              
            if(kk(i)==1)
                plot(X1,Y1);               
            end  
            axis square;
            hold off;       
    end
end