% Sep 4, 2012
%
% This program illustrates Alpert quadrature for discretizing
% BIEs in two dimensions.
%
% We solve Laplace and Helmholtz BVPs on a domain exterior to a cornered/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
%
% 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 Alpert correction
%   's6'    Laplace single layer potential with  6th order Alpert correction
%   'su'    Laplace single layer potential with 10th order Alpert correction
%   'st'    Laplace single layer potential with 16th order Alpert correction
%   'ds'    Laplace single + double layer potential with 6th order Alpert
%   'h2'    Helmholtz double layer potential with 2nd order Alpert correction
%   'h6'    Helmholtz double layer potential with 6th order Alpert correction   
%   'hu'    Helmholtz combined field [single + i*kh*double], with 10th order Alpert
%   'ht'    Helmholtz combined field [single + i*kh*double], with 16th order Alpert


function main09_Alpert()
% rng('default')
% rng(5)

main_laplace_exterior
main_helmholtz_exterior

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function main_laplace_exterior

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

%%% Set the geometry.
flag_geom    = 'star';

fprintf(1,'flag_pot  err_submat       errmax       errmsq      int_qq  sum(qqsrc)\n');
FLAGS = ['s2';'s6';'su';'st';'ds'];

for ipot = 1:size(FLAGS,1)
  flag_pot = FLAGS(ipot,:);
  
  %%% Create the "correction matrix" and the full matrix.
  [A_corr,C,xxsrc,xxtrg,curvelen] = LOCAL_construct_A_corr(ntot,nsrc,ntrg,kh,flag_geom,flag_pot);
  [A,     ~,    ~,    ~,       ~] = LOCAL_construct_A_diag(ntot,nsrc,ntrg,kh,flag_geom,flag_pot);
  params = [curvelen,kh];
  
  %%% 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);
 
  % 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,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

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

return

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

function main_helmholtz_exterior

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

%%% Set the geometry.
flag_geom    = 'star';

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,C,xxsrc,xxtrg,curvelen] = LOCAL_construct_A_corr(ntot,nsrc,ntrg,kh,flag_geom,flag_pot);
  [A,     ~,    ~,    ~,       ~] = LOCAL_construct_A_diag(ntot,nsrc,ntrg,kh,flag_geom,flag_pot);
  params = [curvelen,kh];
  
  %%% 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);

  % 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,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

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

return

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

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

ntot = length(tt);

if strcmp(flag_geom,'star')
  r        = 0.3;
  k        = 5;
  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,C,xxint,xxext,curvelen] = LOCAL_construct_A_corr(ntot,nint,next,kh,flag_geom,flag_pot)

tt       = linspace(0,2*pi*(1 - 1/ntot),ntot);
[C,curvelen,xxint,xxext] = LOCAL_construct_cont(tt,nint,next,flag_geom);
h        = curvelen/ntot;


%%% Check what order corrections are required,
%%% and prepare the relevant index vectors for A_corr.
if ((strcmp(flag_pot,'s2')) || (strcmp(flag_pot,'h2')))
  [xx, ww, nskip] = alpert_log(2);
  ninterp = 2 + 3;

elseif ((strcmp(flag_pot,'s6')) || (strcmp(flag_pot,'ds')) ||...
       (strcmp(flag_pot,'h6')))
  [xx, ww, nskip] = alpert_log(6); 
  ninterp = 6 + 3;

elseif ((strcmp(flag_pot,'su')) || (strcmp(flag_pot,'hu')))
  [xx, ww, nskip] = alpert_log(10);
  ninterp = 10 + 3;

elseif((strcmp(flag_pot,'st')) || (strcmp(flag_pot,'ht')))  
  [xx, ww, nskip] = alpert_log(16);
  ninterp = 16 + 3;
end

xx = [-xx(end: -1: 1), xx];
ww = [ ww(end: -1: 1), ww];

