function dist = parallelogram_point_dist_3d ( p1, p2, p3, p )

%*****************************************************************************80
%
%% parallelogram_point_dist_3d(): distance ( parallelogram, point ) in 3D.
%
%  Diagram:
%
%            P2.................
%            /                 .
%           /                 .
%          /                 .
%        P1---------------->P3
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real P1(3,1), P2(3,1), P3(3,1), determine the
%    parallelogram, generated by the vectors from P1 to P2
%    and from P1 to P3.
%
%    real P(3,1), the point which is to be checked.
%
%  Output:
%
%    real DIST, the distance from the point to the
%    parallelogram.  DIST is zero if the point lies exactly on the
%    parallelogram.
%

%
%  Compute PP, the unit normal to X2-X1 and X3-X1:
%
  pp(1,1) = ( p2(2,1) - p1(2,1) ) * ( p3(3,1) - p1(3,1) ) ...
          - ( p2(3,1) - p1(3,1) ) * ( p3(2,1) - p1(2,1) );
  pp(2,1) = ( p2(3,1) - p1(3,1) ) * ( p3(1,1) - p1(1,1) ) ...
          - ( p2(1,1) - p1(1,1) ) * ( p3(3,1) - p1(3,1) );
  pp(3,1) = ( p2(1,1) - p1(1,1) ) * ( p3(2,1) - p1(2,1) ) ...
          - ( p2(2,1) - p1(2,1) ) * ( p3(1,1) - p1(1,1) );

  temp = sqrt ( sum ( pp(1:3,1).^2 ) );

  if ( temp == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARALLELOGRAM_POINT_DIST_3D - Fatal error!\n' );
    fprintf ( 1, '  The normal vector is zero.\n' );
    error ( 'PARALLELOGRAM_POINT_DIST_3D - Fatal error!' );
  end

  pp(1:3,1) = pp(1:3,1) / temp;
%
%  Find PN, the nearest point to P in the plane.
%
  t = pp(1:3,1)' * ( p(1:3,1) - p1(1:3,1) );

  pn(1:3,1) = p(1:3,1) - pp(1:3,1) * t;
%
%  If P lies WITHIN the parallelogram, we're done.
%
  inside = parallelogram_contains_point_3d ( p1, p2, p3, p );

  if ( inside )
    dist = sqrt ( sum ( ( pn(1:3,1) - p(1:3,1) ).^2 ) );
    return
  end
%
%  Otherwise, find the distance between P and each of the
%  four line segments that make up the boundary of the parallelogram.
%
  p4(1:3,1) = p2(1:3,1) + p3(1:3,1) - p1(1:3,1);

  dis13 = segment_point_dist_3d ( p1, p3, p );
  dis34 = segment_point_dist_3d ( p3, p4, p );
  dis42 = segment_point_dist_3d ( p4, p2, p );
  dis21 = segment_point_dist_3d ( p2, p1, p );

  dist = min ( dis13, min ( dis34, min ( dis42, dis21 ) ) );

  return
end
