% March, 2012
%
% This program sets up a simple external scattering problem for Helmholtz.

function main_develop

rng('default')
rng(0)

LOCAL_test_single_cont()
LOCAL_test_two_cont()

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test function for scattering from a single contour.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_test_single_cont()

fprintf(1,'\n=== Running a test on a single contour.\n\n')

ntot = 200;

% Set the geometry.
xxc = [1;3];  % The center of the contour.
[C,curvelen,xx_int,xx_ext] = LOCAL_get_geometry_star(ntot,xxc);
h   = curvelen/ntot;
plot(C(1,[1:ntot,1]), C(4,[1:ntot,1]),'k',...
     xx_int(1,:), xx_int(2,:),'r.',...
     xx_ext(1,:), xx_ext(2,:),'b.',...
     'MarkerSize',30)
legend('Contour C','xxint','xxext')
axis equal

% Set the Helmholtz parameter.
kh       = 2;
lencont  = 2*max(sqrt((C(1,:)-xxc(1)).^2 + (C(4,:)-xxc(2)).^2));  % The diameter of the contour.
fprintf(1,'The diameter of the contour in wavelengths = %7.2f\n',kh*lencont/(2*pi))
fprintf(1,'Number of points per wave-length           = %7.2f\n',ntot*2*pi/(curvelen*kh))
params   = [curvelen,kh,h]; % [curvelength, Helmholtz param,h]

% Create a reference problem
A   = LOCAL_get_A_single_diag(C,1:ntot,params);
qqs = randn(1,size(xx_int,2));
u   = LOCAL_pot_from_charges(C([1,4],:),xx_int,qqs,kh);
q   = A\u;

% Check that the fuction for evaluating off-diagonal blocks works
ind1a = 24:26;
ind2a = [33:36,75];
ind1b = 1:10;
ind2b = 11:ntot;
fprintf(1,'The errors in evaluating the off-diagonal blocks are %10.2e and %10.3e\n',...
        norm(A(ind1a,ind2a) - LOCAL_get_A_single_offd(C,ind1a,ind2a,params)),...
        norm(A(ind1b,ind2b) - LOCAL_get_A_single_offd(C,ind1b,ind2b,params)))

% Check the error between the potential generated by the integral equation
% and the reference solution.
v = LOCAL_evalpot(xx_ext, C, q, params);
w = LOCAL_pot_from_charges(xx_ext,xx_int,qqs,kh);
fprintf(1,'The maximum of the error = %10.3e\n', max(abs(v - w)))
fprintf(1,'The mean square error    = %10.3e\n', sqrt(sum((abs(v - w).^2))/length(v)))

% Compute the integral of the boundary potential
ww = h*sqrt(C(2,:).^2 + C(5,:).^2);
fprintf(1,'The integral of q        = %10.3e\n', ww*q)

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test function for scattering from two contours.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_test_two_cont()

fprintf(1,'\n=== Running a test on two contours.\n\n')

n1   = 200;  % The number of points on each contour.
n2   = n1;

%%% Set the geometry.
[C1,len1,xx_int1,~] = LOCAL_get_geometry_star(n1,[0;0]);
[C2,len2,xx_int2,~] = LOCAL_get_geometry_star(n2,[4;0]);
h      = len1/n1;
C      = [C1,C2];
x1min  = min(C(1,:));
x1max  = max(C(1,:));
x2min  = min(C(4,:));
x2max  = max(C(4,:));
xx_ext = [(x1min-1)*ones(1,5),x1min+(x1max-x1min)*sort(rand(1,5)),(x1max+1)*ones(1,5);...
          x2min+(x2max-x2min)*sort(rand(1,5)),(x2max+1)*ones(1,5),linspace(x2min,x2max,5)];
ind1   = 1:n1;
ind2   = n1 + (1:n2);

%%% Plot the geometry
plot(C(1,:), C(4,:),'k.',...
     xx_int1(1,:), xx_int1(2,:),'r.',...
     xx_int2(1,:), xx_int2(2,:),'g.',...
     xx_ext(1,:),  xx_ext(2,:),'b.',...
     'MarkerSize',30)
legend('Contour C','xxint1','xxint2','xxext')
axis equal

