1% Example 9: Helmholtz sound-soft scattering by a peanut shape, using the
2% combined-field (CFIE) representation D - i*k*S.
4mip load --install magland/magland/chunkie;
6tic;
8% chunker discretization of a peanut-shaped domain
9modes = [1.25,-0.25,0,0.5];
10ctr = [0;0];
11chnkr = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr));
13% incident plane wave defines the boundary condition
14kwav = 3*[-1,-5];
15pwfun = @(r) exp(1i*kwav*r(:,:));
16rhs = -pwfun(chnkr.r); rhs = rhs(:);
18% CFIE kernel D - i*k*S
19zk = norm(kwav);
20coefs = [1,-1i*zk];
21kerncfie = kernel('h','c',zk,coefs);
23% matrix discretization of the BIE (with 0.5*I from exterior jump)
24sysmat = chunkermat(chnkr,kerncfie);
25sysmat = sysmat + 0.5*eye(chnkr.npt);
27% solve
28sigma = gmres(sysmat,rhs,[],1e-10,100);
30% grid for plotting (exterior)
31x1 = linspace(-5,5,300);
32[xx,yy] = meshgrid(x1,x1);
33targs = [xx(:).'; yy(:).'];
34in = chunkerinterior(chnkr,{x1,x1});
35uu = nan(size(xx));
37% evaluate scattered field outside, then add incident field
38uu(~in) = chunkerkerneval(chnkr,kerncfie,sigma,targs(:,~in));
39uu(~in) = uu(~in) + pwfun(targs(:,~in)).';
41toc;
43figure(1); clf
44h = pcolor(xx,yy,real(uu)); set(h,'EdgeColor','none'); colorbar
45colormap(redblue); umax = max(abs(uu(:))); clim([-umax,umax]);
46hold on
47plot(chnkr,'k')