% This function sets up and solves an exterior Neumann problem on 
% a domain exterior to a number of contours in the plane.
%
% A single layer Ansatz is used.

function main_hw07()
randn('seed',0)
rand( 'seed',0)

nstar = 100;
nside = 4;

flag_geom    = 'Nell';
[C,curvelen,xx_int,xx_ext] = LOCAL_construct_contour(nstar,nside,flag_geom);
ntot     = size(C,2);
flag_pot = [];
params   = [];

figure(1)
plot(C(1,[1:ntot,1]), C(4,[1:ntot,1]    ), 'k.',...
     xx_int(1,:), xx_int(2,:),'bx',...
     xx_ext(1,:), xx_ext(2,:),'rx')
axis equal
legend('The contour','xx-int','xx-ext')
fprintf(1,'ntot = %d\n',ntot)

%%% Test the single layer potential calculation for an exterior Neumann problem.
fprintf(1,'\n=== Solve an "indirect" BIE based on the single layer Ansatz (exterior Neumann problem).\n')
qq              = randn(size(xx_int,2),1);
[u_exact,~,~]   = LOCAL_pot_from_charges(xx_ext,xx_int,qq);
[~,ud1_C,ud2_C] = LOCAL_pot_from_charges(C([1,4],:),xx_int,qq);
vel             = sqrt(C(2,:).*C(2,:) + C(5,:).*C(5,:));
n1              =  C(5,:)./vel;
n2              = -C(2,:)./vel;
un_C            = n1'.*ud1_C + n2'.*ud2_C;
Dadj            = LOCAL_get_D_adj_op(C,1:ntot);
h               = curvelen/ntot;
A_DT            = -0.5*eye(ntot) + Dadj.*(h*ones(ntot,1)*vel);
sig             = A_DT\un_C;
u_S             = LOCAL_eval_S(xx_ext,C,curvelen,sig,flag_pot,params);
fprintf(1,'    For an exterior point: |u_exact - S*sig|          = %13.5e\n',...
        norm(u_exact - u_S))
fprintf(1,'    Integral of un_C                                  = %13.5e\n',h*vel*un_C)
fprintf(1,'    Integral of sigma                                 = %13.5e\n',(h*vel*sig))
fprintf(1,'    sum(qq)                                           = %13.5e\n',sum(qq));
%ss = svd(A_DT);
%fprintf(1,'    The condition number of the system                = %13.5e\n',ss(1)/ss(end))
fprintf(1,'    The condition number of the system                = %13.5e\n',condest(A_DT))

keyboard
    
return

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

function [C,curvelen,xx_int,xx_ext] = LOCAL_construct_contour(nstar,nside,flag_geom)

%%% Set up a reference contour.
r        = 0.2;
k        = 5;
tt       = linspace(0,2*pi*(1 - 1/nstar),nstar);
D        = zeros(6,nstar );
D(1,:)   =   1.5*cos(tt) + (r/2)*            cos((k+1)*tt) + (r/2)*            cos((k-1)*tt);
D(2,:)   = - 1.5*sin(tt) - (r/2)*(k+1)*      sin((k+1)*tt) - (r/2)*(k-1)*      sin((k-1)*tt);
D(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);
D(4,:)   =       sin(tt) + (r/2)*            sin((k+1)*tt) - (r/2)*            sin((k-1)*tt);
D(5,:)   =       cos(tt) + (r/2)*(k+1)*      cos((k+1)*tt) - (r/2)*(k-1)*      cos((k-1)*tt);
D(6,:)   = -     sin(tt) - (r/2)*(k+1)*(k+1)*sin((k+1)*tt) + (r/2)*(k-1)*(k-1)*sin((k-1)*tt);
D        = 0.25*D;

%%% Make the copies.
C = zeros(6,nstar*nside*nside);
for n1 = 1:nside
  for n2 = 1:nside
    ind      = (n1-1)*nside*nstar + (n2-1)*nstar + (1:nstar);
    C(:,ind) = [n1 - 1 + D(1,:); D(2:3,:); n2 - 1 + D(4,:); D(5:6,:)];
  end
end
curvelen = nside*nside*2*pi;

n_int    = 10;
R_min    = min(sqrt(D(1,:).^2 + D(4,:).^2));
tt       = linspace(0,2*pi*(1 - 1/n_int),n_int);
xx_int   = 0.7*R_min*[cos(tt); sin(tt)];

border   = 0.3;
x1min    = min(C(1,:)) - border;
x1max    = max(C(1,:)) + border;
x2min    = min(C(4,:)) - border;
x2max    = max(C(4,:)) + border;
xx_ext   = [[x1min + (x1max - x1min)*rand(1,10); x2min*ones(1,10)],...
            [x1min + (x1max - x1min)*rand(1,10); x2max*ones(1,10)],...
            [x1min*ones(1,10); x2min + (x2max - x2min)*rand(1,10)],...
            [x1max*ones(1,10); x2min + (x2max - x2min)*rand(1,10)]];

return

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

function [u,ud1,ud2] = LOCAL_pot_from_charges(xx_trg,xx_src,qq)

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

distsq = dist1.*dist1 + dist2.*dist2;
u      = -(1/(4*pi))*log(distsq)*reshape(qq,nsrc,1);
ud1    = (-dist1./(2*pi*distsq))*reshape(qq,nsrc,1);
ud2    = (-dist2./(2*pi*distsq))*reshape(qq,nsrc,1);

return

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

function u = LOCAL_eval_S(xx, C, curvelen, sig, flag_pot, params)

%kh = params(2);

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

vel = sqrt(C(2,:).*C(2,:) + C(5,:).*C(5,:));

X1  = xx(1,:)'*ones(1,n);
X2  = xx(2,:)'*ones(1,n);
Y1  = ones(m,1)*C(1,:);
Y2  = ones(m,1)*C(4,:);
DD  = (X1 - Y1).*(X1 - Y1) + (X2 - Y2).*(X2 - Y2);

A   = -(1/(4*pi))*log(DD);

u   = A*(h*(vel').*sig);

return

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

function D = LOCAL_get_D_adj_op(C,ind)

nloc = length(ind);

vel = sqrt(C(2,ind).*C(2,ind) + C(5,ind).*C(5,ind));
n1  =  C(5,:)./vel;
n2  = -C(2,:)./vel;
N1  = (n1')*ones(1,nloc);
N2  = (n2')*ones(1,nloc);
D1  = C(1,ind)' * ones(1,nloc) - ones(nloc,1) * C(1,ind);
D2  = C(4,ind)' * ones(1,nloc) - ones(nloc,1) * C(4,ind);
DSQ = D1.*D1 + D2.*D2 + eye(nloc);

D   = -(N1.*D1 + N2.*D2)./(2*pi*DSQ);

for icount = 1:nloc
  i                = ind(icount);
  D(icount,icount) = (-C(6,i)*C(2,i) + C(3,i)*C(5,i))/...
                      (4*pi*vel(i)^3);
end

return

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

function S = LOCAL_get_S_op(C,ind)

nloc = length(ind);
D1   = C(1,ind)' * ones(1,nloc) - ones(nloc,1) * C(1,ind);
D2   = C(4,ind)' * ones(1,nloc) - ones(nloc,1) * C(4,ind);
DSQ  = D1.*D1 + D2.*D2 + eye(nloc);
S    = (-1/(4*pi))*log(DSQ);

return