%
disp(' rod_torsional_fn.m   ver 1.3  March 6, 2010 ');
disp(' by Tom Irvine ');
%
clear fn;
clear A;
clear J;
clear G;
clear L;
clear omegan;
%
disp(' ');
disp(' Selected geometry ');
disp('  1=solid circular  ');
disp('  2=pipe');
disp('  3=solid elliptical section ');
%
igeo=input(' ');
%
if(igeo==1)
    disp(' ');
    disp(' Enter diameter ');
    diam=input(' ');
    radius=diam/2;
    A=pi*diam^2/4; 
    Iyy=(pi/64)*(diam^4);
    Izz=Iyy;    
    J=Iyy+Izz;  % perpendicular axis theorem    
end
%
if(igeo==2)
    disp(' ');
    disp(' Enter outer diameter ');
    diam_out=input(' ');
    b=diam_out/2;
    disp(' ');
    disp(' Enter inside diameter ');
    diam_in=input(' ');
    a=diam_in/2;  
%    
    A=pi*(diam_out^2-diam_in^2)/4; 
    Iyy=(pi/64)*((diam_out^4)-(diam_in^4));
    Izz=Iyy;
    J=Iyy+Izz;  % perpendicular axis theorem  
end
%
if(igeo==3)
    disp(' ');
    disp(' Enter major diameter ');
    diam_out=input(' ');
    b=diam_out/2;
    disp(' ');
    disp(' Enter minor diameter ');
    diam_in=input(' ');
    a=diam_in/2;  
%    
    J=pi*(a^3*b^3)/(a^2 + b^2); 
end
%
disp(' Enter length ');
L=input(' ');
%
disp('  ');
disp(' Assume uniform material ');
disp('  ');
imat=input(' Enter material: \n 1=aluminum  2=steel  3=other ');
%
if(imat==1) % aluminum
    rho=0.1;
    E=1.0e+07;
    mu=0.33;
end
if(imat==2) % steel
    rho=0.29;
    E=3.0e+007;
    mu=0.33;
end
%
if(imat ~=1 && imat ~=2)
    disp(' ');   
    E   = input(' Enter the elastic modulus (lbf/in^2) ');
    disp(' ');
    rho = input(' Enter the mass density (lbm/in^3) ');
    disp(' '); 
    mu  = input(' Enter the Poisson ratio ');
end

%
G=E/(2*(1+mu));
rho=rho/386.;
%
Jm=J*(rho*L);
lambda=pi;
%
kr=G*J/L;
%
disp(' ');
disp(' Select boundary condition ');
disp(' 1=fixed-fixed ');
disp(' 2=fixed-free ');
disp(' 3=free-free ');
ibc=input(' ');
%
term=sqrt(kr/Jm);
%
if(ibc==1 || ibc==3)
   omegan(1)=1;
   omegan(2)=2;
   omegan(3)=3;
end
if(ibc==2)
   omegan(1)=1/2;
   omegan(2)=3/2;
   omegan(3)=5/2;
end
omegan=omegan*pi*term;
%
disp(' ');
out1=sprintf('   L=%8.4g in    A=%8.4g in^2   \n',L,A);
out2=sprintf(' kr=%8.4g lbf/in  \n J=%8.4g in^4 \n Jm=%8.4g lbf sec^2/in  \n',kr,J,Jm);
disp(out1); 
disp(out2);
%
fn=omegan/(2*pi);
disp(' ');
disp('    Natural ')
disp(' Frequency Hz) ');
out1=sprintf(' %9.5g \n %9.5g \n %9.5g \n',fn(1),fn(2),fn(3));
disp(out1);