function polvalues=polynomials_exact_order(points,order)
% generates all polynomials of order on points
[numpoints,dim]=size(points); % first some trivial cases
if order==1
    polvalues=ones(numpoints,1);
    return
end
if order==2
    polvalues=points;
    return
end
if dim==1
    polvalues=points(:,1).^(order-1);
    return
end
% What follows is a crude recursive scheme over the DIMENSION.
% Somebody MUST write a better one....
polvalues= points(:,dim).^(order-1);
for i=2:order-1
    pe=polynomials_exact_order(points(:,1:dim-1),i);
    pp=points(:,dim).^(order-i);
    [rpes cpes]=size(pe);
    % pps=size(points(:,dim).^(order-i));
    for j=1:cpes
        polvalues=[polvalues pe(:,j).*pp];
    end
end
polvalues=[polvalues polynomials_exact_order(points(:,1:dim-1),order)];