% Set the Helmholtz parameter.
kh       = 2;
lencont  = 2*max(sqrt(C(1,ind1).^2 + C(4,ind1).^2));  % The diameter of the contour.
fprintf(1,'The diameter of the contour in wavelengths = %7.2f\n',kh*lencont/(2*pi))
fprintf(1,'Number of points per wave-length           = %7.2f\n',n1*2*pi/(len1*kh))
params   = [NaN,kh,h]; % [curvelength, Helmholtz param,h]

% Create the matrix A:
% Note that the off-diagonal blocks are computed via special function
% that does not include the near-diagonal corrections to account for
% the singular kernel.
A11 = LOCAL_get_A_single_diag(C1,1:n1,params);
A22 = LOCAL_get_A_single_diag(C2,1:n2,params);
A12 = LOCAL_get_A_offd_noquad(C,ind1,ind2,params);
A21 = LOCAL_get_A_offd_noquad(C,ind2,ind1,params);
A   = [A11,A12;...
       A21,A22];
   
%%% Solve a potential problem.
xxs = [xx_int1,xx_int2];
qqs = randn(1,size(xxs,2));
u   = LOCAL_pot_from_charges(C([1,4],:),xxs,qqs,kh);
q   = A\u;

% Check the error between the potential generated by the integral equation
% and the reference solution.
v = LOCAL_evalpot(xx_ext, C, q, params);
w = LOCAL_pot_from_charges(xx_ext,xxs,qqs,kh);
fprintf(1,'The maximum of the error = %10.3e\n', max(abs(v - w)))
fprintf(1,'The mean square error    = %10.3e\n', sqrt(sum((abs(v - w).^2))/length(v)))

% Compute the integral of the boundary charge distribution.
ww    = h*sqrt(C(2,:).^2 + C(5,:).^2);
q_int = ww*q;
fprintf(1,'The integral of q        = %10.3e\n', q_int)

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes diagonal blocks of the coefficient matrix.
% Note that this should be used ONLY for an index set ind that is
% entirely contained within one contour.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_get_A_single_diag(C,ind,params)

kh       = params(2);
ntot     = size(C,2);
h        = params(3);
nloc     = length(ind);

