disp(' ');
disp(' SD_SD_cylinder.m   ver 1.0   November 27, 2012');
disp('  ');
disp(' by Tom Irvine   Email: tom@vibrationdata.com ');
disp(' ');
%
disp(' This program calculates the natural frequencies of a ');
disp(' cylindrical shell which is enclosed at each end by a ');
disp(' "shear diaphragm." ');
%
disp(' ');
disp(' Reference:  Leissa, Vibration of Shells, NASA SP-288 ');
%
close all hidden;
%
fig_num=1;
%
disp(' ');
disp(' Select units:  1=English  2=metric ');
iu=input(' ');
%
[E,rho,mu]=materials(iu);
%
if(iu==1)
     diam=input(' Enter the diameter (inch)       ');
        L=input(' Enter the length (inch)         ');
    thick=input(' Enter the wall thickness (inch) ');
else
     diam=input(' Enter the diameter (m)        ');
        L=input(' Enter the length (m)          ');
    thick=input(' Enter the wall thickness (mm) ');
    thick=thick/1000;
end
%
R=diam/2;
h=thick;
%
k=h^2/(12*R^2);
%
M=4;
lambda=zeros(M,1);
L2=zeros(M,1);
L4=zeros(M,1);
%
for m=1:M
    lambda(m)=m*pi*R/L;
    L2(m)=lambda(m)^2;
    L4(m)=lambda(m)^4;    
end    
%
hdR=h/R;
%
out1=sprintf('\n  h/R = %8.4g ',hdR);
disp(out1);
%
a=(1-mu^2)*R^2;
%
n=0;
%
disp(' ');
disp(' Torsional Modes (Love-Timoshenko)');
disp('  m   n  freq (Hz)')
%
for m=1:M
%
    OM2=0.5*(1-mu)*(1+2*k)*L2(m);
%
    om2=OM2*(E/rho)/a;
    om=sqrt(om2);
    f=om/(2*pi);
    out1=sprintf('   %d  %d %8.4g ',m,n,f);
    disp(out1);
end
%
%
disp(' ');
disp(' Torsional Modes (Donnell-Mushtari)');
disp('  m   n  freq (Hz)')
%
for m=1:M
%
    OM2=0.5*(1-mu)*L2(m);
%
    om2=OM2*(E/rho)/a;
    om=sqrt(om2);
    f=om/(2*pi);
    out1=sprintf('   %d  %d %8.4g ',m,n,f);
    disp(out1);
end
%
disp(' ');
disp(' Coupled Axial-Radial modes (Donnell-Mushtari)');
disp('           low        high')
disp('  m   n   freq(Hz)   freq(Hz)')
%
for m=1:M
%
    term1=1+L2(m)+k*L4(m);
    term2=(1-L2(m))^2 + 2*L2(m)*(2*mu^2+k*L2(m)-k*L4(m));
    term2=sqrt(term2);
%
    OM2=0.5*(term1-term2);
%
    om2=OM2*(E/rho)/a;
    om=sqrt(om2);
    flow=om/(2*pi);
%
    OM2=0.5*(term1+term2);    
%
    om2=OM2*(E/rho)/a;
    om=sqrt(om2);
    fhigh=om/(2*pi);
%
    out1=sprintf('   %d  %d  %8.4g  %8.4g',m,n,flow,fhigh);
    disp(out1);
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
n=2;
%
% 
% hdR=0.100
yn2_4 =@(x) -4.14E-4*x^8 +0.00777*x^7 -0.0601*x^6 +0.246*x^5 -...
              0.554*x^4 +0.626*x^3 -0.178*x^2 +0.0383*x +0.0754;
%
% hdR=0.074
yn2_3 =@(x) -8.24E-4*x^8 +0.0143*x^7 -0.103*x^6 +0.391*x^5 -...  
             0.826*x^4 +0.893*x^3 -0.297*x^2 +0.0589*x +0.0564;
%
% hdR=0.0525
yn2_2 =@(x) +1.77E-4*x^8 -0.00207*x^7 +0.00646*x^6 +0.0104*x^5 -...  
             0.0899*x^4 +0.115*x^3 +0.111*x^2 -0.0232*x +0.0419;
%
% hdR=0.010
yn2_1 =@(x) +0.00663*x^4 -0.0687*x^3 +0.211*x^2 +0.0158*x +0.00126;
%
H4=0.1;
H3=0.074;
H2=0.0525;
H1=0.01;
%
if(hdR<=H4)
%
    disp(' ');
    disp(' Coupled Axial-Radial modes (Arnold-Warburton)');
    disp('  m   n   freq(Hz) ')
%
    for m=1:M
%
        OM=0;
        x=lambda(m);
%
        if(x<=4.0)
%
            if(hdR==H4)
                OM=yn2_4(x);
            end
%
            if(H3< hdR && hdR<H4)
                A=[yn2_4(x) yn2_3(x)];
                OMb=max(A);
                OMa=min(A);
                c2=(hdR-H3)/(H4-H3);
                c1=1-c2;
                OM=c1*OMa+c2*OMb;
            end
%
            if(hdR==H3)
                OM=yn2_4(x);
            end
%
            if(H2< hdR && hdR<H3)
                A=[yn2_3(x) yn2_2(x)];
                OMb=max(A);
                OMa=min(A);
                c2=(hdR-H2)/(H3-H2);
                c1=1-c2;
                OM=c1*OMa+c2*OMb;
            end
%
            if(hdR==H2)
                OM=yn2_4(x);
            end
%
            if(H1< hdR && hdR<H2)
                A=[yn2_2(x) yn2_1(x)];
                OMb=max(A);
                OMa=min(A);
                c2=(hdR-H1)/(H2-H1);
                c1=1-c2;
                OM=c1*OMa+c2*OMb;
            end
