function h = xyz_to_h ( x, y, z, r )

%*****************************************************************************80
%
%% xyz_to_h() converts a 3D Hilbert curve lattice point to a linear H index.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 August 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer x, y, z: values in [0,2^r) that are lattice coordinates.
%
%  Output:
%
%    integer h: the corresponding linear H index.
%
  rm = rmin ( x, y, z );
  t = mod ( r - rm, 3 );

  if ( t == 1 )
    xold = x;
    yold = y;
    zold = z;
    x = zold;
    y = xold;
    z = yold;
  elseif ( t == 2 )
    xold = x;
    yold = y;
    zold = z;
    x = yold;
    y = zold;
    z = xold;
  end

  h = 0;
  w = 2^( rm - 1 );

  for k = rm : -1 : 1

    xw = floor ( x / w );
    yw = floor ( y / w );
    zw = floor ( z / w );

    if (      [ xw, yw, zw ] == [ 0, 0, 0 ] )
      o = 0;
      xnew = z;
      ynew = x;
      znew = y;
    elseif ( [ xw, yw, zw ] == [ 1, 0, 0 ] )
      o = 1;
      xnew = y;
      ynew = z;
      znew = x - w;
    elseif ( [ xw, yw, zw ] == [ 1, 0, 1 ] )
      o = 2;
      xnew = y;
      ynew = z - w;
      znew = x - w;
    elseif ( [ xw, yw, zw ] == [ 0, 0, 1 ] )
      o = 3;
      xnew = w - x - 1;
      ynew = y;
      znew = 2 * w - z - 1;
    elseif ( [ xw, yw, zw ] == [ 0, 1, 1 ] )
      o = 4;
      xnew = w - x - 1;
      ynew = y - 2;
      znew = 2 * w - z - 1;
    elseif ( [ xw, yw, zw ] == [ 1, 1, 1 ] )
      o = 5;
      xnew = 2 * w - y - 1;
      ynew = 2 * w - z - 1;
      znew = x - w;
    elseif ( [ xw, yw, zw ] == [ 1, 1, 0 ] )
      o = 6;
      xnew = 2 * w - y - 1;
      ynew = w - z - 1;
      znew = x - w;
    elseif ( [ xw, yw, zw ] == [ 0, 1, 0 ] )
      o = 7;
      xnew = z;
      ynew = w - x - 1;
      znew = 2 * w - y - 1;
    end
%
%  Move to updated position.
%
    x = xnew;
    y = ynew;
    z = znew;

    h = 8 * h + o;
    w = floor ( w / 2 );

  end

  return
end

