% DampedBeams.m  Most recent update 2/11/2009
% This script will animate the motion of a viscously damped
% Bernoulli-Euler elastic beam of length L, stiffness EI, mass
% per unit length mu.  It either has an initial deflection curve or
% is driven by a load that is of the form Ao[f(x)g(t)] where f(x) can
% be either a point load at location x = a or a uniform load.
% The temporal function g(t) can be either an impulse, a step or a cosine.
% The constant Ao will have different units depending on the nature
% of f(x) and g(t). The solution is found for four distinct sets of
% boundary conditions. In the case of no forcing function the initial
% deformation may be specified as a static deflection curve or the modal 
% expansion coefficients. The damping is specified as the first modal
% damping ratio, zeta1=c/(2*mu*omega1)where c is the per unit'(
% length damping coefficient
% This m-file was written at the University of Wyoming in the Electrical
% and Computer Engineering Department and is to be distributed without
% cost.
clc
clear all
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
disp('  This script will animate the motion of a viscously damped ')
disp('  Bernoulli-Euler elastic beam of length L, stiffness EI, mass ')
disp('  per unit length mu.  It either has an initial deflection curve or')
disp('  is driven by a load that is of the form Ao[f(x)g(t)] where f(x)')
disp('  can be either a point load at location x = a or a uniform load.')
disp('  The temporal function g(t) can be either an impulse,a step oa a cosine.')
disp('  The constant Ao will have different units depending on the nature')
disp('  of f(x) and g(t).  The solution is found for four distinct sets of')
disp('  boundary conditions. In the case of no forcing function the initial')
disp('  deformation may be specified as a static deflection curve or the')
disp('  modal expansion coefficients. The damping is specified as the first')
disp('  modal damping ratio, zeta1=c/(2*mu*omega1)where c is the per unit')
disp('  length damping coefficient')
disp(' ')
disp('  The four types of boundary conditions are:')
disp('      1-Pinned-Pinned (Simply Supported)')% BC=1 indicates pinned-pinned
disp('      2-Clamped-Free (Cantilever)')% BC=2 indicates clamped-free
disp('      3-Clamped-Pinned')% BC=3 indicates clamped-pinned
disp('      4-Clamped-Clamped')% BC=4 indicares clamped-clamped
disp(' ')
BC=-1;
while((BC~=1)&(BC~=2)&(BC~=3)&(BC~=4))|(k1==1)
    BC=input('  Select an integer 1-4 ');
    k1=isempty(BC);
end
if BC==1
    BCtext='Pinned-Pinned Boundary Conditions';
elseif BC==2
    BCtext='Clamped-Free Boundary Conditions';
elseif BC==3
    BCtext='Clamped-Pinned Boundary Conditions';
else
    BCtext='Clamped-Clamped Boundary Conditions';
end
disp(' ')
disp('  The beam can be unforced or have one of two spatial load distributions:')
disp('      1-No load, initial deflection only')% LT=1 indicates initial conditions only
disp('      2-Point load at x = a')% LT=2 indicates point load in space
disp('      3-Uniform load along the length of the beam')% LT=3 indicates uniform load in space 
LT=-1;
disp(' ')
while((LT~=1)&(LT~=2)&(LT~=3))|(k2==1)
    LT=input('  Select an integer 1-3 ');
    k2=isempty(LT);
end
if LT==2
    disp(' ')
    disp('  The load is located at dimensionless location a/L')
    disp(' ')
    aoverL=-1;
    if BC==2
        while((aoverL<=0)|(aoverL>1)|(k4==1))
            aoverL=input('      a/L = ');
            k4=isempty(aoverL);
        end
    end
    if((BC==1)|(BC==3)|(BC==4))
        while((aoverL<=0)|(aoverL>=1)|(k4==1))
            aoverL=input('      a/L = ');
            k4=isempty(aoverL);
        end
    end 
aoverL1=fix(100*(aoverL))+1;
end

FF=-1;
if(LT~=1)
    disp(' ')
    disp('  The three types of temporal (time) forcing functions are:')
    disp('     1-Impulse function')% FF=1 indicates temporal impulse
    disp('     2-Step function')% FF=2 indicates temporal step
    disp('     3-Cosine function')% FF=3 indicates temporal cosine
    disp(' ')
    while((FF~=1)&(FF~=2)&(FF~=3))|(k3==1)
        FF=input('  Select an integer 1-3 ');
        k3=isempty(FF);
    end
    disp(' ')
    if ((FF==2)|(FF==3))&(LT==2)
        ylabeltext='Displacement, \muL\omega_1^2y(x,t)/A_0';
        disp('  Ao has units of force.')
    elseif ((FF==2)|(FF==3))&(LT==3)
        ylabeltext='Displacement, \mu\omega_1^2y(x,t)/A_0';
        disp('  Ao has units of force/length.')
    elseif (FF==1)&(LT==2)
        ylabeltext='Displacement, \muL\omega_1y(x,t)/A_0';
        disp('  Ao has units of force*time.')
    elseif (FF==1)&(LT==3) 
        ylabeltext='Displacement, \mu\omega_1y(x,t)/A_0';
        disp('  Ao has units of force*time/length.')
    end
