You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: delaunay_neighbors.m
    path of file: ./repos/gptoolbox/mesh
    
    the code of the file until where you have to start completion: function N = delaunay_neighbors(q,P)
  % DELAUNAY_NEIGHBORS Find the Delaunay neighbors of q in a list of points P.
  % 
  % N = delaunay_neighbors(q,P)
  % 
  % Experimentally for uniformly randomly distributed points the number of
  % iterations is O(|P|) and specifically seems to do 1.355*|P| "iterations".
  % There's a sort on the angles so at best this implementation is O(|P|log|P|),
  % but it's very simple... so maybe the constants and overheads are much
  % smaller than doing a full Delaunay triangulation.
  %
  % Inputs:
  %   q  2D position of query point
  %   P  #P by 2 list of nearby points
  % Outputs:
  %   N  #N list of indices into P of neighboring points
  %

  % Example:
  %   q = [0 0];
  %   P = rand(60,2)*2-1;
  %   N = delaunay_neighbors(q,P);
  %   Nm = find(adjacency_matrix(delaunay([P;q]))*sparse(size(P,1)+1,1,1,size(P,1)+1,1));
  %   assert(isempty(setxor(N,Nm)));

  
  % 2----1
  % |θ₂ / \
  % |  /   \
  % | /     \
  % |/θq   θ₀\
  % q---------0
  %% Angle at p0-q-p1
  %thq = th(p1)-th(p0);
  %% Law of sines: θq/|01| = θ₀/|q1|
  %% θ₀ = |q1|*θq/|01|
  %th0 = l(p1)*thq/sqrt(sum((P(p0,:)-P(p1,:)).^2,2));
  % Vector from p0 to p1

  % Inputs:
  %   P  #P by 2 list of nearby points
  %   l  #P by 2 list of distances to origin: l = normrow(P)
  %   p0  first 2d point
  %   p1  second 2d point
  %   p2  third 2d point
  % Outputs:
  %   flag  whether edge from origin to p1 forms a Delaunay edge between points
  %     p1 and p2.
  % 
  isdel = @(P,l,p0,p1,p2) ...
    acos((-P(p0,:)/l(p0))*(((P(p1,:)-P(p0,:))/sqrt(sum((P(p1,:)-P(p0,:)).^2,2)))')) + ...
    acos((-P(p2,:)/l(p2))*(((P(p1,:)-P(p2,:))/sqrt(sum((P(p1,:)-P(p2,:)).^2,2)))')) <= pi;

  %% Or perhaps: http://stackoverflow.com/a/8523979

    
function N = delaunay_neighbors(q,P)
  % DELAUNAY_NEIGHBORS Find the Delaunay neighbors of q in a list of points P.
  % 
  % N = delaunay_neighbors(q,P)
  % 
  % Experimentally for uniformly randomly distributed points the number of
  % iterations is O(|P|) and specifically seems to do 1.355*|P| "iterations".
  % There's a sort on the angles so at best this implementation is O(|P|log|P|),
  % but it's very simple... so maybe the constants and overheads are much
  % smaller than doing a full Delaunay triangulation.
  %
  % Inputs:
  %   q  2D position of query point
  %   P  #P by 2 list of nearby points
  % Outputs:
  %   N  #N list of indices into P of neighboring points
  %

  % Example:
  %   q = [0 0];
  %   P = rand(60,2)*2-1;
  %   N = delaunay_neighbors(q,P);
  %   Nm = find(adjacency_matrix(delaunay([P;q]))*sparse(size(P,1)+1,1,1,size(P,1)+1,1));
  %   assert(isempty(setxor(N,Nm)));

  
  % 2----1
  % |θ₂ / \
  % |  /   \
  % | /     \
  % |/θq   θ₀\
  % q---------0
  %% Angle at p0-q-p1
  %thq = th(p1)-th(p0);
  %% Law of sines: θq/|01| = θ₀/|q1|
  %% θ₀ = |q1|*θq/|01|
  %th0 = l(p1)*thq/sqrt(sum((P(p0,:)-P(p1,:)).^2,2));
  % Vector from p0 to p1

  % Inputs:
  %   P  #P by 2 list of nearby points
  %   l  #P by 2 list of distances to origin: l = normrow(P)
  %   p0  first 2d point
  %   p1  second 2d point
  %   p2  third 2d point
  % Outputs:
  %   flag  whether edge from origin to p1 forms a Delaunay edge between points
  %     p1 and p2.
  % 
  isdel = @(P,l,p0,p1,p2) ...
    acos((-P(p0,:)/l(p0))*(((P(p1,:)-P(p0,:))/sqrt(sum((P(p1,:)-P(p0,:)).^2,2)))')) + ...
    acos((-P(p2,:)/l(p2))*(((P(p1,:)-P(p2,:))/sqrt(sum((P(p1,:)-P(p2,:)).^2,2)))')) <= pi;

  %% Or perhaps: http://stackoverflow.com/a/8523979
  %function f = isdel(P,l,p0,p1,p2)
  %  f = acos((-P(p0,:)/l(p0))*(((P(p1,:)-P(p0,:))/sqrt(sum((P(p1,:)-P(p0,:)).^2,2)))')) + ...
  %  acos((-P(p2,:)/l(p2))*(((P(p1,:)-P(p2,:))/sqrt(sum((P(p1,:)-P(p2,:)).^2,2)))')) <= pi;
  %  % Shewchuk's incircle test via determinant
  %  h(1,:) = [P(p0,:) 1];
  %  h(2,:) = [P(p1,:) 1];
  %  h(3,:) = [P(p2,:) 1];
  %  d = det([ ...
  %    q(1) q(2) q(1)^2+q(2)^2 1 ; ...
  %    h(1,1) h(1,2) h(1,1)^2+h(1,2)^2 h(1,3) ; ...
  %    h(2,1) h(2,2) h(2,1)^2+h(2,2)^2 h(2,3) ; ...
  %    h(3,1) h(3,2) h(3,1)^2+h(3,2)^2 h(3,3) ; ...
  %    ]);
  %  if f ~= (d<0)
  %    fprintf('%d %d\n',f,d<0);
  %    pause
  %  end