m = n/6;
k = 1;
for r = 0:m
  for j = 0:6*r
    if (r==0) | (j>0)
      a = max(r,1);
      x(k) = 2*r/m*cos(2*pi*j/(6*a));
      y(k) =   r/m*sin(2*pi*j/(6*a));
      z(k) = 0;
      c1(k) = eval1(x(k),y(k),xx,yy,mu);
      c2(k) = c1(k) - log( (2+x(k))^2 + (1+y(k))^2 );
      k = k+1;
    end
  end
end

% make a contour plot of c1 or c2
tri = delaunay(x,y);
trisurf(tri,x,y,z,c1);
shading interp;
view(0,90);
axis equal;

% plot the mesh on top of the contour plot
if (n<=96)
  hold on;
  h = trimesh(tri,x,y,z,z,'EdgeColor',[0 0 0],'FaceAlpha',0);
  hold off;
end

% add a legend
colorbar('horiz')