clear all
clc
close all

syms rho real positive
syms phi real positive

z = rho*(cos(phi)+1i*sin(phi));
p = 1+z+0.5*z^2;
q = simplify(real(p*conj(p)));

dPhi = 2*pi/1440;
Phi = 0:dPhi:2*pi-dPhi;

hf = figure; hold on; axis equal; grid on;
for j=1:length(Phi)
  Q = subs(q,'phi',Phi(j));
  RHO{j} = solve(Q-1, 'Real', true);
  figure(hf);
  for k=1:length(RHO{j})
    plot(RHO{j}(k)*cos(Phi(j)),RHO{j}(k)*sin(Phi(j)),'k.');
  end
end