function d = d_measure ( dim_num, n, z, ns, sample_routine )

%*****************************************************************************80
%
%% d_measure() determines the pointset quality measure D.
%
%  Discussion:
%
%    The D measure of point distribution quality for a set Z of
%    N points in an DIM_NUM-dimensional region is defined as follows:
%
%    For each point Z(I) in the pointset, let V(I) be the subregion
%    defined by the intersection of the region with the Voronoi
%    region associated with Z(I).
%
%    Let D(I) be the determinant of the deviatoric tensor associated with
%    the region V(I).
%
%    Then D = maximum ( 1 <= I <= N ) D(I).
%
%    This quantity can be estimated using sampling.  A given number of
%    sample points are generated in the region, assigned to the nearest
%    element of the pointset, and used to approximate the Voronoi regions
%    and the deviatoric tensors.
%
%    In an ideally regular mesh, each deviatoric tensor would have a
%    zero determinant, and hence D would be zero.  In general, the smaller
%    D, the better.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 July 2021
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Max Gunzburger, John Burkardt,
%    Uniformity Measures for Point Samples in Hypercubes.
%
%  Input:
%
%    integer DIM_NUM, the spatial dimension.
%
%    integer N, the number of points.
%
%    real Z(N,DIM_NUM), the points.
%
%    integer NS, the number of sample points.
%
%    character SAMPLE_ROUTINE, the name of the function to be used to
%    return sample points in the region.  The function must have the form
%      x = sample_routine ( dim_num, n )
%
%  Output:
%
%    real D, the D quality measure.
%
  centroid(1:n,1:dim_num) = 0.0;
  hit(1:n) = 0;
  moment(1:dim_num,1:dim_num,1:n) = 0.0;

  for k = 1 : ns

    x = feval ( sample_routine, 1, dim_num );

    closest = find_closest ( dim_num, n, 1, x, z );

    hit(closest(1)) = hit(closest(1)) + 1;

    centroid(closest(1),1:dim_num) = centroid(closest(1),1:dim_num) + x(1,1:dim_num);

    for i1 = 1 : dim_num
      for i2 = 1 : dim_num
        moment(i1,i2,closest(1)) = moment(i1,i2,closest(1)) + x(i1) * x(i2);
      end
    end

  end

  for j = 1 : n

    if ( 0 < hit(j) )

      centroid(j,1:dim_num) = centroid(j,1:dim_num) / hit(j);

      moment(1:dim_num,1:dim_num,j) = moment(1:dim_num,1:dim_num,j) / hit(j);

      for i1 = 1 : dim_num
        for i2 = 1 : dim_num
          moment(i1,i2,j) = moment(i1,i2,j) - centroid(j,i1) * centroid(j,i2);
        end
      end

    end

  end

  tri(1:n) = 0.0;

  for j = 1 : n
    for i = 1 : dim_num
      tri(j) = tri(j) + moment(i,i,j);
    end
  end

  for j = 1 : n
    for i = 1 : dim_num
      moment(i,i,j) = moment(i,i,j) - tri(j) / dim_num;
    end
  end

  d = 0.0;

  for j = 1 : n

    di = det ( moment(1:dim_num,1:dim_num,j) );
    d = max ( d, di );

  end

  return
end
