%   beam_bending_initial_disp.m   ver 1.1  November 20, 2008
%
%   by Tom Irvine   Email: tomirvine@aol.com
%
disp(' ');
disp(' This program finds the response of a free-free beam with a rectangular cross-section');
disp(' to an 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 ');
%
ndof=2*(num+1);
%
disp(' ');
damp=input(' Enter the damping (assume same for each mode) ');
%
disp(' ');
thickness=input(' Enter the thickness (inch) ');
thickness=thickness/39.37;
%
disp(' ');
width=input(' Enter the width (inch) ');
width=width/39.37;
%
inertia=width*(thickness^3)/12;
%
disp(' ');
dur=input(' Enter duration (sec) ');
%
disp(' ');
sr=input(' Enter sample rate ');
dt=1/sr;
%
vi=zeros(ndof,1);
di=zeros(ndof,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(ndof,ndof);
%
k=zeros(ndof,ndof);
%
ml(1,1)=156;
ml(1,2)=22;
ml(1,3)=54;
ml(1,4)=-13;
%
ml(2,2)=4;
ml(2,3)=13;
ml(2,4)=-3;
%
ml(3,3)=156;
ml(3,4)=-22;
%
ml(4,4)=4;
%
kl(1,1)=12;
kl(1,2)=6;
kl(1,3)=-12;
kl(1,4)=6;
%
kl(2,2)=4.;
kl(2,3)=-6;
kl(2,4)=2;
%
kl(3,3)=12;
kl(3,4)=-6;
%
kl(4,4)=4;
%
for(nv=1:num)
%
		i=2*nv-1;
%
		k(i,i)  =k(i,i)  +kl(1,1);
		k(i,i+1)=k(i,i+1)+kl(1,2);
		k(i,i+2)=k(i,i+2)+kl(1,3);
		k(i,i+3)=k(i,i+3)+kl(1,4);
%		
        j=i+1;
%
		k(j,j)  =k(j,j)  +kl(2,2);
		k(j,j+1)=k(j,j+1)+kl(2,3);
		k(j,j+2)=k(j,j+2)+kl(2,4);
%
        j=i+2;
%
		k(j,j)  =k(j,j)  +kl(3,3);
		k(j,j+1)=k(j,j+1)+kl(3,4);
%
        j=i+3;
%
		k(j,j) =k(j,j)+kl(4,4);
%
		m(i,i)  =m(i,i)  +ml(1,1);
		m(i,i+1)=m(i,i+1)+ml(1,2);
		m(i,i+2)=m(i,i+2)+ml(1,3);
		m(i,i+3)=m(i,i+3)+ml(1,4);
%		
        j=i+1;
%
		m(j,j)  =m(j,j)  +ml(2,2);
		m(j,j+1)=m(j,j+1)+ml(2,3);
		m(j,j+2)=m(j,j+2)+ml(2,4);
%
        j=i+2;
%
		m(j,j)  =m(j,j)  +ml(3,3);
		m(j,j+1)=m(j,j+1)+ml(3,4);
%
        j=i+3;
%
		m(j,j)  =m(j,j)+ml(4,4);
%  
end
%
for(i=1:ndof)
    for(j=(i+1):ndof)
        m(j,i)=m(i,j);
        k(j,i)=k(i,j);
    end
end
%
k=k*E*inertia/(h^3);
m=m*h*mpv*thickness*width/420.;
%
% 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:ndof)
    ff(i)=fn(i,i);
end
ff';
fs=sort(ff)
%
disp(' ref 1 ');
%
rdt=1/(10*ff(ndof));
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
%
disp(' ref 2 ');
%
ModeShapes=real(ModeShapes);
MST=ModeShapes';
temp=zeros(ndof,ndof);
temp=m*ModeShapes;
QTMQ=MST*temp;
%
disp(' ref 3 ');
%
clear temp;
%   
for(i=1:ndof)
    scale(i)=1./sqrt(QTMQ(i,i));
    ModeShapes(:,i) = ModeShapes(:,i)*scale(i);   
    om(i)=omega(i,i);
    om2(i)=om(i)^2;
end 
%
disp(' ref 4 ');
%
om(1)=0.;
ff(1)=0.;
om2(1)=0.;
%
om(2)=0.;
ff(2)=0.;
om2(2)=0.;
%
MST=ModeShapes';
%
% disp(' ');
% disp('  Modes Shapes (column format) =');
% disp(' ');
% ModeShapes
%
%%%%%%%%%%%%%%%%%%%%%
%
disp(' ref 5 ');
%
nn=round(dur/dt);
%
ndi=MST*m*di;
nvi=MST*m*vi;
%
omegad=om*sqrt(1-damp^2);
%
domega=damp*om;
%
disp(' ref 6 ');
%
progressbar % Create figure and set starting time 
%
for(i=1:nn)
    t=dt*(i-1);
    tt(i)=t;
end
%
for(j=1:ndof)
        if(omegad(j)>1.0e-04)
                WW(j)=domega(j)/omegad(j);
        else
                WW(j)=damp;
        end      
end
%
    sW=size(WW);
%    
    if(sW(1)<sW(2))
        WW=WW';
    end
    WW=WW';
    omd2=omegad.*omegad; 
%
for(i=1:nn)
%
    progressbar(i/nn) % Update figure  
%
    ac1=2*domega;
%    
    nd=zeros(ndof,1);
    nv=zeros(ndof,1);
    na=zeros(ndof,1);
%
    mmm=ndof;
%
    j=1:mmm;
    expt=exp(-domega(j)*tt(i));  
    cosa=cos(omegad(j)*tt(i));
    sina=cos(omegad(j)*tt(i));
