function cvt_3d_lumping_test ( )

%*****************************************************************************80
%
%% CVT_3D_LUMPING_TEST() tests CVT_3D_LUMPING().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    17 October 2016
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../cvt_3d_lumping' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_3D_LUMPING_TEST():\n' );
  fprintf ( 1, '  CVT_3D_LUMPING() applies a certain number of iterations\n' );
  fprintf ( 1, '  to an initial set of N points in the unit cube,\n' );
  fprintf ( 1, '  using a sampling approach to estimate centroids,\n' );
  fprintf ( 1, '  and a weight function which affects the centroids,\n' );
  fprintf ( 1, '  in an attempt to locate an arrangement of the points\n' );
  fprintf ( 1, '  that satisfies the CVT condition.\n' );
  fprintf ( 1, '\n' );
%
%  We are most interested in values of N that correspond to the number of
%  monomials in X, Y, and Z up to some degree, that is, values like 
%  Degree =  0,  1,  2,  3,  4,  5, ...
%  N      = (1), 4, 10, 20, 35, 56,...
%
  for degree = 1 : 10

    n = ( ( degree + 1 ) * ( degree + 2 ) * ( degree + 3 ) ) / 6;

    it_num = 200;
    s_num = 100;
    init = 1;

    xyz = cvt_3d_lumping ( n, it_num, s_num, @mu_3d_chebyzero, init );
%
%  Write values to a file.
%
    filename = sprintf ( 'generators_%03d.txt', n );
    unit = fopen ( filename, 'wt' );

    for i = 1 : n
      for j = 1 : 3
        fprintf ( unit, '  %g', xyz(i,j) );
      end
      fprintf ( unit, '\n' );
    end
    fclose ( unit );
%
%  Plot the final configuration of generators.
%  Save it as a ".fig" file that can be reopened and manipulated with MATLAB.
%
    figure ( );
    plot3 ( xyz(:,1), xyz(:,2), xyz(:,3), 'b.', 'Markersize', 15 );
    grid ( 'on' );
    title ( 'Generator locations:' );
    filename = sprintf ( 'generators_%03d', n );
    savefig ( filename );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  3D generator location plot saved as "%s.fig"\n', filename );
    fprintf ( 1, '  Reopen it with MATLAB for viewing, printing, or manipulation.\n' );
%
%  The Delaunay tetrahedrons can be plotted by calling
%    T = delaunay ( XYZ )
%    TETRAMESH ( T, XYZ )
%
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_3D_LUMPING_TEST()::\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../cvt_3d_lumping' );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

