%Brad Nelson 2/12/12
%Math 126 HW6 #1 
%Helmholtz BVP!

%m=100; %number of nodes
counter=0;
error=zeros(size(10:10:200));
for m=10:10:200
counter=counter+1;

k=6; %wavenumber

[ker, z, dzdx, w, n]=hlker(m,k);

% %for visualizing kernel
% figure;
% imagesc(real(ker))
% colorbar
% print -depsc 1a.eps;

% now to solve linear system...
A=eye(m); %initializes A matrix
A=A-2*ker; %A matrix for system

soln=@(x) 5i*besselh(0,k*norm(x-[0;2])); %actual solution for system

b=zeros(size(z(1,:)))';
for j=1:numel(z(1,:))
    b(j)=-2*soln(z(:,j)); %solution vector
end

tau=A\b; %solves linear system for tau vector

x=[0.2; 0.1]; %location to evaluate solution potential at

fx=0; % holds answer at x

dphi=@(x,y,n,k) 1i*k/4*n'*(x-y)/norm(x-y)*besselh(1,k*(norm(x-y)));

%nystrom method of reconstructing answer at x
for i=1:m %riemann sum
    fx=fx+dphi(x,z(:,i),n(:,i),k)*tau(i)*dzdx(i)*w(i);
end

error(counter)=abs(fx-soln(x));

end

figure;
loglog(10:10:200,error);
xlabel('number of nodes');
ylabel('error');
%print -depsc 1c.eps;
