%
%  continuous_beam_bc.m  ver 1.2  January 29, 2013
%
function[ModeShape,ModeShape_dd,beta,C,fn,part,emm,ibc]=...
                                    continuous_beam_bc(L,n,EI_term,mass,iu)
%
disp(' ');
disp(' Select boundary conditions ');
disp('  1=fixed-free');
disp('  2=fixed-fixed');
disp('  3=simply supported ');
ibc=input(' ');
%
sq_mass=sqrt(mass);
%
root=zeros(n,1);
beta=zeros(n,1);
omegan=zeros(n,1);
fn=zeros(n,1);
C=zeros(n,1);
part=zeros(n,1);
%
if(ibc==1)  % fixed-free
%
    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
%
    for i=1:n
        C(i)=-(cos(root(i))+cosh(root(i)))/(sin(root(i))+sinh(root(i)));
    end
%
    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*sqrt(mass/L^2);
%
    ModeShape=@(arg,Co,sq_mass)...
         (((cosh(arg)-cos(arg))+Co*(sinh(arg)-sin(arg)))/sq_mass);
%
    ModeShape_dd=@(arg,Co,beta,sq_mass)...
         ((beta^2*((cosh(arg)+cos(arg))+Co*(sinh(arg)+sin(arg))))/sq_mass);
%
end
%
%
if(ibc==2)  % fixed-fixed
%
    root(1)=4.73004;
    root(2)=7.8532;
    root(3)=10.9956;
    root(4)=14.13717;
    root(5)=17.27876;
%
    for i=1:n
        if(i>=6)
            root(i)=(0.5 + i)*pi;
        end  
        beta(i)=root(i)/L;
        omegan(i)=(beta(i))^2*EI_term;
        fn(i)=omegan(i)/(2*pi);
    end
%
    for i=1:n
        bL=root(i);
        C(i)=(sinh(bL)+sin(bL))/(cosh(bL)-cos(bL));
    end
%
    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*sqrt(mass/L^2);
%
    ModeShape=@(arg,Co,sq_mass)...
       (((cosh(arg)-cos(arg))-Co*(sinh(arg)-sin(arg)))/sq_mass);
    ModeShape_dd=@(arg,Co,beta,sq_mass)...
       (((cosh(arg)+cos(arg))-Co*(sinh(arg)+sin(arg)))*(beta^2/sq_mass));
%
end
%
if(ibc==3)  % simply-supported
%
    for i=1:n
    %
        root(i)=i*pi;
    %
        beta(i)=root(i)/L;
        omegan=(beta(i))^2*EI_term;
        fn(i)=omegan/(2*pi);
    end
%
    for i=1:n
        part(i)=(-1/(i*pi))*sqrt(2*mass)*(cos(i*pi)-1);
        C(i)=1;
    end
%
    ModeShape=@(arg,Co,sq_mass)(sqrt(2)*sin(arg)/sq_mass);
    ModeShape_dd=@(arg,Co,beta,sq_mass)(beta^2*(-sqrt(2)*sin(arg))/sq_mass);
%
end
%
emm=part.^2;
if(iu==1)
    emm=emm*386;
end    
%
disp(' ');
disp(' ');
disp('                             Effective     ');
%
if(iu==1)
    disp(' mode   fn(Hz)     PF       Modal Mass (lbm)');
else
    disp(' mode   fn(Hz)     PF       Modal Mass (kg)');    
end
%
for i=1:n
    out1=sprintf(' %d   %7.2f    %8.4f  %8.4f ',i,fn(i),part(i),emm(i));
    disp(out1);
end
%