%            
            if(hdR<=H1)
                OM=yn2_1(x);
            end            
%
            om=OM*sqrt((E/rho)/a);
            f=om/(2*pi);
            out1=sprintf('   %d  %d %8.4g ',m,n,f);
            disp(out1);
%    
        else
            break;
        end
    end
end
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
n=3;
% 
yn3_8 =@(x) +2.47E-5*x^6 -1.99E-5*x^5 -1.66E-4*x^4 -0.0117*x^3 +... 
             0.0981*x^2 -0.0411*x+0.217;
%
yn3_7 =@(x) -0.00139*x^6 +0.0179*x^5 -0.0857*x^4 +0.176*x^3 -... 
             0.0791*x^2 +0.0195*x+0.217;
%
yn3_6 =@(x) +6.46E-4*x^8 -0.0103*x^7 +0.0659*x^6 -0.212*x^5 +... 
             0.357*x^4 -0.304*x^3+0.183*x^2 -0.0402*x +0.162;
%
yn3_5 =@(x) +1.26E-4*x^8 -0.00144*x^7 +0.00381*x^6 +0.0163*x^5 -...
             0.117*x^4 +0.245*x^3 -0.132*x^2 +0.0241*x +0.118;
%
yn3_4 =@(x) -4.86E-4*x^8 +0.00867*x^7 -0.0638*x^6 +0.249*x^5 -... 
             0.549*x^4 +0.649*x^3-0.283*x^2 +0.0433*x +0.0598;
%
yn3_3 =@(x)-3.63E-4*x^8 +0.00646*x^7 -0.0475*x^6 +0.186*x^5 -...
           0.413*x^4 +0.488*x^3 -0.192*x^2 +0.0357*x +0.0423;
%
yn3_2 =@(x) +2.58E-4*x^8 -0.00388*x^7 +0.0227*x^6 -0.0640*x^5 +...
             0.0853*x^4 -0.0588*x^3+0.103*x^2 -0.0125*x +0.0220;
%             
yn3_1 =@(x) -5.68E-4*x^8 +0.00930*x^7 -0.0618*x^6 +0.213*x^5 -...
             0.395*x^4 +0.347*x^3-0.0138*x^2 -0.0156*x +0.00918;
%
H8=0.1;
H7=0.091;
H6=0.074;
H5=0.0525;
H4=0.0243;
H3=0.020;
H2=0.010;
H1=0.001; 
%
if(hdR<=0.100)
%
    disp(' ');
    disp(' Coupled Axial-Radial modes (Arnold-Warburton)');
    disp('  m   n   freq(Hz) ')   
%
    for m=1:M
%
        OM=0;
        x=lambda(m);
%
        if(x<=4.0)
%
            if(hdR==H8)
                OM=yn3_8(x);
            end
%
            if(H7< hdR && hdR<H8)
                A=[yn3_7(x) yn3_8(x)];
                OMb=max(A);
                OMa=min(A);
                c2=(hdR-H7)/(H8-H7);
                c1=1-c2;
                OM=c1*OMa+c2*OMb;
            end  
%
            if(hdR==H7)
                OM=yn3_7(x);
            end
%
            if(H6< hdR && hdR<H7)
                A=[yn3_6(x) yn3_7(x)];
                OMb=max(A);
                OMa=min(A);
                c2=(hdR-H6)/(H7-H6);
                c1=1-c2;
                OM=c1*OMa+c2*OMb;
            end                
%
            if(hdR==H6)
                OM=yn3_6(x);
            end
%
            if(H5< hdR && hdR<H6)
                A=[yn3_5(x) yn3_6(x)];
                OMb=max(A);
                OMa=min(A);
                c2=(hdR-H5)/(H6-H5);
                c1=1-c2;
                OM=c1*OMa+c2*OMb;
            end                        
%
            if(hdR==H5)
                OM=yn3_5(x);
            end
%
            if(H4< hdR && hdR<H5)
                A=[yn3_4(x) yn3_4(x)];
                OMb=max(A);
                OMa=min(A);
                c2=(hdR-H4)/(H5-H4);
                c1=1-c2;
                OM=c1*OMa+c2*OMb;
            end            
%
            if(hdR==H4)
                OM=yn2_4(x);
            end
%
            if(H3< hdR && hdR<H4)
                A=[yn2_4(x) yn2_3(x)];
                OMb=max(A);
                OMa=min(A);
                c2=(hdR-H3)/(H4-H3);
                c1=1-c2;
                OM=c1*OMa+c2*OMb;
            end
%
            if(hdR==H3)
                OM=yn2_4(x);
            end
%
            if(H2< hdR && hdR<H3)
                A=[yn2_3(x) yn2_2(x)];
                OMb=max(A);
                OMa=min(A);
                c2=(hdR-H2)/(H3-H2);
                c1=1-c2;
                OM=c1*OMa+c2*OMb;
            end
%
            if(hdR==H2)
                OM=yn2_4(x);
            end
%
            if(H1< hdR && hdR<H2)
                A=[yn2_2(x) yn2_1(x)];
                OMb=max(A);
                OMa=min(A);
                c2=(hdR-H1)/(H2-H1);
                c1=1-c2;
                OM=c1*OMa+c2*OMb;
            end
%            
            if(hdR<=H1)
                OM=yn2_1(x);
            end            
%
            om=OM*sqrt((E/rho)/a);
            f=om/(2*pi);
            out1=sprintf('   %d  %d %8.4g ',m,n,f);
            disp(out1);
%    
        else
            break;
        end
    end
end
%