% Aug 27, 2012
%
% This program illustrates Kapur-Rokhlin quadrature for discretizing
% BIEs in two dimensions.
%
% We solve Laplace and Helmholtz BVPs on a domain exterior to a smooth
% contour. Dirichlet boundary data is prescribed.
%
% There are four routines for evaluating the matrix:
%
% These functions evaluated entries "on-the-fly".
%    LOCAL_construct_A_diag
%    LOCAL_construct_A_offd
%
% These functions rely on a pre-computation of the corrected entries:
%     LOCAL_construct_A_corr
%     LOCAL_A_submat
%
% Several different kernels are tested:
%   's2'   Laplace single layer potential with  2nd order K-R correction
%   's6'   Laplace single layer potential with  6th order K-R correction
%   'su'   Laplace single layer potential with 10th order K-R correction
%   'dr'   Laplace double layer potential + point source (smooth kernel) 
%   'ds'   Laplace single + double layer potential with 6th order K-R
%   'h2'   Helmholtz double layer potential with 2nd order K-R correction
%   'h6'   Helmholtz double layer potential with 6th order K-R correction
%   'ht'   Helmholtz combined field [single + i*kh*double], with  6th order K-R
%   'hu'   Helmholtz combined field [single + i*kh*double], with 10th order K-R

function main_KR05()
rng('default')
rng(7)

main_laplace_exterior
main_helmholtz_exterior

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function main_laplace_exterior

%%% Set some parameters.
ntot = 200;
nsrc = 3;
ntrg = 20;
kh   = 0; % Set the Helmholtz parameter to zero to signal "Laplace problem".

%%% Set the geometry.
flag_geom    = 'star';
[C,curvelen,xxsrc,xxtrg] = LOCAL_construct_cont(ntot,flag_geom,nsrc,ntrg);
params   = [curvelen,kh];
plot(C(1,[1:ntot,1]), C(4,[1:ntot,1]    ),'k.-',...
     xxsrc(1,:), xxsrc(2,:),'r.',...
     xxtrg(1,:), xxtrg(2,:),'b.')
legend('Contour C','Source points xxt','Target points xxs')
axis equal

%%% Compute the Dirichlet data and the solution.
qqsrc  = randn(nsrc,1);
uu_dir = LOCAL_ref_field(C([1,4],:),xxsrc,qqsrc,params);
uu_ref = LOCAL_ref_field(xxtrg,     xxsrc,qqsrc,params);

fprintf(1,'flag_pot  err_submat       errmax       errmsq      int_qq  sum(qqsrc)\n')
FLAGS = ['s2';'s6';'su';'ds';'dr'];
for ipot = 1:size(FLAGS,1)
  flag_pot = FLAGS(ipot,:);
  
  %%% Create the "correction matrix" and the full matrix.
  A_corr = LOCAL_construct_A_corr(C,flag_pot,params);
  A      = LOCAL_construct_A_diag(C,1:ntot,flag_pot,params);

  % Check that the fuction for evaluating off-diagonal blocks works
  ind1 = [24:26,33];
  ind2 = [33:36,75];
  err_submat = max(max(abs(A(ind1,ind2) - LOCAL_A_submat(C,ind1,ind2,A_corr,flag_pot,params))));
    
  %%% Solve a reference problem.
  qq = A\uu_dir;

  %%% Check the error between the potential generated by the integral equation
  %%% and the reference solution.
  uu = LOCAL_evalpot(xxtrg,C,curvelen,qq,flag_pot,params);
  errmax = max(abs(uu - uu_ref));
  errmsq = sqrt(sum((abs(uu - uu_ref).^2))/length(uu));

  %%% Evaluate the integral of the boundary potential:
  dg = (curvelen/ntot)*sqrt(C(2,:).^2 + C(5,:).^2);
  int_qq = dg*qq;

  fprintf(1,'      %2s  %10.2e   %10.2e   %10.2e   %9.6f   %9.6f\n',...
          flag_pot,err_submat,errmax,errmsq,int_qq,sum(qqsrc))
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function main_helmholtz_exterior

