function fourierViz(varargin)
%FOURIERVIZ Fourier visualization in Matlab
% 
%   By default, the fastest component has a frequency of 1 Hz

    % Handle input
    p = inputParser;
    p.addOptional('frequency', 1);
    p.addOptional('amplitude', 1);
    p.addOptional('phase', 0);
    p.addParamValue('baseMultiply', 2);
    p.addParamValue('export', false);
    p.addParamValue('filename', 'fourierVizLastRun.gif', @ischar);
    p.parse(varargin{:});
    p = p.Results;
    
    % Component frequency scale factor
    freqScaleFactor = max(p.frequency) / p.baseMultiply;
    
    % Normalize magnitudes
    p.amplitude = p.amplitude / sum(p.amplitude);
    
    % Figure
    doPlot = true;
    fh = figure('Position', [200 100 1600 800]);
    opts = struct;
    opts.margin = [0 0 0 0];
    opts.spaceBetweenPlots = 0;
    opts.nX = 1;
    opts.nY = 1;
    opts.iX = 1;
    opts.iY = 1; 
    mk.subplot(opts);
    hold on;
    set(gca,'xtick',[]);
    set(gca,'xticklabel',[]);
    set(gca,'ytick',[]);
    set(gca,'yticklabel',[]);
    xlim([-1.2 3.2]);
    ylim([-1.2 1.2]);
    
    % Time
    dt = 0.016;
    tick = 0;
    
    % Base frequency & period
    baseFrequency = sym(p.frequency(1)/freqScaleFactor);
    for ii = 2:numel(p.frequency)
        baseFrequency = gcd(baseFrequency, sym(p.frequency(ii)/freqScaleFactor));
    end
    baseFrequency = double(baseFrequency);
    basePeriod = 1/baseFrequency;
    
    % Circle
    xCircle = 0:2*pi/100:2*pi;
    
    % Oscilloscope
    osc = nan(6e2,1);
    oscilloscopePrefilled = false;
    
    % Colors
    clrGrey = struct('color', [0.8784 0.8784 0.8784]);
    clrGreen = struct('color', [0.1804 0.4902 0.1961]);
    clrBlue = struct('color', [0.0824 0.3961 0.7529]);
    clrRed = struct('color', [0.7765 0.1569 0.1569]);
    
    while doPlot
        tick = tick + 1;
        
        % Has the figure been closed
        if ~ishandle(fh)
            doPlot = false;
            continue;
        end
        
        % Update time
        t = dt * tick;
        
        % Update oscilloscopePrefilled
        if ~oscilloscopePrefilled && t>=basePeriod
            if ~any(isnan(osc))
                oscilloscopePrefilled = true;
            end
            tick = 1;
            t = dt * tick;
        end
        
        % Init plot coordinates
        x = 0;
        y = 0;
        
        % Plot
        if oscilloscopePrefilled
            cla;
            plot(cos(xCircle), sin(xCircle), '-k', clrGrey);
        end
        
        % Compute components
        for ii = 1:numel(p.frequency)
            if mod(ii,2) == 0
                clr = clrGreen;
            else
                clr = clrBlue;
            end
            alpha = ((t * 2*pi * p.frequency(ii)) / freqScaleFactor) + p.phase(ii);
            x2 = x + (p.amplitude(ii)*cos(alpha));
            y2 = y + (p.amplitude(ii)*sin(alpha));
            if oscilloscopePrefilled
                plot([x x2], [y y2], '-k', clr, 'LineWidth', 2);
                plot(x+cos(xCircle)*p.amplitude(ii), y+sin(xCircle)*p.amplitude(ii), '-k', clr, 'LineWidth', 2);
            end
            x = x2;
            y = y2;
        end
        osc = [y; osc(1:end-1)];
        if oscilloscopePrefilled
            plot(x, y, 'or', clrRed);
            plot([x, 1.5], [y y], '-r', clrRed, 'LineWidth', 1);
            plot(linspace(1.5, 3.2, numel(osc)), osc, '-r', clrRed, 'LineWidth', 2);
        end
        
        % Export to gif
        if oscilloscopePrefilled && p.export && t < basePeriod
            if ishandle(fh)
                frame = getframe(fh);
                im = frame2im(frame);
            end
            [imind, cm] = rgb2ind(im, 256);

            if tick == 1
                imwrite(imind, cm, p.filename, 'gif', 'Loopcount', inf, 'DelayTime', 0.016); 
            else
                imwrite(imind, cm, p.filename, 'gif', 'WriteMode', 'append', 'DelayTime', 0.016); 
            end
        end
        
        % Pause time
        if oscilloscopePrefilled
            pause(dt);
        end
    end

end