%Brad Nelson 2/12/12
%Math 126 HW6 #2 
%Helmholtz scattering problem!


m=100; %number of nodes

k=10; %wavenumber

[ker, z, dzdx, w, n]=hlker2(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


fineness=250; %number of points on each grid axis
x1=linspace(-4,4,fineness); x2=-x1;
[xx yy]=meshgrid(x1,x2); %sets up grid to evaluate soln on

ui=inc(xx,yy); %incident wave

us=zeros(size(xx)); %scattered wave

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

%sets inside 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
        ui(j,h)=0;
        us(j,h)=0;
    end
end
end

u=ui+us; %total wave

%now to plot pretty pictures!
figure;
imagesc(x1,x2,real(ui));
set(gca,'YDir','normal');
hold on;
plot(z(1,:),z(2,:),'LineWidth',2.5);
axis square;
%title('incident wave u^i');
colorbar
hold off;
%print -depsc 2a1.eps;

figure;
imagesc(x1,x2,real(us),[-1.2 1.2]);
set(gca,'YDir','normal')
hold on;
plot(z(1,:),z(2,:),'LineWidth',2.5);
axis square;
%title('scattered wave u^s');
colorbar
hold off;
%print -depsc 2a2.eps;

figure;
imagesc(x1,x2,real(u),[-2.2 2.2]);
set(gca,'YDir','normal')
hold on;
plot(z(1,:),z(2,:),'LineWidth',2.5);
axis square;
%title('u=u^i+u^s');
colorbar
hold off;
%print -depsc 2a3.eps;


