disp(' ');
disp(' cantilever_beam.m  ver 1.7  October 2, 2012');
disp(' by Tom Irvine  Email: tomirvine@aol.com ');
disp(' ');
disp(' Normal Modes & Optional Base Excitation ');
disp(' ');
%
close all hidden;
%
clear root;
clear f;
clear H;
clear HA;
clear beta;
clear fn;
clear YY;
%
iu=1;
fig_num=1;
%
disp(' Enter number of modes ');
n = input(' ');
%
[ccc,A,I,L]=geometry_entry(iu);
%
[E,rho,v]=materials(iu);
%
rho=rho*A;   % mass per unit length
%
EI_term = sqrt(E*I/rho);
%
mass=rho*L;
%
sq_mass=sqrt(mass);
%
root=zeros(n,1);
beta=zeros(n,1);
omegan=zeros(n,1);
fn=zeros(n,1);
%
root(1)=1.87510;
root(2)=4.69409;
root(3)=7.85476;
root(4)=10.99554;
%
for i=1:n
    if(i>=5)
        root(i)=(2*i-1)*pi/2;
    end  
    beta(i)=root(i)/L;
    omegan(i)=(beta(i))^2*EI_term;
    fn(i)=omegan(i)/(2*pi);
end
%
C=zeros(n,1);
for i=1:n
    C(i)=-(cos(root(i))+cosh(root(i)))/(sin(root(i))+sinh(root(i)));
end
%
part=zeros(n,1);
for i=1:n
    arg=root(i);
    p2=(sinh(arg)-sin(arg))+C(i)*(cosh(arg)+cos(arg));
    arg=0;
    p1=(sinh(arg)-sin(arg))+C(i)*(cosh(arg)+cos(arg));
    part(i)=(p2-p1)/beta(i);
end
part=part/L;
emm=part.^2;
emm=emm*mass*386;
part=part*sq_mass;
%
ModeShape=@(arg,Co)(cosh(arg)-cos(arg))+Co*(sinh(arg)-sin(arg));
%
ModeShape_dd=...
     @(arg,Co,beta)(beta^2*((cosh(arg)+cos(arg))+Co*(sinh(arg)+sin(arg))));
%     
%
nb=201;
%
for ijk=1:n
%
    clear x;
    clear y;
    x=zeros(nb,1);   
    y=zeros(nb,1);
    dx=L/(nb-1);
    for i=1:nb
        x(i)=(i-1)*dx;
        arg=beta(ijk)*x(i);
        y(i)=ModeShape(arg,C(ijk));
    end
%
    figure(fig_num);
    fig_num=fig_num+1;
    plot(x,y);
    out1=sprintf('Mode %d   fn=%8.4g Hz',ijk,fn(ijk));
    title(out1);
    xlabel('x (inch) ');
    grid on;
%
end
%
disp(' ');
disp(' ');
disp('                             Effective     ');
disp(' mode   fn(Hz)     PF       Modal Mass (lbm)');
for i=1:n
    out1=sprintf(' %d   %7.2f    %8.4f  %8.4f ',i,fn(i),part(i),emm(i));
    disp(out1);
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Perform base excitation analysis? 1=yes 2=no ');
ibase=input(' ');
%
if(ibase==1)
%    
    [damp]=enter_modal_damping(n);
%    
    disp(' ');
    disp(' Enter distance x to point of interest ');
    x=input(' ');
    if(x>L)
        disp(' Warning: x reset to total length');
        x=L;
    end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
 [fig_num,absolute_accel,relative_disp]=...
 continuous_beam_base_frf(ModeShape,sq_mass,beta,C,part,fn,fig_num,x,damp);
%
    disp(' ');
    disp(' Calculate specific case for sine excitation? ');
    disp(' 1=yes  2=no ');
%
    ib=input(' ');
%    
    if(ib==1)
%
       [~]=continuous_beam_base_sine(ModeShape,ModeShape_dd,fn,beta,C,...
                                                part,sq_mass,damp,x,E,I,L);
%
    end
%
end