end
if FF==3
    xlabeltext='Time, \omegat';
 else
    xlabeltext='Time, \omega_1t';
end
if FF==-1
    ylabeltext='Displacement, y(x,t)/y_m_a_x';
end
if FF==1&LT==2
    fflabel=['A_0\delta(t)\delta(x - ' num2str(aoverL) 'L)'];
elseif FF==1&LT==3
    fflabel='A_0\delta(t)'; 
elseif FF==2&LT==2
    fflabel=['A_0u(t)\delta(x - ' num2str(aoverL) 'L)'];
elseif FF==2&LT==3
    fflabel='A_0u(t)';
elseif FF==3&LT==2
    fflabel=['A_0cos(\omegat)\delta(x - ' num2str(aoverL) 'L)'];
elseif FF==3&LT==3
    fflabel='A_0cos(\omegat)';
end
BetaL=zeros(4,5);
Alpha=zeros(4,5);
BetaL(1,:)=pi*[1 2 3 4 5];
BetaL(2,:)=[1.8751041 4.69409113 7.85475743 10.99554074 14.13716839];
BetaL(3,:)=[3.9266023 7.06858275 10.21017613 13.35176878 16.49336143];
BetaL(4,:)=[4.7300408 7.8532046 10.9956078 14.1371655 17.2787596];
Alpha(2,:)=[.7340955 1.01846644 .99922450 1.000033553 .9999985501];
Alpha(3,:)=[1.0007773 1.00000144 1. 1. 1.];
Alpha(4,:)=[.9825022158 1.000777311 .9999664501 1.00000145 .9999999373];
alpha=Alpha(BC,:);
betaL=BetaL(BC,:);
betaL4=betaL.^4;
fioverf1=(betaL.*betaL)/(betaL(1,1)^2);
xoverL=linspace(0,1,101);
phi=zeros(5,101);
sqrt2=sqrt(2);
f=zeros(1,5);
if BC==1
    for i=1:5
        phi(i,:)=sqrt2*sin(i*pi*xoverL);
    end
else
   for i=1:5
      z=betaL(1,i)*xoverL;
      phi(i,:)=cosh(z)-cos(z)-alpha(1,i)*(sinh(z)-sin(z));
   end
end
% numerically calculate spatial force expansion coefs fi
if LT~=1
    for i=1:5
        if LT==2
            f(1,i)=phi(i,aoverL1);
        elseif LT==3
            f(1,i)=(sum(phi(i,:))-.5*phi(i,101))*.01;
        end
    end
end
zeta1=-1;      
disp(' ')
disp('  The damping is specified in terms of the first mode ')
disp('  damping ratio zeta1.')

while (zeta1<0)|(zeta1>=1)
    disp(' ')
    zeta1=input('  The damping ratio in the first mode is zeta1 = ');
    if isempty(zeta1)==1
        zeta1=0;
    end
end

if FF==3
    disp(' ')
    disp('  The ratio of the forcing frequency to the first natural')
    disp('  frequency is:')
    disp(' ')
    foverf1=input('  f/f1 = ');
end
while FF==3&zeta1==0&((foverf1==fioverf1(1,1))|(foverf1==fioverf1(1,2))|(foverf1==fioverf1(1,3)))
    disp(' ')
    disp('  With a sinusoidal force and zero damping an infinite response ')
    disp('  will rersult, choose another frequency ratio.')
    disp(' ')
    foverf1=input('  f/f1 = ');
    %disp(' ')