%%% Set some parameters.
ntot = 200;
nsrc = 3;
ntrg = 20;
kh   = 2;

%%% Set the geometry.
flag_geom    = 'star';
[C,curvelen,xxsrc,xxtrg] = LOCAL_construct_cont(ntot,flag_geom,nsrc,ntrg);
params   = [curvelen,kh];
plot(C(1,[1:ntot,1]), C(4,[1:ntot,1]    ),'k.-',...
     xxsrc(1,:), xxsrc(2,:),'r.',...
     xxtrg(1,:), xxtrg(2,:),'b.')
legend('Contour C','Source points xxt','Target points xxs')
axis equal

%%% Compute the Dirichlet data and the solution.
qqsrc  = randn(nsrc,1);
uu_dir = LOCAL_ref_field(C([1,4],:),xxsrc,qqsrc,params);
uu_ref = LOCAL_ref_field(xxtrg,     xxsrc,qqsrc,params);

fprintf(1,'flag_pot  err_submat       errmax       errmsq')
fprintf(1,'                        int_qq                   sum(qqsrc)\n')
FLAGS = ['h2';'h6';'ht';'hu'];
for ipot = 1:size(FLAGS,1)
  flag_pot = FLAGS(ipot,:);
  
  %%% Create the "correction matrix" and the full matrix.
  A_corr = LOCAL_construct_A_corr(C,flag_pot,params);
  A      = LOCAL_construct_A_diag(C,1:ntot,flag_pot,params);

  % Check that the fuction for evaluating off-diagonal blocks works
  ind1 = [24:26,33];
  ind2 = [33:36,75];
  err_submat = max(max(abs(A(ind1,ind2) - LOCAL_A_submat(C,ind1,ind2,A_corr,flag_pot,params))));
    
  %%% Solve a reference problem.
  qq = A\uu_dir;

  %%% Check the error between the potential generated by the integral equation
  %%% and the reference solution.
  uu     = LOCAL_evalpot(xxtrg,C,curvelen,qq,flag_pot,params);
  errmax = max(abs(uu - uu_ref));
  errmsq = sqrt(sum((abs(uu - uu_ref).^2))/length(uu));

  %%% Evaluate the integral of the boundary potential:
  dg     = (curvelen/ntot)*sqrt(C(2,:).^2 + C(5,:).^2);
  int_qq = dg*qq;

  fprintf(1,'      %2s  %10.2e   %10.2e   %10.2e   (%9.6f + 1i* %9.6f)   (%9.6f +1i* %9.6f)\n',...
          flag_pot,err_submat,errmax,errmsq,...
          real(int_qq),imag(int_qq),real(sum(qqsrc)),imag(sum(qqsrc)))
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [C,curvelen,xxint,xxext] = ...
             LOCAL_construct_cont(ntot,flag_geom,nint,next)

if strcmp(flag_geom,'star')
  r        = 0.3;
  k        = 5;
  tt       = linspace(0,2*pi*(1 - 1/ntot),ntot);
  C        = zeros(6,ntot );
  C(1,:)   =   1.5*cos(tt) + (r/2)*            cos((k+1)*tt) + (r/2)*            cos((k-1)*tt);
  C(2,:)   = - 1.5*sin(tt) - (r/2)*(k+1)*      sin((k+1)*tt) - (r/2)*(k-1)*      sin((k-1)*tt);
  C(3,:)   = - 1.5*cos(tt) - (r/2)*(k+1)*(k+1)*cos((k+1)*tt) - (r/2)*(k-1)*(k-1)*cos((k-1)*tt);
  C(4,:)   =       sin(tt) + (r/2)*            sin((k+1)*tt) - (r/2)*            sin((k-1)*tt);
  C(5,:)   =       cos(tt) + (r/2)*(k+1)*      cos((k+1)*tt) - (r/2)*(k-1)*      cos((k-1)*tt);
  C(6,:)   = -     sin(tt) - (r/2)*(k+1)*(k+1)*sin((k+1)*tt) + (r/2)*(k-1)*(k-1)*sin((k-1)*tt);
  curvelen = 2*pi;
  rmin     = sqrt(min(C(1,:).^2 + C(4,:).^2));
  rmax     = sqrt(max(C(1,:).^2 + C(4,:).^2));
  ttint    = 2*pi*rand(1,nint);
  xxint    = 0.6*[rmin*cos(ttint);rmin*sin(ttint)];
  ttext    = 2*pi*rand(1,next);
  xxext    = 1.6*[rmax*cos(ttext);rmax*sin(ttext)];
