% Lines of constant phase for the 10th exponential sum (Fig 3)
%
global n phi
n = 10; tol = 1.e-5; sf = 1.5;
clf; axis([-6 6 -1 8]); hold on;
r = roots(1./gamma(n + 1:-1:1));
indices = find(imag(r) >= 0); zero = r(indices)
options= odeset('RelTol',tol,'AbsTol',tol);
for k = 1: length(zero),
   z0 = zero(k);
   for phi = -7/8*pi:pi/8:pi,
      [s,z] = ode45('phase', [0 sf], z0, options);
      plot(z);
   end;
end;