end
if (LT==1)
    disp(' ')
    locindex=51;
    text2='midpoint';
    if BC==2
        text2='tip';
        locindex=101;
    end
    disp('  The initial condition may be specified as the static deflection')
    disp('  curve for a uniform load or as a modal expansion by specification')
    disp('  of the modal coefficients for the expansion')
    disp(' ')
    disp(['      1-Static load deflection curve for point load at the beam ' text2])
    disp('      2-Modal expansion coefficients')
    disp(' ')
    IC=-1;
    while((IC~=1)&(IC~=2))|(k7==1)
        IC=input('  Select an integer 1-2 ');
        k7=isempty(IC);
    end
    if IC==1
        yst=zeros(1,101);
        for i=1:5
            yst=yst+(phi(i,locindex)/betaL4(1,i))*phi(i,:);
        end
        norm=max(yst);
        yst=yst/norm;
        for i=1:5
            q0(1,i)=phi(i,locindex)/(betaL4(1,i)*norm);
        end
     elseif IC==2
        b=zeros(1,5);
        disp(' ')
        disp('  The initial condition is specified in terms of a modal');
        disp('  expansion with coefficients qi(0) i=1,...,5');
        disp(' ')
        
        for i=1:5
            bi=input(['    Coefficient q' num2str(i) '(0) = ']);
            if isempty(bi)
                bi=0;
            end
            q0(1,i)=bi;
        end
    end
end
fdioverf1=[];
zeta=[];
for i=1:5
    zeta(1,i)= zeta1/fioverf1(1,i);
    fdioverf1(1,i)=fioverf1(1,i)*sqrt(1-zeta(1,i)^2);
    fdioverfd1(1,i)=fdioverf1(1,i)/sqrt(1-zeta(1,1)^2);
end
Y=zeros(401,101);
tindex=0;
tau1=[];
if(LT==1)%initial displacement only
    for tau=0:.05:20%time loop
        tau1=[tau1 tau];
        tindex=tindex+1;
        for i=1:5%series loop
            Term=q0(1,i)*exp(-zeta(1,i)*fioverf1(1,i)*tau);
            Term=Term*(cos(fdioverf1(1,i)*tau)+(zeta(1,i)/sqrt(1-zeta(1,i)^2))*sin(fdioverf1(1,i)*tau))*phi(i,:);
            Y(tindex,:)=Y(tindex,:)+Term;
        end 
    end
    Ymax=max(max(abs(Y)));
    Y=Y/Ymax;
end
if FF==1%impulse in time forcing function
    for tau=0:.05:20 %time loop
        tau1=[tau1,tau];
        tindex=tindex+1;
        for i=1:5 %series loop
            Term=f(1,i)*exp(-zeta(1,i)*fioverf1(1,i)*tau)*sin(fdioverf1(1,i)*tau)*phi(i,:)/fdioverfd1(1,i);
            Y(tindex,:)=Y(tindex,:)+Term;
        end
    end
end
if FF==2%step in time forcing function
    for tau=0:.05:20 %time loop
        tau1=[tau1,tau];
        tindex=tindex+1;
        for i=1:5 %series loop
            Term1=(1-exp(-zeta(1,i)*fioverf1(1,i)*tau)*cos(fdioverf1(1,i)*tau));
            Term2=Term1-(zeta(1,i)/(sqrt(1-zeta(1,i)^2)))*exp(-zeta(1,i)*fioverf1(1,i)*tau)*sin(fdioverf1(1,i)*tau);
            Term=Term2*f(1,i)*phi(i,:)/(fioverf1(1,i)^2);
            Y(tindex,:)=Y(tindex,:)+Term;
        end
    end
end
if FF==3%cosine in time forcing function, steady state response only
    for tau=0:.05:20 %time loop
        tau1=[tau1,tau];
        tindex=tindex+1;
        sum=zeros(1,101);
        for i=1:5 %series loop
            Gi=f(1,i)/((fioverf1(1,i)^2)-(foverf1^2)+j*2*zeta(1,i)*fioverf1(1,i)*foverf1);
            Term=phi(i,:)*abs(Gi)*cos(tau+angle(Gi));
            sum=sum+Term;
        end
        Y(tindex,:)=sum;
    end
end
ymax=max(max(abs(Y)));
yclampl=[-1 1 .5 .1 -.3 -.5 -1]*ymax*.45;
yclampr=.49*ymax*[ .4 -.4 -.4 .4 .4];
xclampl=[0 0 -.08 -.06 -.09 -.05 0]*1.6;
xclampr=[1 1 1.15 1.15 1];
xsptr=[1 1.15 1.1 1.04 1 1];
ysptr1=.49*ymax*[.53 .53 .83 .78 .9 .53];
ysptr2=-ysptr1;
xbearing=[1.03 1.12 1.03 1.12];
ybearing=.49*ymax*[.47 .47 -.47 -.47];
if FF==2
    low=-.5*ymax;
else
    low=-1.5*ymax;
