%   beam_long_initial_disp.m   ver 1.3  November 18, 2008
%
%   by Tom Irvine   Email: tomirvine@aol.com
%
disp(' ');
disp(' This program finds the response of a free-free beam to an ');
disp(' initial displacement applied at one end. ');
disp(' The solution is via finite element analysis, modal transient.');
disp(' ');
%
clear k;
clear m;
clear Eigenvalues;
clear ModeShapes;
clear omega;
clear fn;
clear MST;
clear ModalMass;
clear part;
clear QTMQ;
clear r;
%
clear all;
%
tpi=2.*pi;
%
n=2;
%
% disp(' Assume symmetric mass and stiffness matrices. ');
%
disp(' ');
length=input(' Enter the length (inch) ');
length=length*0.0254;
%
disp(' ');
num=input(' Enter the number of elements ');
%
disp(' ');
damp=input(' Enter the damping (assume same for each mode) ');
%
disp(' ');
diameter=input(' Enter the diameter (inch) ');
diameter=diameter/39.37;
%
disp(' ');
dur=input(' Enter duration (sec) ');
%
disp(' ');
sr=input(' Enter sample rate ');
dt=1/sr;
%
vi=zeros(num,1);
di=zeros(num,1);
disp(' ');
din=input(' Enter the initial displacement (in) ');
%
vin=0.;
di(1)=din/39.37;
vi(1)=vin/39.37;
%
h=length/num;
%
disp(' ')
disp(' Enter material:  1=aluminum  2=steel ')
imat=input(' ');
%
if(imat ~=1 & imat ~=2 )
	imat=1;
    disp(' Default to aluminum ');
end
if(imat==1)
    E=6.891e+010;  
    mpv= 2768.; 
end
if(imat==2)
    E=2.067e+011;
    mpv=7887;
end
disp(' ');
%
m=zeros(num,num);
%
k=zeros(num,num);
%
k(1,1)=1.;
m(1,1)=2.;
for(i=2:(num-1))
    k(i,i)=2.;
    m(i,i)=4.;
end
k(num,num)=1.;
m(num,num)=2.;
%
for(i=1:(num-1))
    k(i,i+1)=-1.;
    k(i+1,i)=-1.;
%    
    m(i,i+1)=1.;
    m(i+1,i)=1.;    
end
%
A=diameter^2/4.;
%
k=k*E*A/h;
m=m*h*mpv*A/6.;
%
% disp(' ');
% disp(' The mass matrix is');
% m
% disp(' ');
% disp(' The stiffness matrix is');
% k
%
%  Calculate eigenvalues and eigenvectors
%
[ModeShapes,Eigenvalues]=eig(k,m);
%
omega = sqrt(abs(real(Eigenvalues)));
fn=real(omega/tpi);
%
disp(' ');
disp('  Natural Frequencies = ');
%
for(i=1:num)
    ff(i)=fn(i,i);
end
ff'
%
rdt=1/(10*ff(num));
rsr=1/max(rdt);
if(rsr> sr)
    out1=sprintf('\n Warning: recommended minimum sample rate = %10.5g ',rsr);
    out2=sprintf('\n Change rate?  1=yes 2=no ');
    disp(out1);
    disp(out2);
    icr=input('');
    if(icr==1)
        sr=input(' Enter new rate ');
        dt=1/sr;
    end
end
%
ModeShapes=real(ModeShapes);
MST=ModeShapes';
temp=zeros(num,num);
temp=m*ModeShapes;
QTMQ=MST*temp;
%
clear temp;
%   
for(i=1:num)
    scale(i)=1./sqrt(QTMQ(i,i));
    ModeShapes(:,i) = ModeShapes(:,i)*scale(i);   
    om(i)=omega(i,i);
    om2(i)=om(i)^2;
end 
%
om(1)=0.;
ff(1)=0.;
om2(1)=0.;
%
MST=ModeShapes';
%
% disp(' ');
% disp('  Modes Shapes (column format) =');
% disp(' ');
% ModeShapes
%
%%%%%%%%%%%%%%%%%%%%%
%
nn=round(dur/dt);
%
ndi=MST*m*di;
nvi=MST*m*vi;
%
omegad=om*sqrt(1-damp^2);
%
domega=damp*om;
%
nd(1)=ndi(1);
nv(1)=nvi(1);
na(1)=0.;
%
progressbar % Create figure and set starting time 
%
for(i=1:nn)
    t=dt*(i-1);
    tt(i)=t;
end
%
for(i=1:nn)
%
    progressbar(i/nn) % Update figure  
%
    ac1=2*domega;
%    
    nd=zeros(num,1);
    nv=zeros(num,1);
    na=zeros(num,1);
%    
    for(j=1:num)
%        
         ddd=domega(j)*tt(i);
        arg2=omegad(j)*tt(i);
%    
        expt=exp(-ddd);
%   
        if(omegad(j)>1.0e-04)
                W=(domega(j)/omegad(j));
        else
                W=damp;
        end    
%            
        nd(j) =ndi(j)*expt*(cos(arg2)+W*sin(arg2));
%             
        nv(j) =-domega(j)*nd(j) +       omegad(j)*ndi(j)*expt*(-sin(arg2) + W*cos(arg2));
%            
        na(j) =-domega(j)*nv(j) + ((omegad(j))^2)*ndi(j)*expt*(-cos(arg2) - W*sin(arg2));    
    end    
%
    xd=ModeShapes*nd;
    xv=ModeShapes*nv;
    xa=ModeShapes*na;
%
    for(j=1:num)
        xdm(i,j)=xd(j)*39.37;       
        xvm(i,j)=xv(j)*39.37;        
        xam(i,j)=xa(j)*39.37/386;
    end
%
end
%
for(i=1:num)    
    xm(i)=max(abs(xam(:,i)));
end
% xm'
for(j=1:num)
    p(j)=(j-1)*h/0.0254;
end
sol=[p',xm'];
xm=xm/xm(4);
sol2=[p',xm']
plot(p,xm);
%