function [ pc, normal ] = circle_pppr2imp_3d ( p1, p2, p3, r )

%*****************************************************************************80
%
%% circle_pppr2imp_3d() converts a circle from PPR to implicit form in 3D.
%
%  Discussion:
%
%    The PPPR form of a circle in 3D is:
%
%      The circle of radius R passing through points P1 and P2,
%      and lying in the plane of P1, P2 and P3.
%
%    The implicit form of a circle in 3D is:
%
%      ( P(1) - PC(1) )^2 + ( P(2) - PC(2) )^2 + ( P(3) - PC(3) )^2 = R^2
%      and the dot product of P - PC with NORMAL is 0.
%
%    There may be zero, one, or two circles that satisfy the
%    requirements of the PPPR form.
%
%    If there is no such circle, then PC(1:2,1) and PC(1:2,2)
%    are set to the midpoint of (P1,P2).
%
%    If there is one circle, PC(1:2,1) and PC(1:2,2) will be equal.
%
%    If there are two circles, then PC(1:2,1) is the first center,
%    and PC(1:2,2) is the second.
%
%    This calculation is equivalent to finding the intersections of
%    spheres of radius R at points P1 and P2, which lie in the plane
%    defined by P1, P2 and P3.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 March 2006
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real P1(3,1), P2(3,1), two points on the circle.
%
%    real P3(3,1), a third point.
%
%    real R, the radius of the circle.
%
%  Output:
%
%    real PC(3,2), the centers of the two circles.
%
%    real NORMAL(3,1), the normal to the circles.
%
  p1 = p1(:);
  p2 = p2(:);
  p3 = p3(:);

  dim_num = 3;
%
%  Compute the distance from P1 to P2.
%
  dist = sqrt ( sum ( ( p2(1:dim_num,1) - p1(1:dim_num,1) ).^2 ) );
%
%  If R is smaller than DIST, we don't have a circle.
%
  if ( 2.0 * r < dist )
    for j = 1 : 2
      pc(1:dim_num,j) = 0.5 * ( p1(1:dim_num,1) + p2(1:dim_num,1) );
    end
    return
  end
%
%  H is the distance from the midpoint of (P1,P2) to the center.
%
  h = sqrt ( ( r + 0.5 * dist ) * ( r - 0.5 * dist ) );
%
%  Define a unit direction V that is normal to P2-P1, and lying
%  in the plane (P1,P2,P3).
%
%  To do this, subtract from P3-P1 the component in the direction P2-P1.
%
  v(1:dim_num,1) = p3(1:dim_num,1) - p1(1:dim_num,1);
  dot = v(1:dim_num,1)' * ( p2(1:dim_num,1) - p1(1:dim_num,1) );
  dot = dot / dist;

  v(1:dim_num,1) = v(1:dim_num,1) - dot * ( p2(1:dim_num,1) - p1(1:dim_num,1) ) / dist;

  v_length = sqrt ( sum ( v(1:dim_num,1).^2 ) );

  v(1:dim_num,1) = v(1:dim_num,1) / v_length;
%
%  We can go with or against the given normal direction.
%
  pc(1:dim_num,1) = 0.5 * ( p2(1:dim_num,1) + p1(1:dim_num,1) ) + h * v(1:dim_num,1);

  pc(1:dim_num,2) = 0.5 * ( p2(1:dim_num,1) + p1(1:dim_num,1) ) - h * v(1:dim_num,1);

  normal = plane_exp_normal ( p1, p2, p3 );

  return
end
