function p_roots = parabola_plot(a,b)
%function p_roots = parabola_plot(a,b)
%
% PARABOLA_PLOT: Plots parabolas of the form y = x^2 + ax + b,
%                calculates, the roots, and plots the roots
%                on the complex plane. 
%
% Input Parameters:
%    a, b    --  Coefficients of the parabola y = x^2 + ax + b.  May
%                be vectors to plot multiple parabolas simultaneously.
%
% Output Parameters:
%    p_roots --  Roots of the parabola
%

% Make sure a and b are not matrices
if ~any(size(a) == 1) | ~any(size(b)==1)
    error('a and b cannot be matrices.');
end

% Check input vector sizes
if size(a,1) < size(a,2)
    a = a';
end
if size(b,1) < size(b,2)
    b = b';
end
if ~all(size(b) == size(a))
    if all(size(a) == 1)
        a = repmat(a,size(b));
    elseif all(size(b) == 1)
        b = repmat(b,size(a));
    else
        error('a and b must have the same size.');
    end
end

coefs = [ones(size(b)) a b];

% Generate time axis
t = -2:.01:2;

% Calculate the roots for each parabola
pv = [];
p_roots = [];
for i=1:length(b)
    pv = [pv polyval(coefs(i,:),t)'];
    p_roots = [p_roots; roots(coefs(i,:))'];
end

% Plot the parabolas
subplot(1,2,1)
plot(t,pv);
grid on;
hold on;
plot(t,0,'k-');
hold off;
title('y = x^2 + ax + b');
xlabel('x');
ylabel('y');

% Plot the roots on the complex plane with the unit circle 
subplot(1,2,2)
c = -pi:.01:pi;
h = plot(cos(c),sin(c),'k-',real(p_roots'),imag(p_roots'),'o');
grid on
axis image
set(h,'MarkerSize',10);
title('Root positions');
xlabel('Real(z)');
ylabel('Imag(z)');
ax = axis;
h = line(ax(1:2),[0 0]);
set(h,'color',[0 0 0]);
h = line([0 0],ax(3:4));
set(h,'color',[0 0 0]);
