% An elastic bar fixed at the left end and free at the right end with
% initial deformation of constant strain and no initial velocity.
% Displacement at the at the free end is U0.  Solution via the Fourier
% series solution to the one dimensional 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)=2*(((1/fac)^2)*sin(fac));
end
for i=2:40
    omegaoveromega1=[omegaoveromega1 2*i-1];
end
for i=1:nterm
  cs1=cos(omegaoveromega1(1,i)*omega1t); 
  cs=[cs;cs1];
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

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

figure(2);clf;
axis([0 4*pi -1.2 1.2])
hold on
text(6,1.1,'x/L=')
box on
for j=1:10:101
    plot(omega1t,framedata(:,j))
    hold on
    text(6,framedata(1,j),num2str(xoverL(1,j)))
    hold on
end
xlabel('Dimensionless Time, \omega_1t')
ylabel('Dimensionless Displacement, u(x,t))/U_0')
text(4,-1,'Press Enter to Continue')
pause
hold on
hold off
%now do animation
figure(3);clf;
axis([-.2 1.2 -1.2 1.2])
box on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Displacement, u(x,t)/U_0')
hold on
L=plot(xoverL, zeros(1,101),'k','EraseMode','xor');%,'LineWidth',2);
hold off
texthandl1=text(-.1,1.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
set(texthandl1,'String','Press Enter to Continue');
pause
hold on
hold off

figure(4);clf;
axis([-.2 1.2 -1 1]);
ypatch=[-.5 .5 .3 -.1 -.35 -.5];
xpatch=[0 0 -.15 -.08 -.17 0];
patch(xpatch, ypatch,'r');
hold on
box on
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(.3,.8,'Press Enter to Set Initial Condition');
xlabel('Distance Along the Bar, x/L')
hold on
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

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
    pause(.04)
end
set(t1,'String', 'Press Enter to Continue')
pause
hold on
hold off