else
  fprintf(1,'This option for the geometry is not implemented.\n')
  keyboard
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function constructs a sparse matrix "A_corr" that holds all
% "corrected" matrix entries corresponding to near-diagonal entries.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A_corr = LOCAL_construct_A_corr(C,flag_pot,params)

%%% Extract various parameters.
curvelen = params(1);
kh       = params(2);
ntot     = size(C,2);
h        = curvelen/ntot;

%%% Check what order corrections are required,
%%% and prepare the relevant index vectors for A_corr.
if ((flag_pot == 's2') | (flag_pot == 'h2'))
  MU2 = [ 0.7518812338640025 + 0.1073866830872157e1; ...
         -0.7225370982867850 - 0.6032109664493744];
  JJ = ones(4,1)*(1:ntot);
  II = [(-2):(-1),1:2]' * ones(1,ntot) + JJ;
  II = mod(II-1,ntot)+1;
  KR = MU2([2,1,1,2]) * ones(1,ntot);
  ii = II(:)';
  jj = JJ(:)';
  kr = KR(:)';
elseif ((flag_pot == 's6') | (flag_pot == 'ds')  | (flag_pot == 'h6') | (flag_pot == 'ht'))
  MU6 = [ 0.2051970990601250e1 + 0.2915391987686505e1;...
         -0.7407035584542865e1 - 0.8797979464048396e1;...
          0.1219590847580216e2 + 0.1365562914252423e2;...
         -0.1064623987147282e2 - 0.1157975479644601e2;...
          0.4799117710681772e1 + 0.5130987287355766e1;...
         -0.8837770983721025   - 0.9342187797694916];
  JJ = ones(12,1)*(1:ntot);
  II = [(-6):(-1),1:6]' * ones(1,ntot) + JJ;
  II = mod(II-1,ntot)+1;
  KR = [MU6(end:(-1):1);MU6] * ones(1,ntot);
  ii = II(:)';
  jj = JJ(:)';
  kr = KR(:)';
elseif ((flag_pot == 'hu') | (flag_pot == 'su'))
  MU10 = [ 0.3256353919777872D+01 + 0.4576078100790908D+01;...
          -0.2096116396850468D+02 - 0.2469045273524281D+02;...
           0.6872858265408605D+02 + 0.7648830198138171D+02;...
          -0.1393153744796911D+03 - 0.1508194558089468D+03;...
           0.1874446431742073D+03 + 0.1996415730837827D+03;...
          -0.1715855846429547D+03 - 0.1807965537141134D+03;...
           0.1061953812152787D+03 + 0.1110467735366555D+03;...
          -0.4269031893958787D+02 - 0.4438764193424203D+02;...
           0.1009036069527147D+02 + 0.1044548196545488D+02;...
          -0.1066655310499552D+01 - 0.1100328792904271D+01];
  JJ = ones(20,1)*(1:ntot);
  II = [(-10):(-1),1:10]' * ones(1,ntot) + JJ;
  II = mod(II-1,ntot)+1;
  KR = [MU10(end:(-1):1);MU10] * ones(1,ntot);
  ii = II(:)';
  jj = JJ(:)';
  kr = KR(:)';
