%Solves and animates the problem of transient motion of a uniform cantilever
%forced by an impulse of intensity I0 at location x=a.
%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);
aoverL=-1;
while ((aoverL<=0)|(aoverL>1))
    aoverL=input('Input the location of the  impulse, a/L = ');
end
phiofL=[2 -2 2 -2 2];
boverL=1-aoverL;
pt=(100*aoverL)+1;
betaL=[1.8751041 4.69409113 7.85475743 10.99554074 14.13716839];
gamma=(betaL/betaL(1,1)).^2;
betaL4=betaL.^4;
alpha=[.7340995 1.0184664 .9992245 1.000033553 .9999985501];
betan=betaL4/betaL4(1,1);
x1=linspace(0,1,101);
phiofx=[];
omega1t=linspace(0,8*pi,1001);
for i=1:5
    arg1=betaL(1,i)*x1;
    phiofx=[phiofx;cosh(arg1)-cos(arg1)-alpha(1,i)*(sinh(arg1)-sin(arg1))];
end
phiofa=phiofx(:,pt)';
Y=[];

for t=1:1001
    sum=zeros(1,101);
    for i=1:5
        sum=sum+(phiofa(1,i)*phiofx(i,:)/gamma(1,i))*sin(gamma(1,i)*omega1t(1,t));
    end
    Y=[Y;sum];
end
Y=[zeros(1,101);Y];
Y=[zeros(1002,1) Y];



figure(1);clf;% Plot y(x,t) at x/L=0.2,0.4,0.6,0.8,1
h3=axes('YDir','reverse');
hold on
axis([0 8*pi -5 5]);
hold on
for j=22:20:102
    plot(omega1t,Y(2:1002,j))
    hold on
end
for j=22:20:102
    plot([omega1t(1,431) 13],[Y(432,j) -(2+0.5*((j-2)/20))],'k');
    hold on
    text(13.2,-(2+0.5*((j-2)/20)),['x/L = ' num2str(x1(1,j-1))]);
    hold on
end
box on
xlabel('Dimensionless Time, \omega_1t')
ylabel('Beam Displacement, y(x,t)\muL\omega_1/I_0')
text(3,4,['a/L = ' num2str(aoverL)])
text(2,-4,'Press Enter to Continue');
pause

figure(2);clf;
h3=axes('YDir','reverse');
hold on
axis([-.2 1.3 -5 5]);
hold on
box on
xpatch1=[0 0 -.1 -.08 -.09 0];
ypatch1=[.2 -.2 -.15 -.05 .05 .2]*6;
xlabel('Dimensionless Distance, x/L')
ylabel('Beam Displacement, y(x,t)\muL\omega_1/I_0')
patch(xpatch1,ypatch1,'r')
hold on
plot([-.05 x1], Y(1,:),'k','LineWidth',[3.5]);
hold on
for i=13:12:253
    plot([-.05 x1], Y(i,:),'k','LineWidth',[3.5]);
    hold on
end
text(0,4,'y(x,0)=0');
text(0,4.5,['Impulse Location a/L = ' num2str(aoverL)])
text(0,-4,'Press Enter to Continue');
pause


figure(3);clf;% plot motion of point where load is applied
h3=axes('YDir','reverse');
hold on
axis([0 8*pi -5 5]);
hold on
plot(omega1t,Y(2:1002,pt+1));
xlabel('Dimensionless Time, \omega_1t')
ylabel('Forced Point Displacement, y(a,t)\muL\omega_1/I_0')
text(2,-4,'Press Enter to Continue')
text(2,4,['a/L = ' num2str(aoverL)])
box on
pause
hold off;

figure(4);clf;% Animate the beam
h3=axes('YDir','reverse');
hold on
axis([-.2 1.3 -5 5])
hold on
xpow=[1.1 1.02 1.05 .95 .88 .92 .9 1.01 1.04 1.03 1.1];
ypow=-10*[.1 .14 .17 .13 .15 .07 .04 .06 .05 .08 .1];
xpatch1=[0 0 -.1 -.08 -.09 0];
ypatch1=[.2 -.2 -.15 -.05 .05 .2]*6;
xlabel('Dimensionless Distance, x/L')
ylabel('Beam Displacement, y(x,t)\muL\omega_1/I_0')
patch(xpatch1,ypatch1,'r')
L=plot([-.05 x1], Y(1,:),'k','EraseMode','xor','LineWidth',[3.5]);
box on
texthandl1=text(0,-4,'Press Enter to Set Initial Condition');
text(0,4,'y(x,0)=0');
text(0,4.5,['Impulse Location a/L = ' num2str(aoverL)])
pause
xpow=xpow+(-1+x1(1,pt))*ones(1,11);
POW=patch(xpow,ypow,'m');
set(texthandl1,'String','Press Enter to Animate');
pause
set(POW,'Visible','off')
set(texthandl1,'String','                       ');
for i=3:1:1001
    set(L,'Ydata',Y(i,:))
    pause(.005)
end
set(texthandl1,'String','Press Enter to Continue');
pause
hold off;