%
disp(' ');
disp(' cant_beam_sine_force.m   ver 1.3   April 3, 2013 ');
disp(' ');
disp(' Transverse vibration of a fixed-free beam, ');
disp(' subjected to a sine force at its free end.     ');
disp(' ');
%
close all;
%
clear beta;
clear bL;
clear omegan;
clear D;
clear H;
clear HM;
clear freq;
clear damp;
clear cc;
clear df;
clear fn;
%
iu=1;
%
[E,I,rho,c,L,A]=geometry_materials();
%
damp=input(' Enter the uniform modal damping ');
%
QM=1/(rho*A*L);
SQM=sqrt(QM);
%
max_nmodes=20;
%
beta=zeros(max_nmodes,1);
%
beta(1)=1.87510;
beta(2)=4.69409;
%
for i=3:max_nmodes
    beta(i)=(2*i-1)*pi/2;
end
%
disp(' ');
nmodes=input(' Enter number of modes to include in analysis: ');
if(nmodes>max_nmodes)
    nmodes=max_nmodes;
end
%
  beta=beta/L;
    bL=beta*L;
%
omegan=beta.*beta;
%
omegan=omegan*sqrt(E*I/(rho*A));
%
clear a;
D=zeros(nmodes,1);
%
for i=1:nmodes
    a=bL(i);
    D(i)=( cos(a)+cosh(a))/(sin(a)+sinh(a) );
end
%
disp(' ');
disp(' Natural Frequency (Hz) ');
fn=zeros(nmodes,1);
for k=1:nmodes
    fn(k)=omegan(k)/(2*pi);
    out1=sprintf(' %8.4g ',fn(k));
    disp(out1);
end
%
%%%%%%%%%%%%%%%%%%%%%%%%
  Y = inline('SQM*((cosh(arg)-cos(arg))-Di*(sinh(arg)-sin(arg)))','Di','arg','SQM')
DDY = inline('beta^2*SQM*((cosh(arg)+cos(arg))-Di*(sinh(arg)+sin(arg)))','Di','arg','SQM','beta')
%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
disp(' ');
B=input(' Enter force amplitude (lbf) ');
%
disp(' ');
f=input(' Enter force frequency (Hz) ');
%
T=0.5/f;
%
alpha=pi/T;
alpha2=alpha^2;
%
%%%
%
srs=max([ 20*max(fn) 24/T ]);
%
disp(' ');
out1=sprintf(' Recommended sample rate > %8.4g samples/sec',srs);
disp(out1);
disp(' ');
sr=input(' Enter the sample rate (samples/sec) ');
%
dt=1/sr;
%
durr=20*T;
%
out1=sprintf('\n Suggested analysis duration = %8.4g sec \n',durr);
disp(out1);
%
dur=input(' Enter the analysis duration (sec) ');
%
%
nt=round(dur/dt);
%
clear n;
clear nd;
clear nv;
clear na;
clear nvelox;
clear naccel;
clear x;
clear v;
clear accel;
clear force;
clear xL;
clear vL;
clear aL;
%
x=zeros(nt,2);
v=zeros(nt,2);
accel=zeros(nt,2);
%
omegan2=zeros(nmodes,1);
domegan=zeros(nmodes,1);
omegad=zeros(nmodes,1);
omegad2=zeros(nmodes,1);
%
AA=zeros(nmodes,1);
AAT=zeros(nmodes,1);
%
U1=zeros(nmodes,1);
U2=zeros(nmodes,1);
%
V1=zeros(nmodes,1);
V2=zeros(nmodes,1);
%
den=zeros(nmodes,1);
C=zeros(nmodes,1);
YL=zeros(nmodes,1);
%
Fdd=zeros(nmodes,1);
AO=zeros(nmodes,1);
%
for j=1:nmodes
%
     omegad(j)=omegan(j)*sqrt(1-(damp)^2);
    omegad2(j)=omegad(j)^2;
    domegan(j)=damp*omegan(j);
    omegan2(j)=omegan(j)^2;
%
    den(j)=(alpha2-omegan2(j))^2+(2*domegan(j)*alpha)^2;
%
    arg=beta(j)*L;
    YL(j)=Y(D(j),arg,SQM);   
%
    Fdd(j)=B*YL(j)/den(j);
