%Brad Nelson 2/12/12
%Math 126 HW6 #2b
%Helmholtz scattering problem!
%Error at (-2,0)

ms=60:20:280; %ms to check values at
u=zeros(size(ms));
counter=0; %to be used for indexing each loop
for m=ms; %number of nodes
counter=counter+1; %keeps track of loop number
    
k=10; %wavenumber

[ker, z, dzdx, w, n]=hlker(m,k); %makes kernel, parameterizes curve

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

inc=@(x,y) exp(1i*10*(x*cos(0.2)+y*sin(0.2))); %incident wave

b=-2*inc(z(1,:),z(2,:)).'; %solution vector

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

%dphi=@(x,y,n,k) 1i*k/4*n'*(x-y)/norm(x-y)*besselh(1,k*(norm(x-y))); %kernel
dphi=@(x,y,z,n,k) 1i*k/4*(n(1)*(x-z(1))+n(2)*(y-z(2)))./...
    sqrt((x-z(1)).^2+(y-z(2)).^2).*besselh(1,k*...
    (sqrt((x-z(1)).^2+(y-z(2)).^2))); %vectorized kernel function

x=[-2;0]; %point to check convergence at
ui=inc(x(1),x(2)); %incident wave at x

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

u(counter)=ui+us; %total wave

end

figure;
plot(ms,real(u),'.-')
xlabel('N'); ylabel('u');
print -depsc 2b1.eps;

uact=u(numel(u));
figure;
semilogy(ms,abs(real(u-uact)));
xlabel('N'); ylabel('error');
print -depsc 2b2.eps;
