% Example 9: Helmholtz sound-soft scattering by a peanut shape, using the
% combined-field (CFIE) representation D - i*k*S.

mip load --install magland/magland/chunkie;

tic;

% chunker discretization of a peanut-shaped domain
modes = [1.25,-0.25,0,0.5];
ctr = [0;0];
chnkr = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr));

% incident plane wave defines the boundary condition
kwav = 3*[-1,-5];
pwfun = @(r) exp(1i*kwav*r(:,:));
rhs = -pwfun(chnkr.r); rhs = rhs(:);

% CFIE kernel D - i*k*S
zk = norm(kwav);
coefs = [1,-1i*zk];
kerncfie = kernel('h','c',zk,coefs);

% matrix discretization of the BIE (with 0.5*I from exterior jump)
sysmat = chunkermat(chnkr,kerncfie);
sysmat = sysmat + 0.5*eye(chnkr.npt);

% solve
sigma = gmres(sysmat,rhs,[],1e-10,100);

% grid for plotting (exterior)
x1 = linspace(-5,5,300);
[xx,yy] = meshgrid(x1,x1);
targs = [xx(:).'; yy(:).'];
in = chunkerinterior(chnkr,{x1,x1});
uu = nan(size(xx));

% evaluate scattered field outside, then add incident field
uu(~in) = chunkerkerneval(chnkr,kerncfie,sigma,targs(:,~in));
uu(~in) = uu(~in) + pwfun(targs(:,~in)).';

toc;

figure(1); clf
h = pcolor(xx,yy,real(uu)); set(h,'EdgeColor','none'); colorbar
colormap(redblue); umax = max(abs(uu(:))); clim([-umax,umax]);
hold on
plot(chnkr,'k')