% traveling.m
%
% This Matlab script creates a simple waveguide string model that can
% be used to view traveling-wave components as they move along the
% string.
%
% by Gary P. Scavone
% McGill University, 2016-2020.

clear all;
fprintf('\nSelect waveform type:\n\n');
fprintf('  1. Pluck\n');
fprintf('  2. Sinusoid\n');
fprintf('  3. Impulse\n');
fprintf('  4. Helmholtz\n\n');
waveform = input('Enter option number [1]: ');
if isempty(waveform)
  waveform = 1;
end
if (waveform < 1) || (waveform > 4)
  waveform = 1;
end

if waveform == 1
  pp = input('\nSet pluck point (0-1.0) [0.5]:\n');
  if isempty(pp)
    pp = 0.5;
  end
  if pp < 0.0 || pp > 1.0
    pp = 0.5;  % pluck point between bridge (0) and nut (1.0)
  end
end

if waveform == 2
  harmonic = input('\nSet harmonic number [1]:\n');
  if isempty(harmonic)
    harmonic = 1;
  end
  if harmonic < 1
    harmonic = 1;
  end
end

anechoic = input('\nInfinite string length? Y/[N]:\n', 's');
if anechoic ~= 'Y'
  anechoic = 'N';
end

if waveform ~= 4
  gradientOrder = input('\nPlot displacement (0), velocity (1) or acceleration (2)? [0], 1, or 2:\n');
  if isempty(gradientOrder) || gradientOrder < 0 || gradientOrder > 2
    gradientOrder = 0;
  end
else
  gradientOrder = 0;
end

pausetime = input('\nPause time between time steps (default=1s): \n');
if isempty(pausetime) || pausetime < 0
  pausetime = 1;
end

% ******* Initialize delay lines ******* %

% Create the initial physical shape of the string. Then each traveling-wave
% component is 1/2 the initial shape. To avoid double endpoints, each
% delayline is length D-1 and the "missing" element is added when plotting
% from the corresponding element on the opposite delayline. We flip
% left-right the upper delayline elements so we can use a single pointer
% incrementing in the same direction for both delaylines.

D = 100; % this is arbitrary
if waveform == 1 % triangular pluck
  ppsamp = round(pp * D);
  wshape = [(0:ppsamp-1)/(ppsamp-1), (D-ppsamp-1:-1:0)/(D-ppsamp)];

elseif waveform == 2 % sinusoidal
  wshape = sin(harmonic*pi*(0:D-1)/(D-1));

elseif waveform == 3 % pulse
  plen = round(0.2 * D);
  if ~mod(plen, 2)
    plen = plen+1;
  end
  wshape = [0.5 * sin(pi*(0:plen-1)/(plen-1)), zeros(1,D-plen)];
else % helmholtz
  wshape = linspace(1, 0, D);
end

if gradientOrder == 1
  ytext = 'Vertical Velocity';
  wshape = -D/2*gradient(wshape);
  dlines = 0.5 * [fliplr(wshape(2:D)); -wshape(1:D-1)];
elseif gradientOrder == 2
  ytext = 'Vertical Acceleration';
  wshape = -D/2*gradient(gradient(wshape));
  dlines = 0.5 * [fliplr(wshape(2:D)); -wshape(1:D-1)];
else
  ytext = 'Vertical Displacement';
  if waveform == 4
    wshape = wshape.^2;
    dlines = 0.5 * [fliplr(wshape(2:D)); -wshape(1:D-1)];
  else
    dlines = 0.5 * [fliplr(wshape(2:D)); wshape(1:D-1)];
  end
end

ptr = 1;

% ******* Plot initial string shape and traveling-wave components ******* %

figure(1);
h = gca;
x = (0:D-1)/(D-1);
upline = fliplr([dlines(1, ptr:D-1), dlines(1, 1:ptr-1), -dlines(2, ptr)]);
dwnline = [dlines(2, ptr:D-1), dlines(2, 1:ptr-1), -dlines(1, ptr)];
plot(x, upline + dwnline)
hold on
plot(x, upline, 'r')
plot(x, dwnline, 'g')
hold off
axis([0 1 -1.5 1.5])
grid
ylabel(ytext)
xlabel('Normalized String Length')
legend('Physical String', 'Traveling Wave to Right', 'Traveling Wave to Left')

disp('Paused ... hit any key to start and then close the figure to quit.');
pause

% *************************************** %
%                                         %
%             Delay Line(d)               %
%    |-------------------------------|    %
%    ^                                    %
%  pointer                                %
%                                         %
%     >>--- pointer increments --->>      %
%                                         %
% *************************************** %
%
% The pointer initially points to the delay line output.
% We can take the output and calculate a new input value
% which is placed where the output was taken from.  The
% pointer is then incremented and the process repeated.

% ******* Run Loop Start ******* %

while( ishandle(h) )

  temp = dlines(2, ptr);

  if anechoic == 'Y'
    dlines(2, ptr) = 0;
    dlines(1, ptr) = 0;
  else % rigid terminations and lossless
    dlines(2, ptr) = -dlines(1, ptr);
    dlines(1, ptr) = -temp;
  end
  
  % ****** Increment Pointers & Check Limits ****** %
  ptr = ptr + 1;
  if ptr >= D
    ptr = 1;
  end 

  % ****** Plot string shape and traveling-wave components  ****** %
  if anechoic == 'Y'
    upline = fliplr([dlines(1, ptr:D-1), dlines(1,1:ptr-1), 0]);
    dwnline = [dlines(2, ptr:D-1), dlines(2, 1:ptr-1), 0];
  else
    upline = fliplr([dlines(1, ptr:D-1), dlines(1, 1:ptr-1), -dlines(2, ptr)]);
    dwnline = [dlines(2, ptr:D-1), dlines(2, 1:ptr-1), -dlines(1, ptr)];
  end

  plot(x, upline + dwnline, 'b')
  hold on
  plot(x, dwnline, 'g')
  plot(x, upline, 'r')
  hold off
  axis([0 1 -1.5 1.5])
  grid
  ylabel(ytext)
  xlabel('Normalized String Length')
  legend('Physical String', 'Traveling Wave to Right', 'Traveling Wave to Left')
  
  % Uniform pause between each tick
  pause(pausetime)
end
