function values = ModalBaseValues2D(a,b,i,j)

% evaluates 2D orthonormal polynomial of order (i,j) on simplex at (a,b)

values1 = JacobiP(a,0,0,i);
values2 = JacobiP(b,2 * i + 1,0,j);
values = sqrt(2) * values1 .* values2 .* ((1 - b) .^ i);

end