% beamvibration.m
% 7/28/2008
% Animate the free vibration of a uniform Bernoulli-Euler cantilever 
% beam with initial conditions composed of quadratic and
% cubic factors in x/L i.e. y(x,0)=.667(x/L)^2+.333*(x/L)^3
% The initial velocity is zero.
% This m-file was written at the University of Wyoming in the Electrical
% and Computer Engineering Department and is to be distributed without
% cost.

clear all
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
betaL=[1.875 4.694 7.8547 10.995 14.1371];
betaLsq=betaL.*betaL;
alpha=[.7340995 1.018664 .9992245 1.00003355 .99999855];
xoverL=linspace(0,1,201);
omega1t=linspace(0,4*pi,401);
sn=[];
y0=.667*(xoverL.^2)+.333*(xoverL.^3);% initial deflection
b=zeros(1,5);
nterm=5;
for i=1:nterm
    z=betaL(1,i)*xoverL;
    phi(i,:)=cosh(z)-cos(z)-alpha(1,i)*(sinh(z)-sin(z));%eigenfunctions
    cs(i,:)=cos(betaLsq(1,i)*omega1t/betaLsq(1,1));
    b(1,i)=0.005*sum(y0.*phi(i,:));%calc. Fourier coef. of initial def.
end
xoverL=[-.06 xoverL];
framedata=[];

% calculate the response data
for k=1:401
    sum=zeros(1,201);
    for i=1:nterm
        sum=sum+b(1,i)*cs(i,k)*phi(i,:);
    end
    framedata=[framedata;sum];
end
lim=max(framedata(1,:));
framedata=[zeros(401,1) framedata];

figure(1);clf;
h1=axes('YDir','reverse');
axis([-.2 1.2 -1.2*lim 1.2*lim])
hold on
xpatch1=[0 0 -.1 -.08 -.09 0];
ypatch1=[.2 -.2 -.15 -.05 .05 .2]*lim/.4;
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Displacement, y(x,t)/y_0')
patch(xpatch1,ypatch1,'r')
hold on
S=plot(xoverL,framedata(1,:),'k','LineWidth',[3.5]);
hold on
box on
text(0,1.05*lim,'y(x,0)=y_0[0.667(x/L)^2+0.333(x/L)^3]');
texthandl=text(.1 ,-1.1*lim,'Press Enter to Animate One Frame at a Time');
pause
for k=11:20:201
    plot(xoverL,framedata(k,:),'k','LineWidth',[3.5]);
    set(texthandl,'String','Press Enter to Continue');
    pause
end
hold off
set(gca,'Box','on')

figure(2);clf;%plot y(x,t) vs t for various x's
%h2=axes('YDir','reverse');
hold on
axis([0 4*pi -1.2*lim 1.2*lim]);
hold on
plot(omega1t,framedata(:,1));
hold on
for j=42:40:202
    plot(omega1t,framedata(:,j));
    hold on
    plot([omega1t(1,201) 8],[framedata(201,j) .15+.2*(j/40)],'k')
    text(8.1,.15+.2*(j/40),['x/L = ' num2str((j-2)/200)]);
    hold on
end
xlabel('Dimensionless Time, \omega_1t')
ylabel('Dimensionless Displacement, y(x,t)/y_0')
box on
text(1,-1.1*lim,'Press Enter to Continue')
hold off
pause
%now do animation

figure(3);clf;
h3=axes('YDir','reverse');
axis([-.2 1.2 -1.2*lim 1.2*lim])
hold on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Displacement, y(x,t)/y_0')
hold on
patch(xpatch1,ypatch1,'r')
L=plot(xoverL, zeros(1,202),'k','EraseMode','xor','LineWidth',[3.5]);
plot([xoverL(1,1) xoverL(1,2)],[0 0],'k','Linewidth',[3.5]);
hold on
box on
texthandl1=text(.1,-1.1*lim,'Press Enter to Set Initial Condition');
text(0,1.05*lim,'y(x,0)=y_0[0.667(x/L)^2+0.333(x/L)^3]');
pause
set(L,'Ydata',framedata(1,:))
set(texthandl1,'String','Press Enter to Animate');
pause
set(texthandl1,'String','                       ');
for i=2:401
    set(L,'Ydata',framedata(i,:))
    pause(.05)
end
hold off
set(texthandl1,'String','Press Enter to Continue');
pause

figure(4);clf;
h4=axes('ZDir','reverse');
hold on
axis([0 1 0 4*pi -1.05*lim 1.05*lim])
hold on
[X,Y]=meshgrid(xoverL(1,2:202),omega1t);
mesh(X,Y,framedata(:,2:202))
xlabel('Distance, x/L','Rotation',-31)
ylabel('Dimensionless Time, \omega_1t','Rotation',11)
zlabel('Dimensionless Displacement, y(x,t)/y_0')
view(60,30)
grid on
text(0,4*pi,-.9,'Press Enter to Continue','Rotation',-31)