disp(' ')
disp(' infinite_cylinder.m  ver 1.3   November 14, 2011');
disp(' By Tom Irvine  Email: tomirvine@aol.com ');
%
disp(' ');
disp(' Select Units:  1=English  2=metric ');
iu=input(' ');
%
disp(' ')
if(iu==1)
    h=input(' Enter thickness (inch) ');
    disp(' ')
    D=input(' Enter diameter (inch) ');
else
    h=input(' Enter thickness (m) ');
    disp(' ')
    D=input(' Enter diameter (m) ');    
end
%
R=D/2.;
k=h^2/(12*R^2);
%
%  Material
%
disp(' ');
disp(' Select material '); 
disp(' 1=aluminum  2=steel  3=other '); 
imat = input(' '); 
% 
v=0.33;
if(imat==1)      % aluminum 
    if(iu==1)
        E=1.0e+07; 
        rho=0.1;
    else
        E=6.891e+010; 
        rho=2768;        
    end
end  
% 
if(imat==2)      % steel 
    if(iu==1)
        E=3.0e+07; 
        rho=0.285; 
    else
        E=2.067e+011;
        rho=7887;
    end
end 
% 
if(imat~=1 && imat ~=2) 
    disp(' ');
    if(iu==1)
        disp(' Enter elastic modulus (lbf/in^2)');
    else
        disp(' Enter elastic modulus (N/m^2)');        
    end
    E=input(' '); 
%
    disp(' '); 
    if(iu==1)
        disp(' Enter mass density (lbm/in^3)');
    else
        disp(' Enter mass density (kg/m^3)');      
    end
    rho=input(' '); 
%
    disp(' '); 
    disp(' Enter Poisson Ratio '); 
    v=input(' ');     
end 
% 
if(iu==1)
    rho=rho/386.;
end
%
a=E/((rho*(1-v^2))*R^2);
%
disp(' ')
disp(' n    fn1(Hz)    fn2(Hz) ');
n=0;
f1=0;
OMEGA2=1+k;
omega2=OMEGA2*E/((rho*(1-v^2))*R^2);
f2=sqrt(omega2)/(2*pi);
out1=sprintf(' %d \t %11.2f \t %11.2f ',n,f1,f2);
disp(out1)
%
B=1;
%
for(i=1:21)
%    
    n=i;
    n2=n^2;
%
    b=1+n2+k*(1-n2)^2;
    c=-4*k*n2*(1-n2)^2;
%
    OM2a(i)=0.5*(b-sqrt(b^2+c));
    OM2b(i)=0.5*(b+sqrt(b^2+c));
%
    Ca(i)=-B*(n^2-OM2a(i))/n;
    Cb(i)=-B*(n^2-OM2b(i))/n;
%
    om2a(i)=OM2a(i)*a;
    om2b(i)=OM2b(i)*a;
%
    fa(i)=sqrt(om2a(i))/(2*pi);
    fb(i)=sqrt(om2b(i))/(2*pi);
%
    out1=sprintf(' %d \t %11.2f \t %11.2f ',n,fa(i),fb(i));
    disp(out1)
%
end
%
disp(' ')
%
p=1;
while(p==1)
%
    ss=1;
    n=input(' Enter mode number by n value ');
    disp(' ')
    m=input(' Enter 1=low 2=high ');
%
    N=1000;
    delta=2*pi/N;
    theta=0:delta:2*pi;
    arg=n*theta;
%
    while(ss==1)
    clear V;
    clear Wa;
    clear Wb;
    clear Ra;
    clear und;
    Wa(1:length(theta))=0.;
    Wb(1:length(theta))=0.;    
%
    V=B*sin(arg);
    R=4*abs(Ca(n));
    if(n>=1)
        Wa=R+Ca(n)*cos(arg);
        Wb=R+Cb(n)*cos(arg); 
    else
        if(m==1)  % This case needs further explanation
            V(1:length(theta))=0.; 
            Wa(1:length(theta))=0.+R;
        else
            V(1:length(theta))=0.;  
            Wb(1:length(theta))=R*1.2;        
        end
    end
%    
    Ra=sqrt(V.^2+Wa.^2);   
    Rb=sqrt(V.^2+Wb.^2);
%
    if(m==1)
        alpha=theta+atan2(V,Wa);
    else
        alpha=theta+atan2(V,Wb);        
    end
%
    und(1:length(theta))=R;
%
    scale=1;
    figure(1);
%
        if(m==1)
            Ra=Ra*scale;
            polar(alpha,Ra);
            if(n>=1)
                out5 = sprintf(' n=%d  fn= %9.4g Hz',n,fa(n));
            else
                out5 = sprintf(' n=%d  fn= %9.4g Hz',n,f1);            
            end    
        else
            Rb=Rb*scale;
            polar(theta,Rb);
            if(n>=1)  
            out5 = sprintf(' n=%d  fn= %9.4g Hz',n,fb(n));
            else
                out5 = sprintf(' n=%d  fn= %9.4g Hz',n,f2);    
            end
        end
        title(out5);
        hold on;
        polar(theta,und,'-.');
        disp(' ')
        ss=input(' Rescale mode?  1=yes  2=no ');
        if(ss==1)
            hold off;
            clear figure(1);
            scale=input(' Enter new scale factor ');
            Ca=Ca*scale;
            Cb=Cb*scale;
        end
    end
%
    disp(' ');
    p=input('Make another plot?  1=yes  2=no ');
    if(p==1)
        hold off;
    end
end