% Elastic bar fixed at the left end and free at the right end with
% suddenly applied constant force F with no initial deformation and 
% no initial deformation rate.  Solution via the Fourier series solution
% to the wave equation.
% 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);
xoverL=linspace(0,1,101);
omega1t=linspace(0,4*pi,401);
cs=[];
nterm=40;%40 modes
omegaoveromega1=[1];
for i=1:nterm
    fac=((2*i)-1)*pi/2;
    sine(i,:)=sin(fac*xoverL);
    B(1,i)=((-1)^(i+1))/((2*i-1)^2);
end
for i=2:40
    omegaoveromega1=[omegaoveromega1 2*i-1];
end
for i=1:nterm
  term=ones(1,401)-cos(omegaoveromega1(1,i)*omega1t); 
  cs=[cs;term];
end
framedata=[];

for k=1:401
    sum=zeros(1,101);
    for i=1:nterm
        sum=sum+B(1,i)*cs(i,k)*sine(i,:);
    end
    framedata=[framedata;sum];
end
framedata=(8/(pi^2))*framedata;

figure(1);clf;
axis([-.2 1.2 -.2 2.2])
hold on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Displacement, AEu(x,t)/FL')
S=plot(xoverL,framedata(1,:),'k','EraseMode','xor');%,'LineWidth',2);
text(-.1,2.1,'Press Enter to Continue')
box on
pause
for k=11:10:201
    set(S,'Ydata',framedata(k,:));
    text(-.1,2.1,'Press Enter to Continue')
    pause
end
hold off

figure(2);clf;
axis([0 4*pi -.2 2.2])
hold on
text(9,2.1,'x/L=')
box on
for j=1:10:101
    plot(omega1t,framedata(:,j))
    hold on
    text(9,framedata(301,j),num2str(xoverL(1,j)))
    hold on
end
xlabel('Dimensionless Time, \omega_1t')
ylabel('Dimensionless Displacement, AEu(x,t)/FL')
text(6,8,'Press Enter to Continue')
%axis([0 6 -1.2 1.2]);
pause
hold off

%now do animation
figure(3);clf;
axis([-.2 1.2 -.2 2.2])
hold on
box on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Displacement, AEu(x,t)/FL')
hold on
L=plot(xoverL, zeros(1,101),'k','EraseMode','xor');%,'LineWidth',2);
hold off
texthandl1=text(-.1,2.1,'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:401
    set(L,'Ydata',framedata(i,:))
    pause(.04)
end
hold on
%plot(xoverL,framedata(1,:),'k')
%hold on
set(texthandl1,'String','Press Enter to Continue');
pause
hold off

figure(4);clf;
xarrow=[0 .18 .16 .16 .18];
yarrow=[0 0 .02 -.02 0];
axis([-.2 1.4 -1 1]);
ypatch=[-.5 .5 .3 -.1 -.35 -.5];
xpatch=[0 0 -.15 -.04 -.17 0];
patch(xpatch, ypatch,'r');
hold on
box on
xlabel('Dimensionless Distance, x/L');
xbar=[0 1 1 0];
ybar=[.2 .2 -.2 -.2];
bar=plot(xbar, ybar,'Color', 'k','LineWidth',[2],'EraseMode','xor');
ytic=[-.18 .18];
for i=1:19
    tic(1,i)=plot([i*.05 i*.05],ytic,'EraseMode','xor');
end

t1=text(.5,.8,'Press Enter to Set Initial Condition');
pause
set(bar, 'Xdata', [0 1+.1*framedata(1,101) 1+.1*framedata(1,101) 0]);
for i=1:19
    set(tic(1,i), 'Xdata',[i*.05+.1*framedata(1,i*5+1) i*.05+.1*framedata(1,i*5+1)]);
end
arrow=plot(xarrow+ones(1,5),yarrow,'LineWidth',[2]);
text(1.22,.1,'F');
set(t1,'String','Press Enter to Animate               ');
pause
set(t1,'String','                                      ');
for k=2:401
    set(bar,'Xdata', [0 1+.1*framedata(k,101) 1+.1*framedata(k,101) 0]);
    for i=1:19
        g=i*.05+.1*framedata(k,i*5+1);
        set(tic(1,i), 'Xdata',[ g g]);
    end
    set(arrow,'Xdata',xarrow+(0.1*framedata(k,101)+1)*ones(1,5))
    pause(.04)
end
set(t1,'String', 'Press Enter to Continue')
pause