clear
k = 2*pi; % so that the wavelength is 1
lambda = 2*pi/k; % the wavelength
R = 0.5 % so that the radius of the circle is 0.5 wavelengths
x = [R/2,R/2] % picking a position at which we will compute the pressure
x0 = [0,0] % putting the source at the centre of the circle
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 = pressure_grad(xx,A*ones(N(m),1),x0,k,N(m),v,x) % computing pressure gradient at x by BEM
    u = pressure_exact_grad(R,k,x) % the exact pressure gradient
    Rel_error(m) = norm(u-u_bem)/norm(u); % the relative error in the BEM pressure gradient value at x
end
plot(h/lambda,Rel_error)
xlabel('Element length as a fraction of the wavelength')
ylabel('Relative error')
title('Error in BEM calculation of pressure gradient')
grid
disp(' ')
disp('    h/lambda  Relative error')
disp([h',Rel_error'])