disp(' ');
disp(' clamped_pinned_clamped.m   ver 1.0  January 27, 2012 ');
disp(' ');
disp(' by Tom Irvine  Email: tomirvine@aol.com');
disp(' ');
disp(' This program calculates the natural frequencies ');
disp(' and mode shapes of a clamped-pinned-clamped beam. ');
%
close all
clear C;
clear x;
clear y;
clear XS;
clear YS;
clear length;
clear beta;
clear beta_L;
clear fn;
clear ModeShapes;
%
tpi=2*pi;
%
k=1;
%
disp(' ');
icross=input(' Enter the cross-section: \n 1=rectangle  2=solid cylinder  3=other ');
disp(' ');
if(icross==1)
			width=input(' Enter the width (inch) ');
			thick=input(' Enter the thickness (inch) ');
			A=width*thick;
			I=(1./12.)*width*(thick^3.);
end
if(icross==2)
			diam=input(' Enter the diameter (inch) ');
			A=(pi/4.)*(diam^2.);
			I=(pi/64.)*(diam^4.);
end
if(icross ~=1 & icross ~=2)
            disp(' ');
            I = input(' Enter the area moment of inertia (in^4) ');
            %
            disp(' ');
            A = input(' Enter the cross-section area (in^2) ');
end
%
disp(' ');
imat=input(' Enter material: \n 1=aluminum  2=steel  3=other ');
%
if(imat==1)
    rho=0.1;
	E=1.0e+07;
end
if(imat==2)
    rho=0.29;
    E=3.0e+007;
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) ');           
 end
%
mass_per_length=rho*A/386.;
%
disp(' ');
L = input(' Enter the total length (inch) ');
%
disp(' ');
disp(' Enter the distance of the pinned location ');
disp(' relative to the clamped end. ');
%
a=input(' ');
a_orig=a;
%
b=(L-a);
%
x=0.1+linspace(0,24,1001);
%
a=a/L;
b=1-a;
LL=1;
%
%  beam_type =  1 clamped_pinned_free
%               2 clamped_pinned_clamped
%
beam_type=2;
%
[beta_L]=find_root_mbeam(x,a,b,beam_type);
%
beta_L=beta_L'
%
n=length(beta_L);
%
disp(' ');
disp(' Natural Freq (Hz) ');
%
fig_num=1;
%
clear x;
%
x=linspace(0,L,501);
%
numx=length(x);
%
y=zeros(numx,n);
%
a=a_orig;
b=L-a;
%
for i=1:5
    beta=beta_L(i)/L;
    omega=beta^2*sqrt(E*I/mass_per_length);
    fn(i)=omega/tpi;
%
    beta_a=beta_L(i)*a/L;
    beta_b=beta_L(i)*b/L;
    [q,C,B]=cpc_C_coefficients(beta_a,beta_b);
%
    for j=1:numx
        if(x(j)<=a)
            W=beta_a*x(j)/a;
            y(j,i)=q(1)*sinh(W)+q(2)*cosh(W)+q(3)*sin(W)+q(4)*cos(W);
        else
            W=beta_b*(x(j)-a)/b;
            y(j,i)=q(5)*sinh(W)+q(6)*cosh(W)+q(7)*sin(W)+q(8)*cos(W);            
        end
    end
%
    Z=max(abs(y(:,i)));
    y(:,i)=y(:,i)/Z;
%
    out1=sprintf(' %8.4g ',fn(i));
    disp(out1);
%
    figure(fig_num);
    plot(x,y(:,i));
    out1=sprintf(' Mode %i  Freq=%8.4g Hz',i,fn(i));
    title(out1);
    ylabel('Modal Displacement');
    xlabel('x');
    grid on;
    axis([0,max(XLim),-1.2,1.2]);
    fig_num=fig_num+1;
    ModeShapes(:,i+1)=y(:,i);
end
ModeShapes(:,1)=x(:);