% Level curves r = 1 for the first 21 exponential sums (Figure 2)
global n
nmin = 1;  nmax = 21;  tol = 3.e-8;
axis equal;  hold on; % arc = [];
options = odeset('RelTol',tol,'AbsTol',tol);
for n = nmin:nmax,
   sf = 0.94574*(n +.5*log(2*pi*n) + 3); 
% 0.94574=(1+pi/2)*exp(-1)
   [s,z] = ode45('level4', [0 sf], [0; 0], options);
   indices = find(z(:, 2) >= 0); 
   l = length(indices);  l1 = l + 1;
   w = [-z(l1, 2), z(l, 2)]; w = w/sum(w); 
   z(l1,:) = w*z(l:l1,:); 
   plot(z(1:l1, 1), z(1:l1, 2));
% approximate arclengths and bounds sf
%  arc = [arc; w*s(l:l1), sf]; 
end; 

