function [u,v,w]=cnodes(N)
%
% function [u,v,w]=cnodes(N)
%
% coordinates of 2*(N+1)^2 points on an octahedron
%
% u+v+|w|=1 on two faces of the octahedron
%
% 'computational' points are uniformly distributed on the faces
%
% used by other functions that find maps of these that are
% approximately uniformly distributed on the sphere
%

if nargin < 1 help cnodes; return; end

M=N+1;

w=(1/N)*(0:N)'*ones(1,M);           % N+1 cols of 0:1
v=w; u=w';                          % same old tricks
v(u<v)=u(u<v);                      % 0:1 on diag decreasing toward edges
u=reshape(v(M^2:-1:1),size(v));     % faster than fliplr(flipud(v))
w=w'-w;                             % -1:1 on diagonals, increasing on goofys