1mip load --install magland/magland/chunkie;
2tic;
4% planewave definitions
6kvec = 20*[1;-1.5];
7zk = norm(kvec);
8planewave = @(kvec,r) exp(1i*sum(bsxfun(@times,kvec(:),r(:,:)))).';
10% discretize domain
12narms = 5;
13amp = 0.5;
14chnkr = chunkerfunc(@(t) starfish(t,narms,amp),struct('maxchunklen',4/zk));
16% build CFIE and solve
18fkern = kernel('helm','c',zk,[1,-zk*1i]);
19sysmat = chunkermat(chnkr,fkern);
20sysmat = 0.5*eye(chnkr.k*chnkr.nch) + sysmat;
22rhs = -planewave(kvec(:),chnkr.r(:,:));
23sol = gmres(sysmat,rhs,[],1e-13,100);
25% evaluate at targets
27x1 = linspace(-3,3,400);
28[xxtarg,yytarg] = meshgrid(x1,x1);
29targets = [xxtarg(:).';yytarg(:).'];
31in = chunkerinterior(chnkr,targets);
32out = ~in;
34uscat = chunkerkerneval(chnkr,fkern,sol,targets(:,out));
36uin = planewave(kvec,targets(:,out));
37utot = uscat(:)+planewave(kvec,targets(:,out));
39% plot
41maxu = max(abs(utot(:)));
42figure()
43zztarg = nan(size(xxtarg));
44zztarg(out) = utot;
45h=pcolor(xxtarg,yytarg,imag(zztarg));
46set(h,'EdgeColor','none')
47hold on
48plot(chnkr,'LineWidth',2)
49axis equal tight
50colormap(redblue)
51caxis([-maxu,maxu])
53toc;