% kernel function: fills kernel for laplace or helm eqn
%               also gives nodes and weights
% inputs k : wavenumber
%        m : number of nodes
%
% outputs ker: kernel matrix
%         [z dzdx, w, n] : parameterization


function [ker, z, dzdx, w, n]=hlker2(m,k)


%first, set up nodes and weights

[z,dzdx,n,kap]=param(m); %parameterizes equation

w=zeros(1,numel(dzdx))+2*pi/m; %periodic trapezoid weights


%now to set up kernel matrix
ker=zeros(m,m); %initalizes kernel matrix

%first, to fill non-diagonal entries:
if k==0 %laplace-type equation

%function for evaluating dphi/dny
dphi=@(x,y,n) -1/(2*pi)*(n'*(y-x))/((y-x)'*(y-x));

for i=1:m
    for j=1:m
       if ~(i==j)
        ker(i,j)= dphi(z(:,i),z(:,j),n(:,j))*dzdx(j)*w(j);
       end
    end
end

    
else %helmholtz-type equation

[xa xb]=meshgrid(z(1,:),z(1,:)); %sets up x and y grids
[ya yb]=meshgrid(z(2,:),z(2,:));
n1=(zeros(size(n(1,:)))'+1)*n(1,:);
n2=(zeros(size(n(2,:)))'+1)*n(2,:);

    
dphi=@(x1,y1,x2,y2,n1,n2,k) 1i*k/4*(n1'.*(x1-x2)+n2'.*(y1-y2))./...
    sqrt((x1-x2).^2+(y1-y2).^2).*besselh(1,k*...
    (sqrt((x1-x2).^2+(y1-y2).^2))); %vectorized kernel function

ker=(dphi(xa,ya,xb,yb,n1,n2,k).*((dzdx.*w)'*(zeros(size(dzdx))+1))).';

end
   
for i=1:m; %does diagonal entries- same for laplace or helm eqns
    ker(i,i)=-kap(i)/(4*pi)*dzdx(i)*w(i);
end
