% solve interior/exterior BVP with double-layer potential: Laplace and Helmholtz
% barnett 10/31/08

w = 3; a = 0.3;
R = @(t) 1 + a*cos(w*t); Rp = @(t) -a*w*sin(w*t);   % wobbly shape
Rpp = @(t) -a*w*w*cos(w*t);
%for k=0.1:0.1:4 % ................ uncomment for k-sweep
k = 6;                       % k=0 laplace (int only), k>0 helmholtz (int,ext)
ext = 1;                                            % ext=0,1: interior,exterior
g = -1.3:0.02:1.3;                                  % viewing grid
%g = 1;                                             % single pt (for N-loop)
N = 50;                                             % quad pts
%for N=20:20:300 % -------------- uncomment for N-loop
verb = 0; tauone = 0; % verb=0:quiet, 1:show stuff. set tauone=1 to test tau=1
if k==0
  f = @(x1,x2) cos(x1).*exp(x2);                    % Dirichlet data, laplace
else   % H0 source at (0,2) or (0,.5) if need exterior case; Helmholtz
  f = @(x1,x2) 5i*besselh(0,k*sqrt(x1.^2 + (x2 - (2-1.5*ext)).^2));
end
   
if tauone, f = @(x1,x2) -ones(size(x1)); end      % interior val for tau=1
t = (1:N)/N*2*pi; sp = sqrt(Rp(t).^2 + R(t).^2);  % speed func
w = ones(size(t))*2*pi/N;                         % weights, peri trap quad
y = [R(t).*cos(t); R(t).*sin(t)];
% compute n_y * speed:
nysp = [0 1;-1 0]*[Rp(t).*cos(t) - R(t).*sin(t); Rp(t).*sin(t) + R(t).*cos(t)];
ypp = [-R(t).*cos(t) - 2*Rp(t).*sin(t) + Rpp(t).*cos(t); ...
       -R(t).*sin(t) + 2*Rp(t).*cos(t) + Rpp(t).*sin(t)];
ny = nysp./repmat(sp, [2 1]);           % unit vector normal n_y
kappa = -sum(ny.*ypp,1)./sp.^2;
if verb                                 % show geometry, normals
  figure; plot(y(1,:), y(2,:), 'k.'); hold on; l=0.1; axis equal;
  plot([y(1,:); y(1,:)+l*ny(1,:)], [y(2,:); y(2,:)+l*ny(2,:)],'-'); 
  figure; plot(t, kappa, '-');
end
% The next 9 lines are in lieu of a separate function for filling kernel...
y1 = repmat(y(1,:), [N 1]); y2 = repmat(y(2,:), [N 1]);
ny1 = repmat(ny(1,:), [N 1]); ny2 = repmat(ny(2,:), [N 1]);
r = sqrt((y1-y1').^2+(y2-y2').^2); r(diagind(r)) = 1; % dummy diag
costh = -(ny1.*(y1-y1') + ny2.*(y2-y2'))./r;
if k==0
  K = (1/2/pi)* costh ./ r;
else
  K = (k*1i/4)* costh .* besselh(1, k*r);
end
K(diagind(K)) = -kappa/(4*pi);
K = K * diag(w.*sp);                        % scale cols, weights and speed
if verb, figure; imagesc(t, t, real(K)); colorbar; drawnow; end
if ext, A = eye(N) + 2*K; else A = eye(N) - 2*K; end    % matrix
fprintf('N=%d, k=%g, A filled: cond # = %g\n', N, k, cond(A))
b = -2*f(y(1,:),y(2,:)).'; if ext, b = -b; end   % rhs
tau = A \ b;                                    % linear solve
if tauone, tau = ones(size(tau)); end              % use const tau

[xx1 xx2] = meshgrid(g);     % plot solution (slow)
xx = [xx1(:) xx2(:)]';
u = 0*xx1(:)'; n = numel(u);
for j=1:N           % sum density * weight * fund sol w/o prefac
  d = xx - repmat(y(:,j), [1 n]);           % array of displacements 2-by-n
  if k==0 
    u = u + tau(j) * w(j) * sum(repmat(nysp(:,j), [1 n]).*d,1) ./ sum(d.^2,1) / 2/pi;
  else
    r = sqrt(sum(d.^2,1));
    costh = sum(repmat(nysp(:,j), [1 n]).*d,1)./r; 
    u = u + tau(j) * w(j) * (k*1i/4) * costh .* besselh(1,k*r);
  end
end
u = reshape(u,size(xx1));          % make rect array

if numel(u)==1, fprintf('\tRe u = %.16g\n', real(u)); end
%end  % ----------- uncomment to end N loop

uex = f(xx1,xx2); if verb
figure; imagesc(g, g, real(uex)); caxis([-3 3]); axis equal; colormap(jet(256))
xlabel('x_1'); ylabel('x_2'); title('exact u'); colorbar;
figure; imagesc(g, g, real(u)); caxis([-3 3]); axis equal; colormap(jet(256))
xlabel('x_1'); ylabel('x_2'); title('DLP soln'); colorbar;
end
figure; imagesc(g, g, log10(abs(u-uex))); caxis([-16 0]); axis equal;
colormap(jet(256)); colorbar;
xlabel('x_1'); ylabel('x_2'); title('DLP soln: log_{10} error');
if numel(g)>1
  fprintf('err @ an interior pt = %g\n', abs((u(70,70)-uex(70,70))/uex(70,70)));
end