% Plucked string animation via the Fourier series solution to the wave 
% Equation. Input is the pluck location a/L, a fraction between zero 
% and one.
% 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
aoverL=-1;
while ((aoverL<=0)|(aoverL>=1))
    aoverL=input('Input the point where the string is plucked as fraction of the length (a/L) = ');
end
figure(1);clf
figure(2);clf
figure(3);clf
figure(4);clf
xoverL=linspace(0,1,101);
ctoverL=linspace(0,6,181);
sn=[];
nterm=40;%40 modes
for i=1:nterm
    sn(i,:)=sin(i*pi*xoverL);
    cs(i,:)=cos(i*pi*ctoverL);
    b(1,i)=2*sin(i*pi*aoverL)/(((i*pi)^2)*aoverL*(1-aoverL));
end
sn(:,101)=0;
framedata=[];
for k=1:181
    sum=zeros(1,101);
    for i=1:nterm
        sum=sum+b(1,i)*cs(i,k)*sn(i,:);
    end
    framedata=[framedata;sum];
end

figure(1)
xpatch1=[0 0 -.1 -.08 -.09 0];
ypatch1=[.2 -.2 -.15 -.05 .05 .2];
patch(xpatch1,ypatch1,'r')
xpatch2=1-fliplr(xpatch1);
ypatch2=fliplr(ypatch1);
patch(xpatch2,ypatch2,'r');
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Displacement, y(x,t)/y_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])
pause
for k=7:6:61
    set(S,'Ydata',framedata(k,:));
    axis([-.2 1.2 -1.2 1.2])
    text(-.1,1.1,'Press Enter to Continue')
    pause
end
figure(2)
axis([0 6 -1.2 1.2])
hold on
text(2.8,1.1,'x/L=')
for j=11:10:91
    plot(ctoverL,framedata(:,j))
    hold on
    text(2.8,framedata(1,j),num2str(xoverL(1,j)))
    hold on
end
xlabel('Dimensionless Time, ct/L')
ylabel('Dimensionless Displacement, y(x,t)/y_0')
text(.5,1.1,'Press Enter to Continue')
axis([0 6 -1.2 1.2]);
pause
%now do animation
figure(3)
patch(xpatch1,ypatch1,'r')
patch(xpatch2,ypatch2,'r');
axis([-.2 1.2 -1.2 1.2])
box on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Displacement, y(x,t)/y_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=1:181
    set(L,'Ydata',framedata(i,:))
    pause(.04)
end
hold on
set(texthandl1,'String','Press Enter to Continue');
pause

figure(4)
[X,Y]=meshgrid(xoverL,ctoverL);
mesh(X,Y,framedata)
xlabel('Dimensionless Distance, x/L','Rotation',-36)
ylabel('Dimensionless Time, ct/L','Rotation',11)
zlabel('Dimensionless Displacement, y(x,t)/y_0')
view(60,30)
axis([0 1 0 6 -1 1])