%
%    continuous_beam_base_sine.m  ver 1.2   October 15, 2012
%
function[M]=...
continuous_beam_base_sine(ModeShape,ModeShape_dd,fn,beta,C,part,...
                                                   sq_mass,damp,x,E,I,L,iu)
%
n=length(fn);
%
disp(' ');
disp(' Enter base frequency option ');
disp(' 1=fundamental frequency  2=other frequency (Hz) ');
ibf=input(' ');
f_base=fn(1);
if(ibf==2)
   disp(' Enter base excitation frequency (Hz) ');
   f_base=input(' ');
end
%
disp(' ');
disp(' Enter amplitude (G) ');
amp=input(' '); 
%
arg=beta*x;
%
sH=0;
sH_relA=0;
om=2*pi*f_base;
%
for i=1:n
%
   pY=part(i)*ModeShape(arg(i),C(i),sq_mass);
   omn=2*pi*fn(i);
   num=-pY;
   den=(omn^2-om^2)+(1i)*2*damp(i)*om*omn;
   sH=sH+num/den;
%
   num=om^2*pY;
   sH_relA=sH_relA+(num/den); 
%
end
%
sH=abs(sH);
sHA=abs(sH_relA+1);        
%
disp(' '); 
out1 =sprintf('  Base excitation = %8.4g G at %8.4g Hz',amp,f_base);
%
if(iu==1)
    out2 =sprintf('     Response at x=%8.4g inch \n',x);
else
    out2 =sprintf('     Response at x=%8.4g m \n',x);    
end
%
out3 =sprintf(' Transfer Magnitude = %8.4g G/G',sHA);
out4 =sprintf('              Accel = %8.4g G',sHA*amp);
%
if(iu==1)
    out5 =sprintf('           Rel Disp = %8.4g in \n',sH*amp*386);
else
    out5 =sprintf('           Rel Disp = %8.4g m \n',sH*amp*9.81);    
end
%        
disp(out1);  
disp(out2);      
disp(out3);
disp(out4);
disp(out5);
% 
xx=[0 0.5*L L];
%
for j=1:3
%
    M(j)=0;
%
    for i=1:n
       omn=2*pi*fn(i);
       den=(omn^2-om^2)+(1i)*2*damp(i)*om*omn;
%
       arg=beta(i)*xx(j);
       num=-part(i)*ModeShape_dd(arg,C(i),beta(i),sq_mass);
       M(j)=M(j)+(num/den);        
    end
%
end
%
M=abs(M);
M=M*E*I*amp;
%
if(iu==1)
    M=M*386;
else
    M=M*9.81;
end
%
disp('  x      Bending Moment ');
%
if(iu==1)
    disp(' (in)    (in-lbf)  ');
else
    disp(' (in)    (N-m)  ');    
end
%
for j=1:3
    out1=sprintf(' %5.2f   %8.1f',xx(j),M(j));
    disp(out1);
end