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



m=50; %number of nodes

k=6; %wavenumber

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

soln=@(x,y) 5i*besselh(0,k*sqrt(x.^2+(y-2).^2));%actual solution

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

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

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

x1=-1.3:0.02:1.3; x2=-x1;
[xx yy]=meshgrid(x1,x2); %sets up grid to evaluate soln on
num=numel(x1); %number of data points on grid axis

fx=zeros(size(xx)); %solution grid

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

act=soln(xx,yy); %actual solution

errors=log(abs(real(act-fx)));

%plotting time!
figure;
contourf(x1,x2,real(errors));
%title('log error')
colorbar;
print -depsc 1b.eps;

figure;
imagesc(x1,x2,real(act));
set(gca,'YDir','normal');
%title('act')
colorbar;
print -depsc 1b2.eps;

%sets outside of boundary to 0
rad=0;
bdry=0;
num=numel(x1); %number of data points on grid axis
for h=1:num
for j=1:num
    rad=norm([xx(j,h),yy(j,h)]); %radius
    bdry=1+0.3*cos(3*acos(xx(j,h)/rad)); 
    if rad>bdry %checks to see if point is inside boundary
        fx(j,h)=0;
    end
end
end

figure;
imagesc(x1,x2,real(fx),[-2 2]);
set(gca,'YDir','normal');
%title('comp')
colorbar;
print -depsc 1b3.eps;

