%% clean up
clear
close all

%% load data
fid=fopen('asteroid.grid','r');
num=fscanf(fid,'%i',2);
np=num(1);
nt=num(2);
ptmp=fscanf(fid,'%f',[2,np]);
ttmp=fscanf(fid,'%f',[3,nt]);
fclose(fid);
% Change for Matlab
nodes=ptmp';
triangles=ttmp'+1;
clear ptmp ttmp
% Load Neumann boundary data
bdryData = load('normal_derivative.dat'); % here, the node numbering does not start from 0!
ns=size(bdryData,1);
% Plot the mesh
triplot(triangles,nodes(:,1),nodes(:,2),'r');
hold on;
x1 = nodes(bdryData(:,1),:);
x2 = nodes(bdryData(:,2),:);
% for i=1:ns
%     plot([x1(i,1),x2(i,1)],[x1(i,2),x2(i,2)]);
% end
vec = x2 - x1;
quiver(x1(:,1),x1(:,2),vec(:,1),vec(:,2));

y = [15,-5];
y = [1000,0];

%% computation of u on the boundary
% assemble V and K
V = zeros(ns,ns);
K = V;
for i = 1:ns
    for j = 1:ns
        yi = .5*(nodes(bdryData(i,1),:)+nodes(bdryData(i,2),:));
        xa = nodes(bdryData(j,1),:);
        xb = nodes(bdryData(j,2),:);
        lVec = xa - xb;
        aVec = xa - yi;
        l = norm(lVec,2);
        a = norm(aVec,2);
        if ( i == j )
            V(i,i) = l * ( 1 + log( 2 / l ) );
        else
            alpha = acos( lVec * aVec' / ( l * a ) );
            if ( ( alpha > 1.e-5 ) && ( alpha < ( pi - 1.e-5 ) ) )
               V(i,j) = -(-l+a*sin(alpha)*atan((l-a*cos(alpha))/(a*sin(alpha)))+.5*(l-a*cos(alpha))*log(l^2+a^2-2*a*l*cos(alpha))...
                        -a*sin(alpha)*atan(-cot(alpha))+.5*a*cos(alpha)*log(a^2));
            elseif ( alpha <= 1.e-5 ) 
                V(i,j) = -(-l+.5*(l-a)*log(l^2+a^2-2*a*l)+a*log(a));
            else
                V(i,j) = -(-l+.5*(l+a)*log(l^2+a^2+2*a*l)-a*log(a));
            end                
            nu = -[lVec(2),-lVec(1)];
            A = aVec * nu';
            C = a^2;
            D = -aVec * lVec';
            E = l^2;
            sqrtDiscr = sqrt(C*E-D^2);
            if ( sqrtDiscr < 1.e-5 )
                K(i,j) = -A/l*(1/(a+l)-1/a);
            else
                K(i,j) = -A/sqrtDiscr*(atan((E+D)/sqrtDiscr)-atan(D/sqrtDiscr));
            end
        end
    end
end

%u = (1/(2*pi)*K+.5*eye(ns))\(1/(2*pi)*V*bdryData(:,3));
op = (-1/(2*pi)*K+.5*eye(ns));
rhs = -(1/(2*pi)*V*bdryData(:,3));
u = op\rhs

%% computation of u outside
% assemble V and K
Vvec = zeros(1,ns);
Kvec = Vvec;
for j = 1:ns
    xa = nodes(bdryData(j,1),:);
    xb = nodes(bdryData(j,2),:);
    lVec = xa - xb;
    aVec = xa - y;
    l = norm(lVec,2);
    a = norm(aVec,2);
    alpha = acos( lVec * aVec' / ( l * a ) );
    if ( ( alpha > 1.e-5 ) && ( alpha < ( pi - 1.e-5 ) ) )
       Vvec(1,j) = -(-l+a*sin(alpha)*atan((l-a*cos(alpha))/(a*sin(alpha)))+.5*(l-a*cos(alpha))*log(l^2+a^2-2*a*l*cos(alpha))...
                -a*sin(alpha)*atan(-cot(alpha))+.5*a*cos(alpha)*log(a^2));
    elseif ( alpha <= 1.e-5 ) 
        Vvec(1,j) = -(-l+.5*(l-a)*log(l^2+a^2-2*a*l)+a*log(a));
    else
        Vvec(1,j) = -(-l+.5*(l+a)*log(l^2+a^2+2*a*l)-a*log(a));
    end                
    nu = -[lVec(2),-lVec(1)];
    A = aVec * nu';
    C = a^2;
    D = -aVec * lVec';
    E = l^2;
    sqrtDiscr = sqrt(C*E-D^2);
    if ( sqrtDiscr < 1.e-5 )
        Kvec(1,j) = -A/l*(1/(a+l)-1/a);
    else
        Kvec(1,j) = -A/sqrtDiscr*(atan((E+D)/sqrtDiscr)-atan(D/sqrtDiscr));
    end
end

u_y = ( -Vvec*bdryData(:,3)+Kvec*u )/(2*pi)

