% Animate the free vibration of a uniform Bernoulli-Euler Clamped-Clamped 
% beam with initial conditions composed of quadratic and
% cubic factors in x/L i.e. y(x,0)=2(x/L)^2-(x/L)^3-4(x/L)^4+3(x/L)^5
% 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=[4.7300408 7.8532046 10.9956078 14.1371655 17.2787596];
betaLsq=betaL.*betaL;
alpha=[.9825022158 1.000777311 .9999664501 1.000001450 .9999999373];
xoverL=linspace(0,1,101);
omega1t=linspace(0,4*pi,401);
sn=[];
% initial condition function
y0=2*(xoverL.^2)-1*(xoverL.^3)-4*(xoverL.^4)+3*(xoverL.^5);
y0max=max(abs(y0));
nterm=5;
b=zeros(1,5);
for i=1:nterm
    z=betaL(1,i)*xoverL;
    phi(i,:)=cosh(z)-cos(z)-alpha(1,i)*(sinh(z)-sin(z));
    cs(i,:)=cos(betaLsq(1,i)*omega1t/betaLsq(1,1));
    b(1,i)=.01*sum(phi(i,:).*y0);
end
xoverL=[-.06 xoverL 1.06];
framedata=[];
for k=1:401
    sum=zeros(1,101);
    for i=1:nterm
        sum=sum+b(1,i)*cs(i,k)*phi(i,:);
    end
    framedata=[framedata;sum];
end
lim=1;
framedata=[zeros(401,1) framedata zeros(401,1)]/y0max;

figure(1);clf;
h1=axes('YDir','reverse');
axis([-.2 1.2 -1.2*lim 1.2*lim])
xpatch1=[0 0 -.1 -.08 -.09 0];
ypatch1=[.2 -.2 -.15 -.05 .05 .2]*lim/.4;
xpatch2=ones(1,6)-xpatch1;
xlabel('Dimensionless Distance, x/L')
ylabel('Dimensionless Displacement, y(x,t)/y_0_m_a_x')
patch(xpatch1,ypatch1,'r');
patch(xpatch2,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)=2(x/L)^2-(x/L)^3-4(x/L)^4+3(x/L)^5');
texthandl=text(.1 ,-1.1*lim,'Press Enter to Animate One Frame at a Time');
pause
for k=11:10:201
    plot(xoverL,framedata(k,:),'k','LineWidth',[3.5]);
    set(texthandl,'String','Press Enter to Continue');
    pause
end
set(gca,'Box','on')

figure(2);clf;
h2=axes('YDir','reverse');
hold on
axis([0 4*pi -1.2*lim 1.2*lim]);
plot(omega1t,framedata(:,2));
hold on
for j=12:10:102
    plot(omega1t,framedata(:,j));
    hold on
end
xlabel('Dimensionless Time, \omega_1t')
ylabel('Dimensionless Displacement, y(x,t)')
text(1,-1.1*lim,'Press Enter to Continue')
pause

%now do animation
figure(3);clf;
h3=axes('YDir','reverse');
axis([-.2 1.2 -1.2*lim 1.2*lim])
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Displacement, y(x,t).y_0_m_a_x')
hold on
patch(xpatch1,ypatch1,'r')
patch(xpatch2,ypatch1,'r')
L=plot(xoverL, zeros(1,103),'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)=2(x/L)^2-(x/L)^3-4(x/L)^4+3(x/L)^5');
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(.04)
end
hold on
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:102),omega1t);
mesh(X,Y,framedata(:,2:102))
xlabel('Distance, x/L','Rotation',-31)
ylabel('Dimensionless Time, \omega_1t','Rotation',11)
zlabel('Dimensionless Temperature, y(x,t)/y_0_m_a_x')
view(60,30)
grid on
text(0,4*pi,-.9,'Press Enter to End','Rotation',-31)
pause