else
  A_corr = sparse(ntot,ntot);
  return
end

speed = sqrt(C(2,jj).*C(2,jj) + C(5,jj).*C(5,jj));
nn1   =  C(5,jj)./speed;
nn2   = -C(2,jj)./speed;
dd1   = C(1,jj) - C(1,ii);
dd2   = C(4,jj) - C(4,ii);

if ((flag_pot == 's2') | (flag_pot == 's6') | (flag_pot == 'su')) % Laplace: single 
  aa   = -(h/(4*pi))*log(dd1.*dd1 + dd2.*dd2).*speed;
elseif (flag_pot == 'ds') % Laplace: single + double
  ddsq  = dd1.*dd1 + dd2.*dd2;
  aa    = (-h/(2*pi))*((nn1.*dd1 + nn2.*dd2)./(ddsq) + 0.5*log(ddsq)).*speed;
elseif ((flag_pot == 'h2') | (flag_pot == 'h6')) % Helmhotz: double 
  dd    = sqrt(dd1.*dd1 + dd2.*dd2);
  aa    = h*(nn1.*dd1 + nn2.*dd2).*(1./dd).*(-kh*besselh(1, kh*dd)).*speed;
elseif ((flag_pot == 'ht') | (flag_pot == 'hu')) % Helmoltz: double + i*kh*single
  dd    = sqrt(dd1.*dd1 + dd2.*dd2);
  aa    = h*((nn1.*dd1 + nn2.*dd2).*(1./dd).*(-kh*besselh(1, kh*dd)) + ...
             1i*kh*besselh(0,kh*dd)                                 ).*speed;
end

aakr   = (1+kr).*aa;
A_corr = sparse(ii,jj,aakr,ntot,ntot);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A12 = LOCAL_A_submat(C,ind1,ind2,A_corr,flag_pot,params)

%%% Extract various parameters.
n1       = length(ind1);
n2       = length(ind2);
curvelen = params(1);
kh       = params(2);
ntot     = size(C,2);
h        = curvelen/ntot;

%%% Create the index vectors.
II  = reshape(ind1,n1,1) * ones(1,n2);
JJ  = ones(n1,1) * reshape(ind2,1,n2);
ii  = II(:)';
jj  = JJ(:)';

%%% Find the index vectors marking the different types of entries.
aa      = reshape(full(A_corr(ind1,ind2)),1,n1*n2);
indnear = find(aa);      % Entries stored in A_corr
inddiag = find(ii==jj);  % Diagonal entries (must be zeroed out explicitly for K-R).
indreg  = LOCAL_index_complement([indnear,inddiag],n1*n2);

%%% Compute some vectors used in almost all kernels.
speed = sqrt(C(2,jj(indreg)).*C(2,jj(indreg)) + ...
             C(5,jj(indreg)).*C(5,jj(indreg)));
nn1   =  C(5,jj(indreg))./speed;
nn2   = -C(2,jj(indreg))./speed;
dd1   = C(1,jj(indreg)) - C(1,ii(indreg));
dd2   = C(4,jj(indreg)) - C(4,ii(indreg));

if ((flag_pot == 's2') | (flag_pot == 's6') | (flag_pot == 'su')) % single layer laplace
  aa(inddiag) = zeros(1,length(inddiag));
  aa(indreg)  = -(h/(4*pi))*log(dd1.*dd1 + dd2.*dd2).*speed;
elseif (flag_pot == 'dr') % The regular double layer potential for Laplace.
  %%% Construct the regular entries.
  aa(indreg) = (-h/(2*pi))*((nn1.*dd1 + nn2.*dd2)./(dd1.*dd1 + dd2.*dd2)).*speed;
  %%% Construct the diagonal entries.
  dg1  = C(2,ii(inddiag));
  dg2  = C(5,ii(inddiag));
  ddg1 = C(3,ii(inddiag));
  ddg2 = C(6,ii(inddiag));
  aa(inddiag) = (h/(4*pi))*(-ddg2.*dg1 + ddg1.*dg2)./(dg1.^2 + dg2.^2) + ...
                0.5*ones(1,length(inddiag));
  %%% Add the contribution from a mono-pole at xxc.
  xxc = [0;0];
  aa  = aa + (-h/(4*pi))*log((C(1,ii) - xxc(1)).^2 + (C(4,ii) - xxc(2)).^2).*sqrt(C(2,jj).^2 + C(5,jj).^2);
