function [ point, weight ] = clenshaw_curtis_compute_nd ( dim_num, order_1d )

%% CLENSHAW_CURTIS_COMPUTE_ND returns a multidimensional Clenshaw-Curtis rule.
%
%  Discussion:
%
%    The value ORDER_ND, used to dimension the output arrays, is
%    simply the product of the entries in ORDER_1D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension of the points.
%
%    Input, integer ORDER_1D(DIM_NUM), the order of the Clenshaw-Curtis
%    rule in each dimension.
%
%    Output, real POINT(DIM_NUM,ORDER_ND), the points in
%    the grid.
%
%    Output, real WEIGHT(ORDER_ND), the integration weights
%    associated with the points.
%
  order_nd = prod ( order_1d(1:dim_num) );

  weight(1:order_nd) = 1.0;

  order = -1;

  for dim = 1 : dim_num

    order_old = order;

    order = order_1d(dim);
%
%  For efficiency's sake, we reuse the 1D rule if we can.
%
    if ( order ~= order_old )

      [ x1d, w1d ] = clenshaw_curtis_compute ( order );

    end

    p = 0;
    n1 = prod ( order_1d(1:dim-1) );
    n2 =        order_1d(dim);
    n3 = prod ( order_1d(dim+1:dim_num) );

    for k = 1 : n1
      for j = 1 : n2
        point(dim,p+1:p+n3) = x1d(j);
        weight(p+1:p+n3) = weight(p+1:p+n3) * w1d(j);
        p = p + n3;
      end
    end

  end

  x1d = [];
  w1d = [];

  return
end