%
disp(' ');
disp('  Longitudinal vibration of a ');
disp('  Fixed-Free beam subjected to initial velocity ');
disp(' ');
disp(' By Tom Irvine,  November 10, 2008 ');
disp(' ');
%
clear omega;
clear fn;
%
G =386.;
tpi=2.*pi;
%
c=200000.;  % speed of sound in aluminum or steel 
%
L=input(' Enter beam length (inch) ');
disp(' ');
%
disp(' The natural frequencies (Hz) are ');
%
for(k=1:5)
    j=(2*k-1); 
    omega(j)=j*pi*c/(2*L);
    fn(j)=omega(j)/tpi;
    out1=sprintf(' %9.5g Hz',fn(j));
    disp(out1)
end
%
disp(' ');
h=input(' Enter drop height (inches) ');
v=sqrt(2*G*h);
%
out5 = sprintf('\n The initial velocity is %8.4g in/sec',v);
disp(out5);  
%
damp = 0.05;
%
sr=100000.;
dt=1/sr;
dur=0.030;
%
n=round(dur/dt);
%
B=(-8*v*L)/(pi^2*c);
%
clear acc;
clear tt;
%
tt=linspace(0,n*dt,(n-1)); 
tt=tt-0.004;
m=size(tt);
%
for(i=1:m(2))
    t=tt(i);
    acc(i)=0.;
%    
    if(t>=0)
        for(k=1:5)
             j=(2*k-1);
            acc(i)=acc(i)+((omega(j)^2)/(j^2))*exp(-damp*omega(j)*t)*sin(j*pi/2)*sin(omega(j)*t);   
        end
    end
    acc(i)=B*acc(i)/G;
end
%
plot(tt,acc);
ylabel(' Accel (G)');
xlabel(' Time (sec)');
title('ACCELERATION TIME HISTORY  BEAM  x=L')
%
out5=sprintf('\n The peak acceleration is %10.5g G',max(abs(acc)));
disp(out5);
%
disp(' ');
disp(' The output name is: time_history ');
%
time_history=[tt' acc'];