% Cantilever beam animation via the steady state solution to the 
% Bernoulli-Euler beam equation.  Beam is forced at the tip by a force 
% of amplitude P0.Input is the ratio of driving freq. to the first 
% natural freq.
% 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);
string1='The ratio of excitation freq. to first natl freq. = ';
fratio=input(string1);
xoverL=linspace(0,1,201);
betaL=sqrt(fratio)*1.8751041;
betax=betaL*xoverL;
cL=cos(betaL);
sL=sin(betaL);
chL=cosh(betaL);
shL=sinh(betaL);
c1=cos(betax);
s1=sin(betax);
ch1=cosh(betax);
sh1=sinh(betax);
delta=-2*(1+chL*cL);
A=(sL+shL)*(c1-ch1)-(cL+chL)*(s1-sh1);
amplitude=(3/((betaL^3)*delta))*A;
lim=max(abs(amplitude));

figure(1);clf;
plot(xoverL,abs(amplitude))
hold on
plot([0 1],[0 0]);
axis([0 1 -.2*lim 1.2*lim]);
text(.8,1.1*lim,['f/f_1 = ',num2str(fratio)])
xlabel('Dimensionless distance from the fixed end, x/L')
ylabel('Dimensionless Amplitude, y(x)/y_s_t')
text(.1,.6*lim,'y_s_t = P_0L^3/3EI')
text(.1,1.1*lim,'Press Enter to Continue')
%xp1=[0 -.1 -.1 0 0];
%yp1=[-.2*lim -.2*lim .2*lim .2*lim -.2*lim];
amplitude=[0 amplitude];
xoverL=[-.07 xoverL];
pause

figure(2);clf;
set(gca,'Box','on')
axis([-.3 1.3 -2*lim 2*lim]);
%P2=plot(xp1,yp1)
xlabel('Distance, x/L')
ylabel('Displacement, y(x,t)/y_s_t')
hold on
texthandle2=text(-.1,-1.6*lim,['f/f_1 = ',num2str(fratio)]);
text(.6,-1.6*lim,'y_s_t = P_0L^3/3EI')
ys=amplitude*0;
plot(xoverL,ys,'k','LineWidth',[1.5]);
hold on
ydata=[ys];
texthandle2=text(0,1.5*lim,'Press Enter to Animate one Frame at a Time');
xlabel('Distance , x/L');
pause
t=linspace(0,8*pi,257);
sine=sin(t);
for k=1:2:65
    ys=amplitude*sine(1,k);
    plot(xoverL, ys,'k','LineWidth',[2.5]);
    hold on
    set(texthandle2,'String', 'Press Enter to Animate One Frame at a Time');
    pause
end
set(texthandle2,'String','Press Enter to Continue');
pause

figure(3):clf;
xpatch=[0 0 -.15 -.08 -.12 -.07 -.1 0];
ypatch=3*lim*[-.2 .2 .14 .05 -.01 -.14 -.18 -.2];
axis([-.3 1.3 -2*lim 2*lim]);
hold on
set(gca,'Box','on');
hold on
P=patch(xpatch,ypatch,'r');
hold on
xlabel('Distance, x/L');
ylabel('Displacement, y(x,t)/y_s_t');
hold on
x1=[-.3 1.3];y1=[0 0];
plot(x1,y1);
hold on
texthandle2=text(-.1,-1.6*lim,['f/f_1 = ',num2str(fratio)]);
text(.6,-1.6*lim,'y_s_t = P_0L^3/3EI')
ys=amplitude*0;
L=plot(xoverL,ys,'k','EraseMode','xor','LineWidth',[2.5]);
ydata=[ys(1,2:202)];
texthandl=text(0,1.5*lim,'Press Enter to Animate');
xlabel('Distance , x/L');
ylabel('Displacement, y(x,t)/y_s_t');
pause
set(texthandl,'String','   ');
t=linspace(0,8*pi,257);
sine=sin(t);
for k=2:257
    ys=amplitude*sine(1,k);;
   %yp2=yp1+ones(1,5)*sine(1,k);
    set(L,'Ydata',ys);
    %set(P,'Ydata',yp2);
    pause(.05)
    ydata=[ydata ;ys(1,2:202)];
end
set(texthandl,'String','Press Enter to Continue');
pause

figure(4);clf;
[X,T]=meshgrid(xoverL(1,2:202),t);
mesh(X,T,ydata)
axis([0 1 0 30 -1.2*lim 1.2*lim]);
text(1,5,-lim,'Dimensionless Time, \omegat','Rotation',-10)
text(1,30,-lim,'Dimensionless Distance, x/L','Rotation',32)
zlabel('Displacement, y(x,t)/y_s_t')
text(0,4, lim,'Press Enter to Continue','Rotation',-10)
view(120,30)