clear
k = 2*pi; % so that the wavelength is 1
R = 0.5 % so that the radius of the circle is 0.5 wavelengths
x = [2*R,2*R] % picking a position outside the circle where u should be 0
x0 = [R/4,0] % putting the source off-centre
N = [8 16 32 64 128 256];
for m = 1:length(N)
    [xx,A] = circ(N(m),R); % dividing the circle into N boundary elements
    h(m) = A; % the length of each boundary element
    v = bem_solve(xx,A*ones(N(m),1),x0,k,N(m)); % solving the integral equation by BEM
    u_bem(m) = pressure(xx,A*ones(N(m),1),x0,k,N(m),v,x); % computing pressure at x by BEM
end
plot(h,abs(u_bem))
xlabel('Element length as a fraction of the wavelength')
ylabel('Relative error')
title('Predicted BEM pressure outside the circle')
grid
disp(' ')
disp('      h/lambda       |u|')
disp([h',abs(u_bem)'])