elseif (flag_pot == 'ds') % The double + single layer for Laplace.
  ddsq        = dd1.*dd1 + dd2.*dd2;
  aa(indreg)  = (-h/(2*pi))*((nn1.*dd1 + nn2.*dd2)./ddsq + 0.5*log(ddsq)).*speed;
  aa(inddiag) = 0.5*ones(1,length(inddiag));
elseif ((flag_pot == 'h2') | (flag_pot == 'h6')) % Double layer Helmholtz
  dd          = sqrt(dd1.*dd1 + dd2.*dd2);
  aa(indreg)  = h*(nn1.*dd1 + nn2.*dd2).*(1./dd).*(-kh*besselh(1, kh*dd)).*speed;
  aa(inddiag) = -2*1i*ones(1,length(inddiag));
elseif ((flag_pot == 'ht') | (flag_pot == 'hu')) % Helmholtz: double + i * k * single
  dd          = sqrt(dd1.*dd1 + dd2.*dd2);
  aa(indreg)  = h*(((nn1.*dd1 + nn2.*dd2)./dd).*(-kh*besselh(1, kh*dd)) + ...
                   1i*kh*besselh(0,kh*dd)                              ).*speed;
  aa(inddiag) = -2*1i*ones(1,length(inddiag));
end

%%% Reshape the vector to the proper form.
A12 = reshape(aa,n1,n2);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Given an index vector ind1, find ind2 = (1:n) - ind1.
% (In other words, [ind1,ind2] forms a permutation of (1:n).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ind2 = LOCAL_index_complement(ind1,n)

indtmp       = 1:n;
indtmp(ind1) = 2*n*ones(1,length(ind1));
indtmp       = sort(indtmp);
ind2         = indtmp(1:(n-length(ind1)));

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_construct_A_diag(C,ind,flag_pot,params)

curvelen = params(1);
kh       = params(2);
ntot     = size(C,2);
h        = curvelen/ntot;
nloc     = length(ind);

%%% Create some index vectors.
II = ind' * ones(1,nloc);
ii = reshape(II, 1,nloc*nloc);
jj = reshape(II',1,nloc*nloc);
inddiag = (1:nloc) + nloc*(0:(nloc-1));
clear II

%%% Check what order corrections are required,
%%% and prepare the relevant index vectors for A_corr.
if ((flag_pot == 's2') | (flag_pot == 'h2'))
  MU2 = [ 0.7518812338640025 + 0.1073866830872157e1; ...
         -0.7225370982867850 - 0.6032109664493744];
  tt = abs(mod(abs(ii - jj)+2,ntot)-2); % This is the distance from the diagonal.
  indnear = find(tt.*(tt < 3)); % Find the near-field.
  kr_corr = (1 + MU2(tt(indnear))');
elseif ((flag_pot == 's6') | (flag_pot == 'ds')  | (flag_pot == 'h6')  | (flag_pot == 'ht'))
  MU6 = [ 0.2051970990601250e1 + 0.2915391987686505e1;...
         -0.7407035584542865e1 - 0.8797979464048396e1;...
          0.1219590847580216e2 + 0.1365562914252423e2;...
         -0.1064623987147282e2 - 0.1157975479644601e2;...
          0.4799117710681772e1 + 0.5130987287355766e1;...
         -0.8837770983721025   - 0.9342187797694916];
  tt = abs(mod(abs(ii - jj)+6,ntot)-6); % This is the distance from the diagonal.
  indnear = find(tt.*(tt < 7)); % Find the near-field.
  kr_corr = 1 + MU6(tt(indnear))';
elseif ((flag_pot == 'su') | (flag_pot == 'hu'))
  MU10 = [ 0.3256353919777872D+01 + 0.4576078100790908D+01;...
          -0.2096116396850468D+02 - 0.2469045273524281D+02;...
           0.6872858265408605D+02 + 0.7648830198138171D+02;...
          -0.1393153744796911D+03 - 0.1508194558089468D+03;...
           0.1874446431742073D+03 + 0.1996415730837827D+03;...
          -0.1715855846429547D+03 - 0.1807965537141134D+03;...
           0.1061953812152787D+03 + 0.1110467735366555D+03;...
          -0.4269031893958787D+02 - 0.4438764193424203D+02;...
           0.1009036069527147D+02 + 0.1044548196545488D+02;...
          -0.1066655310499552D+01 - 0.1100328792904271D+01];
  tt = abs(mod(abs(ii - jj)+10,ntot)-10); % This is the distance from the diagonal.
  indnear = find(tt.*(tt < 11)); % Find the near-field.
  kr_corr = 1 + MU10(tt(indnear))';
end

speed = sqrt(C(2,jj).*C(2,jj) + C(5,jj).*C(5,jj));
nn1   =  C(5,jj)./speed;
nn2   = -C(2,jj)./speed;
dd1   = C(1,jj) - C(1,ii);
dd2   = C(4,jj) - C(4,ii);

if ((flag_pot == 's2') | (flag_pot == 's6') | (flag_pot == 'su')) % single layer Laplace
  ddsq = dd1.*dd1 + dd2.*dd2;
  ddsq(inddiag) = ones(1,nloc);
  aa = -(h/(4*pi))*log(ddsq).*speed;
  aa(indnear) = kr_corr.*aa(indnear);
elseif (flag_pot == 'dr') % This is the "real" double layer potential.
  %%% Construct the regular entries.
  ddsq  = dd1.*dd1 + dd2.*dd2;
  ddsq(inddiag) = ones(1,nloc);
  aa    = (-h/(2*pi))*((nn1.*dd1 + nn2.*dd2)./ddsq).*speed; 
  %%% Construct the diagonal entries.
  dg1  = C(2,ii(inddiag));
  dg2  = C(5,ii(inddiag));
  ddg1 = C(3,ii(inddiag));
  ddg2 = C(6,ii(inddiag));
  aa(inddiag) = (h/(4*pi))*(-ddg2.*dg1 + ddg1.*dg2)./(dg1.^2 + dg2.^2) + ...
                0.5*ones(1,length(inddiag));
  %%% Add the contribution from a monopole at xxc.
  xxc = [0;0];
  dd1 = C(1,ii) - xxc(1);
  dd2 = C(4,ii) - xxc(2);
  aa  = aa + (-h/(4*pi))*log(dd1.*dd1 + dd2.*dd2).*speed;
elseif (flag_pot == 'ds') % This is the "real" double layer potential + a single layer potential.
  %%% Construct the regular entries.
  ddsq  = dd1.*dd1 + dd2.*dd2;
  ddsq(inddiag) = ones(1,nloc);
  aa = (-h/(2*pi))*((nn1.*dd1 + nn2.*dd2)./(ddsq) + 0.5*log(ddsq)).*speed;
  aa(inddiag) = 0.5*ones(1,nloc);
  aa(indnear) = kr_corr.*aa(indnear);
elseif ((flag_pot == 'h2') | (flag_pot == 'h6')) % Double layer Helmholtz
  dd    = sqrt(dd1.*dd1 + dd2.*dd2);
  dd(inddiag) = ones(1,nloc);
  aa = h*(nn1.*dd1 + nn2.*dd2).*(1./dd).*(-kh*besselh(1, kh*dd)).*speed;
  aa(inddiag) = -2*1i*ones(1,nloc);
  aa(indnear) = kr_corr.*aa(indnear);
elseif ((flag_pot == 'ht') | (flag_pot == 'hu')) % Helmoltz kernel:  double +  i kh single
  dd    = sqrt(dd1.*dd1 + dd2.*dd2);
  dd(inddiag) = ones(1,nloc);
  aa = h*((nn1.*dd1 + nn2.*dd2).*(1./dd).*(-kh*besselh(1, kh*dd)) + ...
          1i*kh*besselh(0,kh*dd)                                 ).*speed;
  aa(inddiag) = -2*1i*ones(1,nloc);
  aa(indnear) = kr_corr.*aa(indnear);
end

A = full(sparse(ii,jj,aa));

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function uu = LOCAL_ref_field(xxtrg,xxsrc,qq,params)

nsrc = length(qq);
ntrg = size(xxtrg,2);
DD1  = xxtrg(1,:)'*ones(1,nsrc) - ones(ntrg,1)*xxsrc(1,:);
DD2  = xxtrg(2,:)'*ones(1,nsrc) - ones(ntrg,1)*xxsrc(2,:);
kh   = params(2);

if (kh == 0) % Laplace problem
  AA = -(1/(2*pi))*0.5*log(DD1.*DD1 + DD2.*DD2);
else % Helmholtz problem.
  AA = (1i/4)*besselh(0,kh*sqrt(DD1.*DD1 + DD2.*DD2));
end

uu = AA*reshape(qq,nsrc,1);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function vv = LOCAL_evalpot(xx, C, curvelen, uu, potname, params)

kh = params(2);
n  = size(C,2);
m  = size(xx,2);
h  = curvelen/n;

%%% Create index vectors.
ii = reshape((1:m)' * ones(1,n),1,m*n);
jj = reshape(ones(m,1) * (1:n) ,1,m*n);

%%% Compute some vectors used on almost all kernels.
speed = sqrt(C(2,jj).*C(2,jj) + C(5,jj).*C(5,jj));
nn1   =  C(5,jj)./speed;
nn2   = -C(2,jj)./speed;
dd1   =  C(1,jj) - xx(1,ii);
dd2   =  C(4,jj) - xx(2,ii);

if (potname == 'dr')
  %%% Construct the regular entries.
  ddsq  = dd1.*dd1 + dd2.*dd2;
  ee    = (-h/(2*pi))*((nn1.*dd1 + nn2.*dd2)./ddsq).*speed; 
  %%% Add the contribution from a monopole at xxc.
  xxc = [0;0];
  dd1 = xx(1,ii) - xxc(1);
  dd2 = xx(2,ii) - xxc(2);
  ee  = ee + (-h/(4*pi))*log(dd1.*dd1 + dd2.*dd2).*speed;
elseif (potname == 'ds')
  ddsq  = dd1.*dd1 + dd2.*dd2;
  ee    = (-h/(2*pi))*((nn1.*dd1 + nn2.*dd2)./(ddsq) + 0.5*log(ddsq)).*speed;
elseif ((potname == 's2') | (potname == 's6') | (potname == 'su'))
  ee    = -(h/(4*pi))*log(dd1.*dd1 + dd2.*dd2).*speed;
elseif ((potname == 'h2') | (potname == 'h6'))
  dd    = sqrt(dd1.*dd1 + dd2.*dd2);
  ee    = h*((nn1.*dd1 + nn2.*dd2)./dd).*(-kh*besselh(1,kh*dd)).*speed;
elseif ((potname == 'ht') | (potname == 'hu'))
  dd    = sqrt(dd1.*dd1 + dd2.*dd2);
  ee    = h*(((nn1.*dd1 + nn2.*dd2)./dd).*(-kh*besselh(1, kh*dd)) + ...
             1i*kh*besselh(0,kh*dd)                              ).*speed;
end

EVAL = full(sparse(ii,jj,ee));
vv   = EVAL*uu;

return