% Phasor demo
%
% Tobin Fricke <tobin.fricke@ligo.org> 2011-06-08

% modulation depth
%m = pi/4;
m = 1;

% vector of (normalized) frequencies to include
%freqs = [0 1 -1 2 -2 3 -3 4 -4];
freqs = [0 1 -1 2 -2];

% initial amplitudes of the phasors
amp = besselj(freqs, m);

set(0, 'defaultlinelinewidth', 5);

% Set up the axes
cla reset;
axis equal;
xlim([-1 1.5]);
ylim([-1 1]);

% Draw the unit circle
rectangle('Position', [-1 -1 2 2], 'Curvature', [1 1])

% make some lines
L = [];
L(1) = line([0], [0]);
L(2) = line([0], [0], 'color', 'red');

save_to_disk = 0;

% how many frames to render?
if save_to_disk,
    n_frames = 200;
else
    n_frames = inf;
end

% how much the 1-Omega should rotate in one frame of the animation?
if isinf(n_frames)
    dphi = m/20;
else
    dphi = (2*pi)/(n_frames + 1);
end

frameno = 0;
ret_amp = [];
while (frameno < n_frames),
    frameno = frameno + 1;

    % Add up the phasors cumulatively
    seg = cumsum([0 amp]);

    % Draw both the collection of phasors and the resultant
    set(L(1), 'xdata', real(seg), 'ydata', imag(seg));
    set(L(2), 'xdata', real(seg([1,end])), 'ydata', imag(seg([1,end])));

	ret_amp = [ret_amp, abs(seg([end]))];

    % Update the display
    drawnow;

    if (save_to_disk),
        print('-dpng', '-r72', sprintf('pm%03d.png', frameno));
    else
        pause(0.02);
    end

    % Advance the phasors for the next frame
    amp = amp .* exp(1i * dphi * freqs);
end