s_g1   = C(1,:).'*ones(1, length(xx));
s_g2   = C(4,:).'*ones(1, length(xx));
t      = repmat(tt,1,length(xx));
x      = repmat(xx',1,ntot);
x      = reshape(x',ntot*length(xx),1)';
spt    = t+x*h;
spt    = mod(spt,curvelen);
[C_corr,~,~,~] = LOCAL_construct_cont(spt,nint,next,flag_geom);

t_g1  = reshape(C_corr(1,:)',ntot,length(xx));
t_g2  = reshape(C_corr(4,:)',ntot,length(xx));
t_dg1 = reshape(C_corr(2,:)',ntot,length(xx));
t_dg2 = reshape(C_corr(5,:)',ntot,length(xx));
dt    = sqrt(C_corr(2,:).^2 + C_corr(5,:).^2);
dt    = reshape(dt',ntot,length(xx));

kr = LOCAL_eval_kernel(s_g1,t_g1,s_g2,t_g2,t_dg1,t_dg2,kh,flag_pot);
kr = kr.*dt*h;
  
ind1 = floor(max(xx) - ninterp/2 + 1) + (0 : ninterp - 1);
ind2 = floor(min(xx) - ninterp/2 + 1) + (0 : ninterp - 1);
ind  = union(ind2,ind1);
II = repmat(1:ntot,[1,length(ind)]);
JJ = II + kron(ind,ones(1,ntot));
JJ = mod(JJ-1,ntot)+1;

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((strcmp(flag_pot,'s2'))|| (strcmp(flag_pot,'s6')) ||...
   (strcmp(flag_pot,'su')) || (strcmp(flag_pot,'st')))
    aa = -(h/(4*pi))*log(dd1.*dd1 + dd2.*dd2).*speed;
elseif(strcmp(flag_pot,'ds'))           
    ddsq = dd1.*dd1 + dd2.*dd2;  
    aa   = (-h/(2*pi))*((nn1.*dd1 + nn2.*dd2)./ddsq + 0.5*log(ddsq)).*speed; 
elseif((strcmp(flag_pot,'h2')) || (strcmp(flag_pot,'h6')))    
    ima = sqrt(-1);  
    dd  = sqrt(dd1.*dd1 + dd2.*dd2);
    aa  = h*(nn1.*dd1 + nn2.*dd2).*(1./dd).*(-kh*besselh(1, kh*dd))*ima/4.*speed;
elseif((strcmp(flag_pot,'hu')) || (strcmp(flag_pot,'ht')))  
    ima = sqrt(-1);      
    dd  = sqrt(dd1.*dd1 + dd2.*dd2);
    aa  = h*(((nn1.*dd1 + nn2.*dd2)./dd).*(-kh*besselh(1, kh*dd))*ima/4 - ...
                   ima*kh*besselh(0,kh*dd)                              ).*speed;
end

A_corr = sparse(II,JJ,aa,ntot,ntot);

kk = zeros(ntot,length([-nskip+1:nskip-1]));
indskip = [];
for p=1:ntot
  kk(p,:) = mod(p+[-nskip+1:nskip-1]-1,ntot)+1;
  [~,a1,~] = intersect(1:ntot,kk(p,:));
  indskip = [indskip,a1+(p-1)*ntot];
end

A_corr(indskip) = 0;

for i = 1:length(xx)
  ind = floor(xx(i) - ninterp/2 + 1) + (0 : ninterp - 1); 
  L   = lagrangeInterp(ind, xx(i));
  ii = repmat(1:ntot, [1, ninterp]);
  jj  = ii + kron(ind, ones(1, ntot));
  jj  = mod(jj-1, ntot) + 1;
  aakr = ww(i)*repmat(kr(:, i).', [1, ninterp]).*kron(L, ones(1, ntot));
  iii = sub2ind(size(A_corr), ii, jj);
  A_corr(iii) = A_corr(iii) + aakr;
end


if((flag_pot ~= 's2') |(flag_pot ~= 's6')|...
   (flag_pot ~= 'su') |(flag_pot ~= 'st'))
  A_corr = 0.5*speye(size(A_corr)) + A_corr; 
end

return



function [A,C,xxint,xxext,curvelen] = LOCAL_construct_A_diag(ntot,nint,next,kh,flag_geom,flag_pot)

tt       = linspace(0,2*pi*(1 - 1/ntot),ntot);
[C,curvelen,xxint,xxext] = LOCAL_construct_cont(tt,nint,next,flag_geom);
h        = curvelen/ntot;

if ((strcmp(flag_pot,'s2')) || (strcmp(flag_pot,'h2')))    
  [xx, ww, nskip] = alpert_log(2);
  ninterp = 2 + 3;

elseif ((strcmp(flag_pot,'s6')) || (strcmp(flag_pot,'ds')) ||...
       (strcmp(flag_pot,'h6')))
  [xx, ww, nskip] = alpert_log(6); 
  ninterp = 6 + 3;

elseif ((strcmp(flag_pot,'su')) || (strcmp(flag_pot,'hu')))
  [xx, ww, nskip] = alpert_log(10);
  ninterp = 10 + 3;

elseif((strcmp(flag_pot,'st')) || (strcmp(flag_pot,'ht')))  
  [xx, ww, nskip] = alpert_log(16);
  ninterp = 16 + 3;
end

xx = [-xx(end: -1: 1), xx];
ww = [ ww(end: -1: 1), ww];

[Y_g1,   X_g1   ] = meshgrid(C(1,:), C(1,:));
[Y_g2,   X_g2   ] = meshgrid(C(4,:), C(4,:));
[Y_dg1,  ~      ] = meshgrid(C(2,:), C(2,:));
[Y_dg2,  ~      ] = meshgrid(C(5,:), C(5,:));
dd = sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2);

A = LOCAL_eval_kernel(X_g1,Y_g1,X_g2,Y_g2,Y_dg1,Y_dg2,kh,flag_pot);
A = A.*dd*h;

for j = 1:ntot
  A(j, mod(j+[-nskip+1:nskip-1]-1, ntot)+1) = 0;  
end

s_g1   = C(1,:).'*ones(1, length(xx));
s_g2   = C(4,:).'*ones(1, length(xx));
t      = repmat(tt,1,length(xx));
x      = repmat(xx',1,ntot);
x      = reshape(x',ntot*length(xx),1)';
spt    = t+x*h;
spt    = mod(spt,curvelen);
[C_corr,~,~,~] = LOCAL_construct_cont(spt,nint,next,flag_geom);

t_g1  = reshape(C_corr(1,:)',ntot,length(xx));
t_g2  = reshape(C_corr(4,:)',ntot,length(xx));
t_dg1 = reshape(C_corr(2,:)',ntot,length(xx));
t_dg2 = reshape(C_corr(5,:)',ntot,length(xx));
dt    = sqrt(C_corr(2,:).^2 + C_corr(5,:).^2);
dt    = reshape(dt',ntot,length(xx));

kr = LOCAL_eval_kernel(s_g1,t_g1,s_g2,t_g2,t_dg1,t_dg2,kh,flag_pot);
kr = kr.*dt*h;

for i = 1:length(xx)
  ind = floor(xx(i) - ninterp/2 + 1) + (0 : ninterp - 1); % example index = [-2, -1, 0, 1, 2]
  L   = lagrangeInterp(ind, xx(i));
  ii  = repmat(1:ntot, [1, ninterp]);
  jj  = ii + kron(ind, ones(1, ntot));
  jj  = mod(jj-1, ntot) + 1;
  iii = sub2ind(size(A), ii, jj);
  A(iii) = A(iii) + ww(i)*repmat(kr(:, i).', [1, ninterp]).*kron(L, ones(1, ntot));
end

if((flag_pot ~= 's2') | (flag_pot ~= 's6') |...
   (flag_pot ~= 'su') |(flag_pot ~= 'st'))
  A = 0.5*speye(size(A)) + A; 
end

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';
ii  = II(:)';
jj  = JJ(:)';


%%% Check what order corrections are required,
%%% and prepare the relevant index vectors for A_corr.
%%% Find the index vectors marking the different types of entries.
aa      = reshape(full(A_corr(ind1,ind2)).',n1*n2,1).';
indnear = find(aa);       % Entries stored in A_corr
indreg  = LOCAL_index_complement(indnear,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((strcmp(flag_pot,'s2'))|| (strcmp(flag_pot,'s6')) ||...
   (strcmp(flag_pot,'su')) || (strcmp(flag_pot,'st')))
  aa(indreg)     = -(h/(4*pi))*log(dd1.*dd1 + dd2.*dd2).*speed;
elseif(strcmp(flag_pot,'ds'))           
  ddsq           = dd1.*dd1 + dd2.*dd2;  
  aa(indreg)     = (-h/(2*pi))*((nn1.*dd1 + nn2.*dd2)./ddsq + 0.5*log(ddsq)).*speed;      
elseif((strcmp(flag_pot,'h2')) || (strcmp(flag_pot,'h6')))
  ima            = sqrt(-1);  
  dd             = sqrt(dd1.*dd1 + dd2.*dd2);
  aa(indreg)     = h*(nn1.*dd1 + nn2.*dd2).*(1./dd).*(-kh*besselh(1, kh*dd))*ima/4.*speed;
elseif((strcmp(flag_pot,'hu')) || (strcmp(flag_pot,'ht')))
  ima            = sqrt(-1);      
  dd             = sqrt(dd1.*dd1 + dd2.*dd2);
  aa(indreg)     = h*(((nn1.*dd1 + nn2.*dd2)./dd).*(-kh*besselh(1, kh*dd))*ima/4 - ...
                   ima*kh*besselh(0,kh*dd)                              ).*speed;
end

A12 = reshape(aa.',n2,n1).';

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_eval_kernel(X_g1,Y_g1,X_g2,Y_g2,Y_dg1,Y_dg2,kh,flag_pot)

if((strcmp(flag_pot,'s2'))|| (strcmp(flag_pot,'s6')) ||...
   (strcmp(flag_pot,'su')) || (strcmp(flag_pot,'st'))) 
    
  A  = -1/(4*pi)*log((Y_g1 - X_g1).^2 + (Y_g2 - X_g2).^2);

elseif(strcmp(flag_pot,'ds'))    
    
  nn1  = (Y_dg2./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
  nn2  = (-Y_dg1./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
  ddsq = (Y_g1 - X_g1).^2 + (Y_g2 - X_g2).^2;
  A    = (-1/(2*pi))*((nn1.*(Y_g1 - X_g1) + nn2.*(Y_g2 - X_g2))./ddsq + 0.5*log(ddsq));
    
elseif((strcmp(flag_pot,'h2')) || (strcmp(flag_pot,'h6')))
    
  ima = sqrt(-1);  
  nn1 = ( Y_dg2./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
  nn2 = (-Y_dg1./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
  dd  = sqrt((Y_g1 - X_g1).^2 + (Y_g2 - X_g2).^2);
  A   = (nn1.*(Y_g1 - X_g1) + nn2.*(Y_g2 - X_g2)).*(1./dd).*(-kh*besselh(1, kh*dd))*ima/4;
        
elseif((strcmp(flag_pot,'hu')) || (strcmp(flag_pot,'ht')))
    
  ima = sqrt(-1);  
  nn1 = ( Y_dg2./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
  nn2 = (-Y_dg1./sqrt(Y_dg1.*Y_dg1 + Y_dg2.*Y_dg2));
  dd  = sqrt((Y_g1 - X_g1).^2 + (Y_g2 - X_g2).^2);
  A   = (nn1.*(Y_g1 - X_g1) + nn2.*(Y_g2 - X_g2)).*(1./dd).*(-kh*besselh(1, kh*dd))*ima/4 ...
        - ima*kh*besselh(0, kh*dd);
   
else
  fprintf(1,'This option for the layer potential is not implemented.\n');
  keyboard
end
 
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/(4*pi))*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, uu, potname, params)

curvelen = params(1);
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 (strcmp(potname,'ds'))
  ddsq  = dd1.*dd1 + dd2.*dd2;
  ee    = (-h/(2*pi))*((nn1.*dd1 + nn2.*dd2)./(ddsq) + 0.5*log(ddsq)).*speed;
  
elseif ((strcmp(potname,'s2'))|| (strcmp(potname,'s6')) ||...
   (strcmp(potname,'su')) || (strcmp(potname,'st')))
  ee    = -(h/(4*pi))*log(dd1.*dd1 + dd2.*dd2).*speed;
  
elseif ((strcmp(potname,'h2')) || (strcmp(potname,'h6')))
  ima   = sqrt(-1);  
  dd    = sqrt(dd1.*dd1 + dd2.*dd2);
  ee    = h*((nn1.*dd1 + nn2.*dd2)./dd).*(-kh*besselh(1,kh*dd))*ima/4.*speed;
  
elseif ((strcmp(potname,'hu')) || (strcmp(potname,'ht')))
  ima   = sqrt(-1);  
  dd    = sqrt(dd1.*dd1 + dd2.*dd2);
  ee    = h*(((nn1.*dd1 + nn2.*dd2)./dd).*(-kh*besselh(1, kh*dd))*ima/4 - ...
             ima*kh*besselh(0,kh*dd)                              ).*speed;
end

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

return


%%%%%%%%% Auxillury function%%%%%%%%%
function [xx, ww, nskip] = alpert_log(order)

switch order
    case 2 
    xx = 1.591549430918953D-01;
    ww = 5.000000000000000D-01;
	nskip = 1;
    
    case 6
    xx = [4.004884194926570D-03, 7.745655373336686D-02, ... 
          3.972849993523248D-01, 1.075673352915104D+00, ...
          2.003796927111872D+00];  
    ww = [1.671879691147102D-02, 1.636958371447360D-01, ...
          4.981856569770637D-01, 8.372266245578912D-01, ...
          9.841730844088381D-01];
	nskip = 3;
    
    case 10
    xx = [1.175089381227308D-03, 1.877034129831289D-02, ...
		  9.686468391426860D-02, 3.004818668002884D-01, ...  
		  6.901331557173356D-01, 1.293695738083659D+00, ...  
		  2.090187729798780D+00, 3.016719313149212D+00, ...  
		  4.001369747872486D+00, 5.000025661793423D+00];
    ww = [4.560746882084207D-03, 3.810606322384757D-02, ...
          1.293864997289512D-01, 2.884360381408835D-01, ...
          4.958111914344961D-01, 7.077154600594529D-01, ...
          8.741924365285083D-01, 9.661361986515218D-01, ...
          9.957887866078700D-01, 9.998665787423845D-01];
    nskip = 6;
    
    case 16
    xx = [8.371529832014113E-04, 1.239382725542637E-02,...  
		  6.009290785739468E-02, 1.805991249601928E-01,...  
		  4.142832599028031E-01, 7.964747731112430E-01,...  
		  1.348993882467059E+00, 2.073471660264395E+00,...  
		  2.947904939031494E+00, 3.928129252248612E+00,...  
		  4.957203086563112E+00, 5.986360113977494E+00,...  
		  6.997957704791519E+00, 7.999888757524622E+00....  
		  8.999998754306120E+00];

    ww = [3.190919086626234E-03, 2.423621380426338E-02, ...
         7.740135521653088E-02, 1.704889420286369E-01, ...
         3.029123478511309E-01,...
         4.652220834914617E-01, 6.401489637096768E-01, ...
         8.051212946181061E-01, 9.362411945698647E-01, ...
         1.014359775369075E+00, 1.035167721053657E+00, ...
         1.020308624984610E+00, 1.004798397441514E+00, ...
         1.000395017352309E+00, 1.000007149422537E+00];
	nskip = 10;
   
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function L = lagrangeInterp(xx, t)
% L is a row vector calculating the coefficients of lagrange interpolation
% at t. xx are data

n = length(xx); % here xx = index = [-2, -1, 0, 1, 2]
denom = xx'*ones(1, n) - ones(n, 1)*xx;
numer = t*ones(n, n) - ones(n, 1)*xx;
for i = 1:n
  denom(i, i) = 1;
  numer(i, i) = 1;
end
temp = numer./denom;
L = prod(temp, 2);
L = L';








