% solves string problem,
% u_tt - c^2*u_xx = 0 for 0<x<L, t>0,
% u(0,t) = u(L,t) = 0 for t >= 0,
% u(x,0) = u0(x) for 0 <= x <= L,
% using truncated Fourier series

% typical violin string
rho = 0.0006; % [kg/m]
Tension = 220; % [N = kg*m/s^2]
c = sqrt(Tension/rho); % [m^2/s^2]
L = 0.38; % [m]

N = 100; % number of terms in truncated Fourier series
numx = 8*N; % number of x values
numt = 300; % number of t values
nvec = (1:N); 
tvec = linspace(0,6*L/c,numt)'; % time lenght = three periods
x = linspace(0,L,numx)';

%%% initial position of string (a pluck in the middle)
eps = 0.1*L;
x1 = [linspace(0,L/2-eps,10),linspace(L/2+eps,L,10)];
beta = 1/10;
f = beta*(L/2 - abs(x1-L/2)); 
u0 = spline(x1,f,x);

%%% solve with truncated Fourier series
Phi = sin(pi/L*x*nvec);
T = cos(pi/L*c*tvec*nvec);
UU = (Phi*diag(Phi'*u0)*T')'*2/(numx-1);

%%% plot solution u(x,t) as a surface
figure(1), clf
%mesh(x,tvec,UU)
mesh(linspace(0,L,N),tvec(1:2:end),UU(1:2:end,1:8:end))
xlabel('x','fontsize',16)
ylabel('t','fontsize',16)
zlabel('u(x,t)','fontsize',16)
title(sprintf('string solution when c = %f',c),'fontsize',16)

%figure(3), clf
%plt = waterfall(linspace(0,1,N),tvec(1:4:end),UU(1:4:end,1:8:end));
%set(plt,'edgecolor','b')     
%xlabel('x','fontsize',16)
%ylabel('t','fontsize',16)
%zlabel('u(x,t)','fontsize',16)
%ylim([min(tvec) max(tvec)])
%title(sprintf('string solution when c = %f',c),'fontsize',16)


%%% animated plot of string position
figure(2), clf
ht = plot(x,UU(1,:));
axis([0 L -1.2*max(abs(u0)) 1.2*max(abs(u0))]) 
set(gca,'DataAspectRatio',[1 1 1])
xlabel('x','fontsize',16)
ylabel('t','fontsize',16)
title(sprintf('t = %5.3f',tvec(1)),'fontsize',16)
input('hit return to start\n');
for j = 2:numt
    set(ht,'YData',UU(j,:))
    title(sprintf('t = %5.4f',tvec(j)),'fontsize',16)
    drawnow
%    pause(0.03)
end