function [x,y] = EquidistantNodes2D(N)

% compute (x,y) nodes in equilateral triangle for polynomial of order N

N_p = ((N + 1) * (N + 2)) / 2;

lambda1 = zeros(1,N_p); lambda2 = lambda1; lambda3 = lambda1;
temp = 1;
for i = 0:N
    for j = 0:(N - i)
        lambda1(temp) = i / N;
        lambda3(temp) = j / N;
        temp = temp + 1;
    end
end
lambda2 = 1 - lambda1 - lambda3;

x = -lambda2 + lambda3;
y = (-lambda2  - lambda3 + 2 * lambda1) / sqrt(3);

end