function [ x, y, z ] = h_to_xyz ( h, r )

%*****************************************************************************80
%
%% h_to_xyz() converts a linear Hilbert 3D coordinate to (x,y,z).
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 August 2024
%
%  Author:
%
%    John Burkardt
%

%
%  Determine least significant octal digit of h.
%
  o = mod ( h, 8 );

  if ( o == 0 )
    x = 0;
    y = 0;
    z = 0;
  elseif ( o == 1 )
    x = 1;
    y = 0;
    z = 0;
  elseif ( o == 2 )
    x = 1;
    y = 0;
    z = 1;
  elseif ( o == 3 )
    x = 0;
    y = 0;
    z = 1;
  elseif ( o == 4 )
    x = 0;
    y = 1;
    z = 1;
  elseif ( o == 5 )
    x = 1;
    y = 1;
    z = 1;
  elseif ( o == 6 )
    x = 1;
    y = 1;
    z = 0;
  elseif ( o == 7 )
    x = 0;
    y = 1;
    z = 0;
  end

  w = 2;
  h = floor ( h / 8 );

  while ( 0 < h )

    o = mod ( h, 8 );

    xold = x;
    yold = y;
    zold = z;

    if ( o == 0 )
      x = yold;
      y = zold;
      z = xold;
    elseif ( o == 1 )
      x = zold + w;
      y = xold;
      z = yold;
    elseif ( o == 2 )
      x = zold + w;
      y = xold;
      z = yold + w;
    elseif ( o == 3 )
      x = w - xold - 1;
      y = yold;
      z = 2 * w - zold - 1;
    elseif ( o == 4 )
      x = w - xold - 1;
      y = yold + w;
      z = 2 * w - zold - 1;
    elseif ( o == 5 )
      x = zold + w;
      y = 2 * w - xold - 1;
      z = 2 * w - yold - 1;
    elseif ( o == 6 )
      x = zold + w;
      y = 2 * w - xold - 1;
      z = w - yold - 1;
    elseif ( o == 7 )
      x = w - yold - 1;
      y = 2 * w - zold - 1;
      z = xold;
    end
    h = floor ( h / 8 );
    w = w * w;
  end

  rm = rmin ( x, y, z );
  t = mod ( r - rm, 3 );

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

  return
end