end
xarrow=[0 0 .03 -.03 0];
yarrow=[.4*low 0 .08*low .08*low 0];
xoverL1=[0 .25 .5 .75 1];

figure(1);clf reset;%plot normal modes
axis([0 1 -2.2 2.2])
hold on
box on
for i=1:5
    plot(xoverL,phi(i,:));
    hold on
    plot([.1 .2],[phi(i,11) -(.4+.3*i)],'k')
    text(.21 ,-(.4+.3*i),['i = ' num2str(i)])
end
text(.05,2.05,['Eigenfunctions for ' BCtext])
xlabel('Distance, x/L')
ylabel('\phi_i(x)')
hold off
text(.05,-2.1,'Press Enter to Continue');
pause 

figure(2);clf reset% do animation 
h1=axes('Ydir','reverse');
axis([-.2 1.2 low 1.5*ymax]);
hold on
box on
xlabel('Distance, x/L');
ylabel(ylabeltext);
if BC==1
    patch(xclampl,yclampl,'r');
    patch(xclampr,yclampr,'r');
    patch(xsptr,ysptr1,'r');
    patch(xsptr,ysptr2,'r');
    plot(xbearing, ybearing, 'o','Color','k','LineWidth',[1.5])
    plot([0 1],[0 0],'o','Color','r');
end
if BC==2
    patch(xclampl,yclampl,'r');
    plot([-.05 0],[0 0],'Color','k','LineWidth',[4.5]);
end
if BC==3
    patch(xclampl,yclampl,'r');
    patch(xclampr,yclampr,'r');
    patch(xsptr,ysptr1,'r');
    patch(xsptr,ysptr2,'r');
    plot(xbearing, ybearing, 'o','Color','k','LineWidth',[1.5])
    plot([-.05 0],[0 0],'Color','k','LineWidth',[4.5]);
    plot([ 1],[ 0],'o','Color','r');
end
if BC==4
   patch(xclampl,yclampl,'r');
   patch(xclampr,yclampr,'r');
   patch(xsptr,ysptr1,'r');
   patch(xsptr,ysptr2,'r');
   plot(xbearing, ybearing, 'o','Color','k','LineWidth',[1.5])
   plot([-.05 0],[0 0],'Color','k','LineWidth',[4.5]);
   plot([1 1.05],[0 0],'Color','k','LineWidth',[4.5]);
end 
box on
if LT==2
    plot(xarrow+aoverL*ones(1,5),2*yarrow);
    text(aoverL,2.1*yarrow(1,1),fflabel)
elseif(LT==3)
    for i=1:5
        plot(xarrow+ones(1,5)*xoverL1(1,i),yarrow);
        plot([0 1],[.4*low .4*low]);
    end
    text(.2,1.2*yarrow(1,1),fflabel)
end
if FF==3|LT==1
    data1=Y(1,:);
else
    data1=zeros(1,101);
end
L=plot(xoverL,data1,'k','EraseMode','xor','LineWidth',[3.5]);
thandl=text(-.1,1.2*ymax,'Press Enter to Animate');
text(.6,1.2*ymax,['\omega_1=(' num2str(betaL(1,1)^2) '/L^2)sqrt(EI/\mu)'])
text(.6,1.4*ymax,['\zeta_1 = ' num2str(zeta1)]);
if LT==2
    text(-.1, 1.4*ymax,['a/L = ' num2str(aoverL)]);
end
if FF==3
    text(.9,1.4*ymax,['\omega/\omega_1 = ' num2str(foverf1)]);
end
hold on
pause
clf
h1=axes('Ydir','reverse');
axis([-.2 1.2 low 1.5*ymax]);
hold on
ylabel(ylabeltext);
xlabel('Distance, x/L');
if BC==1
    patch(xclampl,yclampl,'r');
    patch(xclampr,yclampr,'r');
    patch(xsptr,ysptr1,'r');
    patch(xsptr,ysptr2,'r');
    plot(xbearing, ybearing, 'o','Color','k','LineWidth',[1.5])
    plot([0 1],[0 0],'o','Color','r');
end
if BC==2
    patch(xclampl,yclampl,'r');
    plot([-.05 0],[0 0],'Color','k','LineWidth',[4.5]);
end
if BC==3
    patch(xclampl,yclampl,'r');
    patch(xclampr,yclampr,'r');
    patch(xsptr,ysptr1,'r');
    patch(xsptr,ysptr2,'r');
    plot(xbearing, ybearing, 'o','Color','k','LineWidth',[1.5])
    plot([-.05 0],[0 0],'Color','k','LineWidth',[4.5]);
    plot([ 1],[ 0],'o','Color','r');
