%
disp(' ');
disp(' beam_shear_inertia.m  ver 1.3  August 7, 2009 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' This program calculates the natural frequency of beam ');
disp(' undergoing bending and shear, with rotory inertia accounted for. ');
disp(' ');
%
disp(' Enter length (inch) ');
L=input(' ');
%
disp(' Enter cross-section geometry:  1=rectangular  2=circular  3=pipe')
icross=input('');
%
if(icross==1)
%
    disp(' Enter width (inch) ');
    b=input(' ');
%
    disp(' Enter thickness (inch) ');
    h=input(' ');   
%   
    k=2/3;
    A=b*h;
    I=(1/12)*b*h^3;
end
if(icross==2)
%
    disp(' Enter diameter (inch) ');
    d=input(' ');
%
    k=3/4;
    A=pi*d^2/4;
    I=(1/64)*pi*d^4;
end
if(icross==3)
%
    disp(' Enter outer diameter (inch) ');
    d2=input(' ');
%
    disp(' Enter inner diameter (inch) ');
    d1=input(' ');
%
    k=3/4;
    A=pi*(d2^2-d1^2)/4;
    I=(1/64)*pi*(d2^4-d1^4);
end
Area=A;
%
disp(' ');
disp(' Select boundary conditions: ');
disp('  1=simply supported at each end ');
disp('  2=fixed-free');
ibc=input(' ');
%
disp(' ')
disp(' Enter material:  1=aluminum  2=steel  3=other')
imat=input('');
%
iflag=0;
if(imat==1)      % aluminum
    E=1.0e+07;
    rho=0.1;
    mu=0.3;
    iflag=1;
end
%
if(imat==2)      % steel
    E=3.0e+07;
    rho=0.285;
    mu=0.3;
    iflag=1;
end
%
if(iflag==0)      % other
    disp(' ');
    disp(' Enter elastic modulus (lbf/in^2) ');
    E=input(' ');
    disp(' Enter mass density (lbm/in^3) ');
    rho=input(' '); 
    disp(' Enter Poisson ratio ');
    mu=input(' ');  
end
%
G=E/(2*(1+mu));
%
disp(' ');
out1=sprintf('     Elastic modulus = %9.4g lbf/in^2 ',E);
out2=sprintf('       Shear modulus = %9.4g lbf/in^2 ',G);
out3=sprintf('       Poisson ratio = %9.4g  ',mu);
out4=sprintf('        Mass density = %9.4g lbm/in^3',rho);
out5=sprintf('                     = %9.4g lbf sec^2/in^4 ',rho/386);
out6=sprintf('  Cross-section Area = %9.4g in^2',A);
out7=sprintf(' Area Moment Inertia = %9.4g in^4',I);
out8=sprintf('        Shear factor = %9.4g in^4',k);
%
disp(out1);
disp(out2);
disp(out3);
disp(out4);
disp(out5);
disp(' ');
disp(out6);
disp(out7);
disp(out8);
%
rho=rho/386.;
%
a2=E*I/(rho*A);
r2=I/A;
r4=r2^2;
r=sqrt(r2);
%
out9=sprintf('  Beam radius of gyration = %9.4g in',r);
disp(out9);
%
clear fn;
%
EKG=E/(k*G);
EKG_term=1+EKG;
%
if(ibc==1)
    for(n=1:3)
        npL=n*pi/L;
        npL2=npL^2;
        npL4=npL^4;
%
        A=(1+r2*(EKG_term)*npL2)^2;
        B=4*r2*(rho*a2/(k*G))*npL4;
%
        if(A>=B)
        else
            disp(' Warning: A<B ');
        end
%
        R=sqrt(A-B);
%
        om2=(1+r2*(1 + EKG)*npL2) - R;
        om2=k*G*om2/(2*r2*rho);
%
        omega=sqrt(om2);
        fn(n)=omega/(2*pi);
    end
%
    disp(' ');
    disp(' Natural Frequencies ');
    disp(' ');
%
    for(n=1:3)
        out1=sprintf(' f%d = %9.5g Hz ',n,fn(n));
        disp(out1);
    end
%
end
%
if(ibc==2)
%
    disp(' ');
    disp(' Southwell-Dunkerley Approximation (bending & shear only)');
    disp(' ');
%
    fb=sqrt(a2);
    fb=3.5156*fb/(2*pi*L^2);
    fs=sqrt(k*G/rho)/(4*L);
%
    A=(1/fb^2)+(1/fs^2);
    fn=sqrt(1/A);
%    
    out1=sprintf(' fb= %8.4g Hz  bending \n fs= %8.4g Hz  shear \n fn= %8.4g Hz  combined',fb,fs,fn);
    disp(out1);
    disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    fstart=fn*0.6;
    df=0.5;
%
    LLL=round((fn-fstart)/0.5);
%
    disp(' Bruch-Mitchell Method (bending, shear & rotary inertia)');
    disp(' ');
%
    A=Area;
    r2=I/(A*L^2);
    s2=E*I/(k*A*G*L^2);
    b2_initial=rho*A*L^4/(E*I);
%
    M=0;     % no end mass
    khat=0;
%
    term_initial=0.5*M*khat^2/(rho*A*L^3);
    MT_initial=M/(rho*A*L^2);
%
    Mbar=0;
    sigma2=0;
%
    clear ddd;
    clear ap;
    clear ff;
%
    progressbar
    LIMIT=50000;
    for(i=1:LIMIT)
        progressbar(i/LLL)
        fn=(i-1)*df+fstart;
        omega=2*pi*fn;
        om2=omega^2;
        b2=b2_initial*om2;
        b=sqrt(b2);
        bsr2=b2*s2*r2;
        term=term_initial*b2;
%
        clear alpha2;
        clear beta2;
%
        sqterm = sqrt((r2-s2)^2+(4/b^2));
%
        alpha2= -(r2+s2) + sqterm;
        beta2=   (r2+s2) + sqterm;
%
        alpha2= alpha2*0.5;
         beta2=  beta2*0.5;       
%
        ap(i)=alpha2;
        ff(i)=fn;
%
        if( alpha2 > 1.0e-20) 
%
        alpha=sqrt(alpha2);
        beta=sqrt(beta2);
%
        bs2=b2*s2;
        MT=MT_initial*bs2;
%
        balpha=b*alpha;
        bbeta =b*beta;
%
        H1= (1-bs2*(alpha2+r2))*L/balpha;
        H2= (1-bs2*(alpha2+r2))*L/balpha;
        H3=-(1+bs2*(beta2-r2))*L/bbeta;
        H4= (1+bs2*(beta2-r2))*L/bbeta;
%
        ba=b*alpha;
        bb=b*beta;
%     
        coshba =cosh(ba);  
        sinhba =sinh(ba);
        cosbb  =cos(bb);  
        sinbb  =sin(bb);   
%
        MC=Mbar*b2;
        MCC=0.5*Mbar*sigma2*b2;
%        
        alphaC=(alpha2+s2)/alpha;
        betaC =(beta2-s2)/beta;  
%
        R1= (b/alpha)*sinhba + MC*coshba;
        R2= (b/alpha)*coshba + MC*sinhba;
        R3=  (b/beta)*sinbb  + MC*cosbb;
        R4= -(b/beta)*cosbb  + MC*sinbb;
%
        R1P= alphaC*( ba*coshba - MCC*sinhba);
        R2P= alphaC*( ba*sinhba - MCC*coshba);
        R3P= -betaC*( bb*cosbb  - MCC*sinbb);
        R4P=  betaC*(-bb*sinbb  - MCC*cosbb);
%
        M=zeros(4,4);
        M(1,1)=1;
        M(1,3)=1;
        M(2,2)=alphaC;
        M(2,4)=betaC;        
%
        M(3,1)= R1;
        M(3,2)= R2;
        M(3,3)= R3;
        M(3,4)= R4;
%        
        M(4,1)= R1P;
        M(4,2)= R2P;
        M(4,3)= R3P;
        M(4,4)= R4P;
%
        num=det(M);
        ddd(i)=abs(num);
%
        if(i>=2 &&  abs(real(num)) > abs(imag(num)))
            if(num*num_before<=0)
                out1=sprintf(' fn=%8.4g Hz',fn);
                disp(out1);
                break;
            end
        end
%
        num_before=num;
%
    end
    end
    progressbar(1);
    disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
end