%This discretizes a double layer potential on the boundary of the ellipse
%x^2+(y/aa)^2=1. 
%This is output in the matrix Dm
%How many samples?
clear
N=1000;
aa=2;
%Here are the samples
ts=2*pi*linspace(0,1,N+1);
dmat=zeros(N+1);
%Construct source sample points along the ellipse
x2=cos(ts(:));
y2=aa*sin(ts(:));
figure, plot(x2,y2)
%Construct the length element for this parametrization
ds=sqrt(sin(ts(1:N)).^2+aa^2*cos(ts(1:N)).^2)*2*pi/N;
%Build the matrix for the compact part of the double layer:
for j=1:(N+1)
    for k=1:(N+1)
        %Check to see if we're on the diagonal, if not, then
        if ~(k == j) 
            dmat(j,k)=aa*(1-cos(ts(j)-ts(k)));
            dmat(j,k)=dmat(j,k)/((cos(ts(j))-cos(ts(k)))^2+...
                aa^2*(sin(ts(j))-sin(ts(k)))^2);
            else
            %Compute the diagonal accurately
            dmat(j,k)=aa/(2*(sin(ts(j))*sin(ts(j))+...
                aa^2*cos(ts(j))*cos(ts(j))));
        end
         %include the norm of the normal vector and length element
         dmat(j,k)=dmat(j,k)*sqrt((aa^2*cos(ts(k))*cos(ts(k))+...
             sin(ts(k))*sin(ts(k)))/(aa^2*cos(ts(j))*cos(ts(j))+...
             sin(ts(j))*sin(ts(j))));
    end
end

    %Normalize for integration purposes:
    %Here we use the rectangular rule:
    Dm=-eye(N)/2+dmat(1:N,1:N)/N;
    %This shows how close 0 is to being an eigenvalue of this matrix
    [u,dd,v]=svd(Dm);
   % eigerr=(Dm*(1./ds'));
   eigerr=(Dm*v(:,N));
    norm(eigerr)
    [u,dd,v]=svd(Dm);
    dd=diag(dd);
    %Let's check the largest and smallest singular values
   [dd(1:5) dd(N-4:N)]
    %We compute the condition number after adding the "ones" matrix
    Dm1=Dm+dd(N-1)*u(:,N)*v(:,N)';
    %cond(Dm)
    cond(Dm1)
    %We define the outer unit normal vector
    nx=aa*cos(ts(1:N));
    ny=sin(ts(1:N));
    nn=sqrt(nx.*nx+ny.*ny);
    nrm_vec(:,1)=nx./nn;
    nrm_vec(:,2)=ny./nn;
    %Construct target sample points for the simple layer potential:
    Nr=50;
    Nth=100;
    rr=linspace(0,1,Nr+1);
    ths=linspace(0,2*pi,Nth);
    rr=rr(1:Nr);
    %Compute the (x,y) coordinates
    for j=1:Nr
        for k=1:Nth
            x(j,k)=rr(j)*cos(ths(k));
            y(j,k)=rr(j)*aa*sin(ths(k));
        end
    end
    %Construct the single layer potential:
    for j=1:Nr
        for k=1:Nth
            sl(j,k,1:N)=(log(sqrt((x(j,k)-x2(1:N)).^2+...
                (y(j,k)-y2(1:N)).^2)).*ds')/(2*pi);
        end
    end
    
    
    
    

