function [ r, theta ] = cartesian_to_hypersphere ( m, n, c, x )

%*****************************************************************************80
%
%% cartesian_to_hypersphere(): Cartesian to hypersphere coordinate transform.
%
%  Discussion:
%
%    We allow the trivial case M = 1; in that case alone, the value R
%    must be assumed to have a sign.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, the spatial dimension.
%    1 <= M.
%
%    integer N, the number of points to transform.
%
%    real C(M,1), the center of the hypersphere.
%
%    real X(M,N), the Cartesian coordinates of the points.
%
%  Output:
%
%    real R(N,1), the radius of the points on the hypersphere.
%    Except for the trivial case M = 1, R is assumed nonnegative.
%
%    real THETA(M-1,N), the coordinate angles of the points,
%    measured in radians.
%

%
%  Handle special case of M = 1.
%
  if ( m == 1 )
    r(1:n,1) = ( x(1,1:n) - c(1,1) )';
    theta = [];
    return
  end
%
%  Set aside memory.
%
  r = zeros ( n, 1 );
  theta = zeros ( m - 1, n );
%
%  Subtract the center.
%
  x(1:m,1:n) = x(1:m,1:n) - repmat ( c(1:m,1:1), 1, n );
  disp ( '1x' )
  disp ( x )
%
%  Compute R.
%
  r(1,1:n) = sqrt ( sum ( x(1:m,1:n).^2, 1 ) );
  disp ( '2r' )
  disp ( r )
%
%  Compute M-2 components of THETA.
%
  for i = 2 : m - 1
    theta(1:i-1,1:n) = theta(1:i-1,1:n) + repmat ( x(i,1:n).^2, i - 1, 1 );
  end
  disp ( '3theta' )
  disp ( theta )
  theta(1:m-2,1:n) = theta(1:m-2,1:n) + repmat ( x(m,1:n).^2, m - 2, 1 );
  disp ( '4theta' )
  disp ( theta )
  theta(1:m-2,1:n) = atan2 ( sqrt ( theta(1:m-2,1:n) ), x(1:m-2,1:n) );
  disp ( '5theta' )
  disp ( theta )
%
%  Compute last component of THETA.
%
  top(1,1:n) = ( sqrt ( x(m,1:n).^2 + x(m-1,1:n).^2 ) + x(m-1,1:n) );
  disp ( '6top' )
  disp ( top )
  theta(m-1,1:n) = 2.0 * atan2 ( x(m,1:n), top(1,1:n) );
  disp ( '7theta' )
  disp ( theta )

  return
end
