function [x,y] = OptimalNodes2D(N)

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

alpha_opt = [0.0000, 0.0000, 1.4152, 0.1001, 0.2751, 0.9800, 1.0999, 1.2832, 1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258];

if (N < 16)
    alpha_nodes = alpha_opt(N);
else
    alpha_nodes = 5/3;
end;

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);

% jak se tyto vztahy dostaly???
blend1 = 4 * lambda2 .* lambda3;
blend2 = 4 * lambda1 .* lambda3;
blend3 = 4 * lambda1 .* lambda2;

w_edge1 = Warpfactor(N,lambda3 - lambda2);
w_edge2 = Warpfactor(N,lambda1 - lambda3);
w_edge3 = Warpfactor(N,lambda2 - lambda1);

w_final1 = (1 + (alpha_nodes * lambda1) .^ 2) .* blend1 .* w_edge1;
w_final2 = (1 + (alpha_nodes * lambda2) .^ 2) .* blend2 .* w_edge2;
w_final3 = (1 + (alpha_nodes * lambda3) .^ 2) .* blend3 .* w_edge3;

x = x + 1 * w_final1 + cos(2 * pi / 3) * w_final2 + cos(4 * pi / 3) * w_final3;
y = y + 0 * w_final1 + sin(2 * pi / 3) * w_final2 + sin(4 * pi / 3) * w_final3;

end