function m25_faraday()

%% m25_faraday  Play with the Faraday cage.
%    An example of a constrained quadratic programming problem.
%    See Chapman-Hewett-Trefethen, SIAM Review, September 2015.

%% Set the point charge and choose radius of wires:
clf, hold off, axis([-1.2 2.2 -1.7 1.7]), axis manual,
hold on, axis square
zs = 2;
LW = 'linewidth'; MS = 'markersize';
plot(real(zs),imag(zs),'.b',MS,16)
r = input('radius r? (e.g. 0.1 or 0.02) ');
circle = r*exp(1i*linspace(0,2*pi,30));

%% Input the points:
x = []; y = []; button = 1;
disp('input points with mouse, button >= 2 for final point')
while button == 1
  [xx,yy,button] = ginput(1);
  xyc = xx + 1i*yy + circle;
  x = [x; xx]; y = [y; yy]; h = fill(real(xyc),imag(xyc),'r');
  set(h,'edgecolor','r')
z = x + 1i*y; n = length(z);
A = -log(r)*eye(n);
c = ones(n,1); f = log(abs(z-zs));
for j = 1:n
  for i = 1:n
    if i~=j, A(i,j) = A(i,j) - log(abs(z(i)-z(j))); end
  end
end
B = [A c; c' 0]; xs = B\[f;0]; lambda = xs(end); q = xs(1:n);
xp = linspace(-1.2,2.2,70); yp = linspace(-1.7,1.7,70);
[xxp,yyp] = meshgrid(xp,yp); zzp = xxp+1i*yyp;
set(gca,'xtick',[],'ytick',[]), uu = ueval(zzp);
if exist('h2'), delete(h2), end;
[foo,h2] = contour(xxp,yyp,uu,-.31:.02:1.5,LW,.6); colormap([0 0 0])
for j = 1:n
  xyc = z(j) + circle;
  h = fill(real(xyc),imag(xyc),'r');
  set(h,'edgecolor','r')
end
end

function u = ueval(zz)
u = log(abs(zz-zs));
for k = 1:n, u = u + q(k)*log(abs(zz-z(k))); end
end

function u = ueval2(zz)
u = log(abs(zz-zs));
for k = 1:n2, u = u + dz*q2(k)*log(abs(zz-z2(k))); end
end

end