end
if BC==4
   patch(xclampl,yclampl,'r');
   patch(xclampr,yclampr,'r');
   patch(xsptr,ysptr1,'r');
   patch(xsptr,ysptr2,'r');
   plot(xbearing, ybearing, 'o','Color','k','LineWidth',[1.5])
   plot([-.05 0],[0 0],'Color','k','LineWidth',[4.5]);
   plot([1 1.05],[0 0],'Color','k','LineWidth',[4.5]);
end 
box on
text(.6,1.2*ymax,['\omega_1=(' num2str(betaL(1,1)^2) '/L^2)sqrt(EI/\mu)'])
hold on
text(.6,1.4*ymax,['\zeta_1 = ' num2str(zeta1)]);
if LT==2
    text(-.1, 1.4*ymax,['a/L = ' num2str(aoverL)]);
end
if FF==3
    text(.9,1.4*ymax,['\omega/\omega_1 = ' num2str(foverf1)]);
end
L1=plot(xoverL,Y(1,:),'k','EraseMode','xor','LineWidth',[4.5]);
thandl1=text(-.1,1.2*ymax,'                      ');
for i=1:401
    set(L1,'YData',Y(i,:));
    pause(.03)
end
hold off
set(thandl1,'String','Press Enter to Continue');
pause 

figure(3);clf reset% y(x,t) vs t for several x
h1=axes('Ydir','reverse');
hold on
if FF==2
    low=-.2*ymax;
else
    low=-1.5*ymax;
end
axis([0 tau1(1,401) low 1.5*ymax]);
hold on
box on
xlabel(xlabeltext);
ylabel(ylabeltext);
Nplot=3;
if BC==2
    Nplot=4;
end
xoverL1=[.25 0.5 0.75 1];
hist=[26 51 76 101];
for i=1:Nplot
    plot(tau1,Y(:,hist(1,i)));
    hold on
    plot([8 10.9],[Y(161,hist(1,i)),Y(161,hist(1,i))+.2*(ymax+.8*i)],'k');
    text(11,Y(161,hist(1,i))+.2*(ymax+.8*i),['x/L = ' num2str(xoverL1(1,i))]); 
end
if LT==2
    text(2, 1.43*ymax,['\zeta_1 = ' num2str(zeta1) '    a/L = ' num2str(aoverL)]);
else
    text(2,1.43*ymax,['\zeta_1 = ' num2str(zeta1)])
end
hold off
text(2,1.2*ymax,'Press Enter to Continue');
pause 

figure(4);clf reset
h1=axes('Ydir','reverse');
if FF==2
    low=-.2*ymax;
else
    low=-1.5*ymax;
end
axis([0 1 low 1.5*ymax]);
hold on
box on
xlabel(xlabeltext);
ylabel(ylabeltext);
for t=1:24:121
    plot(xoverL, Y(t,:));
    hold on
end
tvals=[0 24 48 72 96 120]/400;
if FF==3
    tstr='\omegat = ';
else
    tstr='\omega_1t = ';
end
for j=1:6
    text(.8,ymax*(1.4-j*.15),[tstr num2str(tvals(1,j))]);
   plot([.44 .79],[Y(1+(j-1)*24,45) ymax*(1.4-j*.15)],'k');
end
xlabel('Distance, x/L');
ylabel(ylabeltext)
if LT==2
    text(.45, 1.35*ymax,['\zeta_1 = ' num2str(zeta1) '   a/L = ' num2str(aoverL)]);
else
    text(.45, 1.35*ymax,['\zeta_1 = ' num2str(zeta1)]) 
end
hold off
text(.05,1.35*ymax,'Press Enter to Continue');
pause   

figure(5);clf reset;%3-D plot of y(x,t) vs x and t.
h4=axes('ZDir','reverse');
axis([0 1 0 20 low 1.5*ymax]);
hold on
tau2=linspace(0,20,201);
[X,T]= meshgrid(xoverL,tau2);
respdata=[];
for i=1:401
    if fix(((i+1)/2))==(i+1)/2;
        respdata((i+1)/2,:)=Y(i,:);
    end
end
box on
grid on
mesh(X,T,respdata);
if FF==2
    ztextloc=1.8*ymax;
else
    ztextloc=2*ymax;
end
view(143,30)
text(.6, 20, ztextloc,'Distance, x/L','Rotation',14)
text(1,2,1.1*ztextloc,xlabeltext,'Rotation',-25)
text(0,2,.9*low,'Press Enter to Continue','Rotation',-25)
zlabel(ylabeltext)
pause