%
    se=size(expt);
    sn=size(ndi); 
    sv=size(nv);  
%
    if(se(1)<se(2))
        expt=expt';
    end
    if(sn(1)<sn(2))
        ndi=ndi';
    end
    if(sv(1)<sv(2))
        nv=ndi';
    end   
%
clear Wsina;
clear Wcosa;
%
    A1=ndi.*expt;   
%
    Wsina=WW.*sina;
    Wcosa=WW.*cosa;    
%
      cws= cosa + Wsina;
     mswc=-sina + Wcosa;
    mcmws=-cosa - Wsina;
%
    A1=A1';  
    nd = A1.*cws; 
    CV=A1.*mswc;
    CA=A1.*mcmws;
    
    nv =-domega.*nd +  omegad.*CV;
%          
    na = domega.*nv +  omd2.*CA;    
%
    xd=ModeShapes*nd';
    xv=ModeShapes*nv';
    xa=ModeShapes*na';
%
    kk=1;
    for(j=1:round(mmm/2))
        jj=2*j-1;
        xdm(i,kk)=xd(jj)*39.37;       
        xvm(i,kk)=xv(jj)*39.37;        
        xam(i,kk)=xa(jj)*39.37/386;
        kk=kk+1;
    end
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear filtered_data;
%
for(nvk=1:round(max(size(MST))/2))
%
disp(' ')
disp(' Butterworth_filter.m ')
disp(' ver 2.1   November 3, 2008 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' Butterworth filter, sixth-order, infinite impulse response,');
disp(' cascade with refiltering option for phase correction ');                 
disp(' ')
disp(' The time history must be in a two-column matrix format: ')
disp(' Time(sec)  & amplitude ')
disp(' ')
%
clear THM;
clear tim;
clear y;
clear yt;
clear a;
clear b;
%
y  =xam(:,nvk);
tim=tt;
n = max(size(y));
ns = n;
%
tmx=max(tim);
tmi=min(tim);
% 
out3 = sprintf('\n start  = %g sec    end = %g sec            ',tmi,tmx);
disp(out3);
dt=(tmx-tmi)/(n-1);
sr=1./dt;
disp(' ')
disp(' Time Step ')
%
out5 = sprintf(' dt     = %8.4g sec  ',dt);
disp(out5)
%
disp(' ')
disp(' Sample Rate ');
out5 = sprintf(' sr     = %8.4g samples/sec  ',1/dt);
disp(out5)
%
ncontinue=1;
%
if(ncontinue==1)
%
fl=0.;
fh=0.;
iband=0;  
%
iband=1;
%
f=0.;
fl=0.;
fh=0.;
%
if(iband == 1)
%          fl=input(' Enter lowpass frequency (Hz) ');
          fl=10000.;
          f=fl;
end
%
disp(' ');
disp('  Refiltering for phase correction? ');
disp(' 1=yes (appropriate for transient signal)');
disp(' 2=no  (appropriate for steady-state signal)');
% iphase=input(' ');
iphase=1;
%
if(iphase ~= 1 & iphase ~=2 )
    iphase=1;
end    
%
freq=f;
%
iflag=1;
%
%****** calculate coefficients *******
%
if(iflag ~= 999 & iband ~=3)
%
    [a,b,iflag] = filter_coefficients(f,dt,iband,iflag);
%    
end
%
if(iflag < 900 )
%
    if(iband == 1 | iband ==2)  % lowpass or highpass
%		 
         if(iphase==1)   % refiltering
			[y]=apply_filter(y,iphase,ns,a,b);
			[y]=apply_filter(y,iphase,ns,a,b); 
	      else
			[y]=apply_filter(y,iphase,ns,a,b);	
          end
%     
    end
%
    if(iband == 3)  % bandpass
 %     
		 f=fh;
		 freq=f;
%
		 disp(' Step 1');
		 iband=2;
         [a,b,iflag] = filter_coefficients(f,dt,iband,iflag);
 %		 
         if(iphase==1) % refiltering
			[y]=apply_filter(y,iphase,ns,a,b);
			[y]=apply_filter(y,iphase,ns,a,b);  
	      else
			[y]=apply_filter(y,iphase,ns,a,b);	
          end
%     
		 f=fl;
		 freq=f;
%
		 disp(' Step 2');
		 iband=1;
         [a,b,iflag] = filter_coefficients(f,dt,iband,iflag);
%		 
         if(iphase==1) % refiltering 
			[y]=apply_filter(y,iphase,ns,a,b);
			[y]=apply_filter(y,iphase,ns,a,b);  
	      else
			[y]=apply_filter(y,iphase,ns,a,b);	
          end
%          
    end
%
disp(' stats ');
%
    mu=mean(y)
    sd=std(y)
    mx(nvk)=max(abs(y));
%
    for(iq=1:max(size(y)))
        filtered_data(nvk,iq)=y(iq);
    end
%
    max(size(y))
    size(filtered_data)
    disp(' End filtered_data copy ');
%   
else
    disp(' ')
    disp('  Abnormal termination.  No output file generated. ');
end
end
end
%
disp(' Begin xm copy ');
%
nn=round(max(size(MST))/2)
%
clear xm;
clear p;
%
for(i=1:nn)    
    xm(i)=max(abs(filtered_data(i,:)));
     p(i)=(i-1)*h/0.0254;
end
%
disp(' End xm copy ');
%
% xm'
%
sol=[p',xm'];
xm=xm/xm(4);
sol2=[p',xm']
figure(1);
plot(p,xm);
figure(2);
plot(tt',filtered_data(round(nn/2),:))
figure(3);
plot(tt',filtered_data(nn,:))