%
    U1(j)=-2*domegan(j)*alpha;
    V1(j)=-(alpha2-omegan2(j));
    U2(j)=2*domegan(j)*omegad(j);
    V2(j)=alpha2+omegan2(j)*(-1+2*damp^2);
%
    AO(j)=alpha/omegad(j);
%
end
%
n=zeros(nmodes,nt);
nvelox=zeros(nmodes,1);
naccel=zeros(nmodes,1);
%
xL=zeros(nt,1);
vL=zeros(nt,1);
aL=zeros(nt,1); 
%
ndT=zeros(nmodes,1);
nvT=zeros(nmodes,1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
t=zeros(nt,1);
force=zeros(nt,1);
%
progressbar;
for i=1:nt
   progressbar(i/nt); 
%  
   t(i)=dt*(i-1);
%   
   at=alpha*t(i);  
%
     sin_at=sin(at);
     cos_at=cos(at);
%
     force(i)=B*sin_at;     
%
     for j=1:nmodes
%
         ee=exp(-domegan(j)*t(i));
        dee=-domegan(j)*ee;
       ddee=(domegan(j))^2*ee;
%
        omdt=omegad(j)*t(i);
        sin_odt=sin(omdt);
        cos_odt=cos(omdt);
%
        n(j,i)=         (U1(j)*cos_at +V1(j)*sin_at)+...
               AO(j)*ee*(U2(j)*cos_odt+V2(j)*sin_odt);
%
       nvelox(j)=   alpha*(-U1(j)*sin_at +V1(j)*cos_at)+...
                AO(j)*dee*(U2(j)*cos_odt+V2(j)*sin_odt)+...
       omegad(j)*AO(j)*ee*(-U2(j)*sin_odt+V2(j)*cos_odt);                         
%
            na1=               alpha2*(-U1(j)*cos_at  - V1(j)*sin_at);
            na2=        omegad2(j)*ee*(-U2(j)*cos_odt - V2(j)*sin_odt)...
                     +2*omegad(j)*dee*(-U2(j)*sin_odt + V2(j)*cos_odt)... 
                     +           ddee*( U2(j)*cos_odt + V2(j)*sin_odt);             
%  
            na2=na2*AO(j);  
%            
            naccel(j)=na1+na2;   
%
        n(j,i)=n(j,i)*Fdd(j);
        nvelox(j)=nvelox(j)*Fdd(j);
        naccel(j)=naccel(j)*Fdd(j);
     end      
%
%  Solve for response at x=L
%
   for j=1:nmodes
          xL(i)=xL(i) +YL(j)*n(j,i);
          vL(i)=vL(i) +YL(j)*nvelox(j);
          aL(i)=aL(i) +YL(j)*naccel(j);   
   end
%   
end
%
pause(0.2);
progressbar(1);
%
aL=aL/386;
%
tmin=0;
tmax=max(t);
%
while(1)
%
    disp(' ');
%
    figure(1);
    plot(t,xL);
    axis([tmin,tmax,min(ylim),max(ylim)]);
    grid on;
    title('Displacement at x=L');
    ylabel('Disp (in)');
    xlabel('Time (sec)');
%
    figure(2);
    plot(t,vL);
    axis([tmin,tmax,min(ylim),max(ylim)]);    
    grid on;
    title('Velocity at x=L');
    ylabel('Vel (in/sec)');
    xlabel('Time (sec)');
%
    figure(3);
    plot(t,aL);
    axis([tmin,tmax,min(ylim),max(ylim)]);   
    grid on;
    title('Acceleration at x=L');
    ylabel('Accel (G)');
    xlabel('Time (sec)');
%
    disp(' ');
    disp(' Replot with different end time?');
    disp('  1=yes  2=no ');
    ijk=input(' ');
%
    if(ijk==2)
        break;
    end
%
    tmax=input(' Enter plot end time(sec) ');
%
end
%
clear FD;
clear FV;
clear FA;
%
FD=[t force xL];
FV=[t force vL];
FA=[t force aL];
%
mx=max(abs(xL));
mv=max(abs(vL));
ma=max(abs(aL));
%
disp(' ');
disp(' Maximum Response Levels ');
out1=sprintf(' Displacement = %8.4g in',mx);
out2=sprintf('     Velocity = %8.4g in/sec',mv);
out3=sprintf(' Acceleration = %8.4g G',ma);
%
disp(out1);
disp(out2);
disp(out3);
%