[Y_g1,   X_g1   ] = meshgrid(C(1,ind), C(1,ind));
[Y_g2,   X_g2   ] = meshgrid(C(4,ind), C(4,ind));
[Y_dg1,  X_dg1  ] = meshgrid(C(2,ind), C(2,ind));
[Y_dg2,  X_dg2  ] = meshgrid(C(5,ind), C(5,ind));

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) + eye(nloc);
A   = ((nn1.*(Y_g1 - X_g1) + nn2.*(Y_g2 - X_g2)).*(1./dd).*(-kh*besselh(1, kh*dd)) + ...
       ima*kh*besselh(0, kh*dd)).*h.*sqrt(Y_dg1.^2 + Y_dg2.^2);
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];
D = abs(ind'*ones(1,length(ind)) - ones(length(ind),1)*ind);
D = abs(mod(D + 6, ntot) - 6);
[ii1, ii2] = find(D.*(D<7));
for j = 1:length(ii1)
   i1 = ii1(j);
   i2 = ii2(j);
   d  = D(i1,i2);
   A(i1,i2) = (1 + MU6(d))*A(i1,i2);
end
for j = 1:nloc
  A(j,j) = - 2*ima;
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the field caused by point charges.
%   xx_trg   = locations of target points.
%   xx_src   = locations of sources points.
%   qq       = the strength of the point charges.
%   kh       = helmholtz parameter.
%
% The output u is the vector
%
%   u(i) = sum_{j} qq(j) * H_0(kh*|xx_trg(:,i) - xx_src(:,j)|)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u = LOCAL_pot_from_charges(xx_trg,xx_src,qq,kh)

ntrg  = size(xx_trg,2);
nsrc  = size(xx_src,2);
dist1 = (xx_trg(1,:)' * ones(1,nsrc)) - ones(ntrg,1)*xx_src(1,:);
dist2 = (xx_trg(2,:)' * ones(1,nsrc)) - ones(ntrg,1)*xx_src(2,:);

dist  = sqrt(dist1.*dist1 + dist2.*dist2);

u = besselh(0,kh*dist)*reshape(qq,numel(qq),1);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function generates a "smooth star" centered at xxc.
% The contour is parameterized as follows:
%
%   C(1,i) =                      x1 coordinate of node i
%   C(2,i) =        derivative of x1 coordinate of node i
%   C(3,i) = second derivative of x1 coordinate of node i
%   C(4,i) =                      x2 coordinate of node i
%   C(5,i) =        derivative of x2 coordinate of node i
%   C(6,i) = second derivative of x2 coordinate of node i
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [C,curvelen,xx_int,xx_ext] = LOCAL_get_geometry_star(ntot,xxc)

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));

%%% Construct interior points.
ttint  = 2*pi*sort(rand(1,3));
xx_int = 0.5*rmin*[cos(ttint); sin(ttint)];

%%% Construct exterior points.
ttext  = 2*pi*sort(rand(1,5));
xx_ext = 1.5*rmax*[cos(ttext); sin(ttext)];

%%% Shift to center at xxc.
C([1,4],:) = C([1,4],:) + xxc*ones(1,size(C,2));
xx_int     = xx_int     + xxc*ones(1,size(xx_int,2));
xx_ext     = xx_ext     + xxc*ones(1,size(xx_ext,2));

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute an off-diagonal block of the matrix discretizing a single contour.
% Note that this INCLUDES corrections near the diagonal for handling the
% singularity in the kernel.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_get_A_single_offd(C,ind1,ind2,params)

ntot     = size(C,2);
kh       = params(2);
h        = params(3);

[Y_g1,  X_g1  ] = meshgrid(C(1,ind2), C(1,ind1));
[Y_g2,  X_g2  ] = meshgrid(C(4,ind2), C(4,ind1));
[Y_dg1, X_dg1 ] = meshgrid(C(2,ind2), C(2,ind1));
[Y_dg2, X_dg2 ] = meshgrid(C(5,ind2), C(5,ind1));

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*kh*besselh(0, kh*dd)).*h.*sqrt(Y_dg1.^2 + Y_dg2.^2);
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];
D = abs(ind1'*ones(1,length(ind2)) - ones(length(ind1),1)*ind2);
D = abs(mod(D + 6, ntot) - 6);
[ii1, ii2] = find(D.*(D<7));
for j = 1:length(ii1)
  i1 = ii1(j);
  i2 = ii2(j);
  d  = D(i1,i2);
  A(i1,i2) = (1 + MU6(d))*A(i1,i2);
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes off-diagonal blocks of the coefficient matrix,
% but without adding any quadrature corrections. This is intended to be 
% used for blocks for which ind1 and ind2 belong to different contours.
% (It can also be used for well-separated blocks regardless of location.)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_get_A_offd_noquad(C,ind1,ind2,params)

kh = params(2);
h  = params(3);

[Y_g1,  X_g1  ] = meshgrid(C(1,ind2), C(1,ind1));
[Y_g2,  X_g2  ] = meshgrid(C(4,ind2), C(4,ind1));
[Y_dg1, X_dg1 ] = meshgrid(C(2,ind2), C(2,ind1));
[Y_dg2, X_dg2 ] = meshgrid(C(5,ind2), C(5,ind1));

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*kh*besselh(0, kh*dd)).*h.*sqrt(Y_dg1.^2 + Y_dg2.^2);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function evaluates a layer potential with the combined field kernel.
% xx     is the matrix of target points, 
% C      is the contour.
% sigma  is the source distribution on C.
% params contains the Helmholtz parameter, the grid spacing "h", etc.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function vv = LOCAL_evalpot(xx, C, sigma, params)

kh = params(2);
h  = params(3);

n = size(C,2);
m = size(xx,2);

X_g1  = xx(1,:)'*ones(1,n);
X_g2  = xx(2,:)'*ones(1,n);
Y_g1  = ones(m,1)*C(1,:);
Y_g2  = ones(m,1)*C(4,:);
Y_dg1 = ones(m,1)*C(2,:);
Y_dg2 = ones(m,1)*C(5,:);

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);
EVAL = ((nn1.*(Y_g1 - X_g1) + nn2.*(Y_g2 - X_g2)).*(1./dd).*(-kh*besselh(1, kh*dd)) + ...
        ima*kh*besselh(0, kh*dd)).*h.*sqrt(Y_dg1.^2 + Y_dg2.^2);
   
vv = EVAL*sigma;

return
