% Plucked string animation via the D'Alembert 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
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
m=fix(100*aoverL);
k=[1:m];
f1=(0.5/m)*k;
k=[m+1:99];
f2=(0.5/(100-m))*[-k+100*ones(1,99-m)];
f0=[f1 f2];
f3=-fliplr(f0);
y0=[0 f0 0 f3 0 f0 0 f3 0 f0 0 f3 0 f0 0 f3 0 f0 0];
xoverL=linspace(0,1,101);
%patch coordinates
xp1=.5*[0 0 -.1 -.3 -.2 -.32 -.2 0 0];
yp1=[0 .4 .3 .2 .1 0 -.2 -.4 0];
xp2=-fliplr(xp1)+ones(1,9);
yp2=-fliplr(yp1);
time=linspace(0,2,201);
xoverL1=linspace(-.3,1.3,161);
ydata=2*y0(1,401:501);
for k=1:200
    ys=y0(1,401+2*k:501+2*k)+y0(1,401-2*k:501-2*k);
    ydata=[ydata ;ys];
end
figure(1);clf;
clf
set(gca,'Box','on')
patch(xp1,yp1,'r')
hold on
patch(xp2,yp2,'r')
hold on
plot([0 1],[0 0],'k')
hold on
k=0;
L1=plot(xoverL1,y0(1,401-30+2*k:501+30+2*k),'r','EraseMode','xor');
hold on
L2=plot(xoverL1,y0(1,401-30-2*k:501+30-2*k),'k','EraseMode','xor');
hold on
L3=plot(xoverL,ydata(k+1,:),'b','EraseMode','xor');
hold on
axis([-.3 1.3 -1.2 1.2]);
texthandl2=text(.3,-1.01,'Press Enter to Continue');
xlabel('Distance , x/L')
ylabel('Displacement, y(x,t)/y0')
pause
for k=2:2:100
    set(L1,'Ydata',y0(1,401-30+2*k:501+30+2*k));
    set(L2,'Ydata',y0(1,401-30-2*k:501+30-2*k));
    set(L3,'Ydata',ydata(k+1,:));
    
    hold on
    xlabel('Distance, x/L')
    ylabel('Displacement, y(x,t)/y_0')
    set(texthandl2,'String','Press Enter to Continue');
    axis([-.3 1.3 -1.2 1.2]);
    pause
end

figure(2);clf;
for i=11:10:91
    plot(time,ydata(:,i))
    hold on
    text(.3,ydata(1,i), ['x/L = ' num2str(xoverL(1,i))]);
end
axis([0 2,-1.1 1.1])
ylabel('Dimensionless Displacement, y(x,t)/y_0')
xlabel('Dimensionless Time, t/T_1')
text(1.2,1,'Press Enter to Continue')
pause

figure(3);clf;
clf
set(gca,'Box','on')
patch(xp1,yp1,'r')
axis([0 1 -1.2 1.2]);
hold on
xlabel('Distance, x/L')
ylabel('Displacement, y(x,t)/y_0')
axis([-.35 1.35 -2 2])
hold on
plot([aoverL aoverL],[-.1 -1.5])
hold on
text(aoverL-.05,-1.7,'a/L')
hold on
patch(xp2,yp2,'r')
plot(xoverL,zeros(1,101))
hold on
texthandl=text(0,1.5,'Press Enter to Set Initial Condition');
pause
ys=2*y0(1,401:501);
L=plot(xoverL,ys,'k','EraseMode','xor','LineWidth',[2.5]);
ydata=[ys];
set(texthandl,'String','Press Enter to Animate');
xlabel('Distance , x/L')
ylabel('Displacement, y(x,t)/y_0')
pause
set(texthandl,'String','   ');
for k=1:200
    ys=y0(1,401+2*k:501+2*k)+y0(1,401-2*k:501-2*k);
    set(L,'Ydata',ys)
    pause(.05)
    ydata=[ydata ;ys];
end
set(texthandl,'String','Press Enter to Continue');
pause

figure(4);clf;
time1=linspace(0,2,201);
[T,X]=meshgrid(time1,xoverL);
mesh(T,X,ydata')
%colormap([0 0 0]);
view([30 15])
axis([0 2 0 1 -1.3 1.3])
text(.2, 0,-1.1 ,'Dimensionless Time, t/T_1','rotation',-7)
text(2, .02 ,-1.1,'Dimensionless Distance, x/L','rotation',20)
zlabel('Dimensionless Displacement, y(x,t)/y_0')
text(.5,1,1,'Press Enter to Quit','Rotation',-7)
pause