% Solves the forced beam vibration problem for a uniform cantilever beam
% subject to a suddenly applied constant force f0 over the length of the
% span and zero initial conditions in displacement and velocity.
% 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
figure(1);clf
figure(2);clf
figure(3);clf
figure(4);clf
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,101);
omega1t=linspace(0,6*pi,181);
sn=[];
nterm=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,:)=ones(1,181)-cos(betaLsq(1,i)*omega1t/betaLsq(1,1));
    b(1,i)=2*alpha(1,i)/(betaL(1,i)^5);%expansion coef for constant load
end
xoverL=[-.06 xoverL];
framedata=[];
for k=1:181
    sum=zeros(1,101);
    for i=1:nterm
        sum=sum+b(1,i)*cs(i,k)*phi(i,:);
    end
    framedata=[framedata;sum];
end
lim=max(abs(framedata(30,:)));
framedata=[zeros(181,1) framedata];
figure(1)
h1=axes('YDir','reverse');
hold on
axis([-.2 1.2 -.5*lim 1.2*lim])
hold on
xpatch1=[0 0 -.1 -.08 -.09 0];
ypatch1=[.2 -.2 -.15 -.05 .05 .2]*lim/.4;
patch(xpatch1,ypatch1,'r')
hold on
xlabel('Dimensionless Distance, x/L')
hold on
xpatch2=[0 0 1 1 0];
ypatch2=[-.0035 -.05 -.05 -.0035 -.0035];
P1=patch(xpatch2,ypatch2,'c','EraseMode','background');
texthandl3=text(.2,-.23*lim,'Load f_0 applied at t = 0');
ylabel('Dimensionless Displacement, y(x,t)EI/f_0L^4')
S=plot(xoverL,framedata(2,:),'k','EraseMode','xor','LineWidth',[3.5]);
texthandl=text(.1,-.45*lim,'Press Enter to Animate One Frame at a Time');
pause
set(texthandl3,'String',' ');
set(P1,'Visible','off');
for k=4:4:60
    set(S,'Ydata',framedata(k,:));
    pause
end
set(texthandl,'String','Press Enter to Continue');
figure(2)%plot y(t) for 10 values of x
h2=axes('YDir','reverse');
hold on
axis([0 6*pi -.2*lim 1.2*lim]);
hold on
for j=2:10:102
    plot(omega1t,framedata(:,j))
    hold on
end
xlabel('Dimensionless Time, \omega_1t')
ylabel('Dimensionless Displacement, y(x,t)EI/f_0L^4')
text(1,-.06*lim,'Press Enter to Continue')
pause
%now do animation
figure(3)%plot animation
h3=axes('YDir','reverse');
hold on
axis([-.2 1.2 -.5*lim 1.2*lim])
hold on
patch(xpatch1,ypatch1,'r')
hold on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Displacement, y(x,t)EI/f_0L^4')
hold on
L=plot(xoverL, zeros(1,102),'k','EraseMode','xor','LineWidth',[3.5]);
hold off
texthandl1=text(.1,-.4*lim,'Press Enter to Set Initial Condition');
pause
set(L,'Ydata',framedata(1,:));
set(texthandl1,'String','Press Enter to Animate');
pause
set(texthandl1,'String','                       ');
for i=2:181
    set(L,'Ydata',framedata(i,:));
    pause(.1)
end
hold on
set(texthandl1,'String','Press Enter to Continue');
pause
figure(4)%3-D y vs. x and t
h4=axes('ZDir','reverse');
hold on
axis([0 1 0 6*pi -.5*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',-36)
ylabel('Dimensionless Time, \omega_1t','Rotation',11)
zlabel('Dimensionless Temperature, y(x,t)EI/f_0L^4')
view(60,30)
grid on
text(0,5,-.08,'Press Enter to Continue')
pause