% scattering from a perfectly conducting disk of radius one, or equally scattering from an
% acoustically hard/soft disk
% k0 - the wavenumber

% x,y - coordinates where the numerical solution is known
% terms - number of positive/negative terms in the series

x = curve0(1,:);
y = curve0(2,:);

terms = 10;

[theta,r0] = cart2pol(x,y); 

%-------acoustically hard object (normal derivative of the total field is zero on the boundary)------- 
uHard=zeros(size(x));
for n=-terms:1:terms
    uexT=-1i^(n)*(k0*besselj(n-1,k0)-n*besselj(n,k0)).*(besselj(n,k0*r0)+1i*bessely(n,k0*r0)).*exp(1i*n.*theta);
    uexN=k0*(besselj(n-1,k0)+1i*bessely(n-1,k0))-n*(besselj(n,k0)+1i*bessely(n,k0));
    uHard =uHard+uexT./uexN;
end
% Note: uHard is the scattered field

%---------acoustically soft object (the total field is zero on the boundary)---------------- 
uSoft=zeros(size(x));
for n=-terms:1:terms
    uSoft =uSoft+(-1i^(n)*besselj(n,k0).*(besselj(n,k0*r0)+1i*bessely(n,k0*r0)).*exp(1i*n.*theta))./(besselj(n,k0)+1i*bessely(n,k0));
end

% Note: uSoft is the scattered field

%Please send an email to christian.engstroem@sam.math.ethz.ch if you find
%errors in the code or have problems using it.