# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <fstream>
# include <cmath>
# include <ctime>

using namespace std;

# include "geometry.H"

//****************************************************************************80

void angle_box_2d ( double dist, double p1[2], double p2[2], double p3[2], 
  double p4[2], double p5[2] )

//****************************************************************************80
//
//  Purpose:
//
//    ANGLE_BOX_2D "boxes" an angle defined by three points in 2D.
//
//  Discussion:
//
//    The routine is given points P1, P2 and P3, determining the two lines:
//      P1 to P2
//    and
//      P2 to P3
//    and a nonnegative distance
//      DIST.
//
//    The routine returns a pair of "corner" points
//      P4 and P5
//    both of which are a distance DIST from both lines, and in fact,
//    both of which are a distance DIST from P2.
//
//                         /  P3
//                        /   /   /
//     - - - - - - - - -P4 - / -P6 - - -
//                      /   /   /
//    P1---------------/--P2-----------------
//                    /   /   /
//     - - - - - - -P7 - / -P5 - - - - -
//                  /   /   /
//
//    In the illustration, P1, P2 and P3 represent
//    the points defining the lines.
//
//    P4 and P5 represent the desired "corner points", which
//    are on the positive or negative sides of both lines.
//
//    The numbers P6 and P7 represent the undesired points, which
//    are on the positive side of one line and the negative of the other.
//
//    Special cases:
//
//    if P1 = P2, this is the same as extending the line from
//    P3 through P2 without a bend.
//
//    if P3 = P2, this is the same as extending the line from
//    P1 through P2 without a bend.
//
//    if P1 = P2 = P3 this is an error.
//
//  Modified:
//
//    17 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double DIST, the nonnegative distance from P1
//    to the computed points P4 and P5.
//
//    Input, double P1[2], P2[2], P3[2].
//    P1 and P2 are distinct points that define a line.
//    P2 and P3 are distinct points that define a line.
//
//    Output, double P4[2], P5[2], points which lie DIST units from
//    the line between P1 and P2, and from the line between P2 and P3.
//
{
# define DIM_NUM 2

  double stheta;
  double temp;
  double temp1;
  double temp2;
  double u[DIM_NUM];
  double u1[DIM_NUM];
  double u2[DIM_NUM];
//
//  If DIST = 0, assume the user knows best.
//
  if ( dist == 0.0 )
  {
    r8vec_copy ( DIM_NUM, p2, p4 );
    r8vec_copy ( DIM_NUM, p2, p5 );
    return;
  }
//
//  Fail if all three points are equal.
//
  if ( r8vec_eq ( DIM_NUM, p1, p2 ) && r8vec_eq ( DIM_NUM, p2, p3 ) )
  {
    cout << "\n";
    cout << "ANGLE_BOX_2D - Fatal error!\n";
    cout << "  Input points P3 = P2 = P3.\n";
    r8vec_print ( DIM_NUM, p1, "  P1:" );
    exit ( 1 );
  }
//
//  If P1 = P2, extend the line through the doubled point.
//
  if ( r8vec_eq ( DIM_NUM, p1, p2 ) )
  {
    u2[0] = p3[1] - p2[1];
    u2[1] = p2[0] - p3[0];
    temp = r8vec_length ( DIM_NUM, u2 );
    u2[0] = u2[0] / temp;
    u2[1] = u2[1] / temp;
    p4[0] = p2[0] + dist * u2[0];
    p4[1] = p2[1] + dist * u2[1];
    p5[0] = p2[0] - dist * u2[0];
    p5[1] = p2[1] - dist * u2[1];
    return;
  }
//
//  If P2 = P3, extend the line through the doubled point.
//
  if ( r8vec_eq ( DIM_NUM, p2, p3 ) )
  {
    u1[0] = p1[1] - p2[1];
    u1[1] = p2[0] - p1[0];
    temp = r8vec_length ( DIM_NUM, u1 );
    u1[0] = u1[0] / temp;
    u1[1] = u1[1] / temp;
    p4[0] = p2[0] + dist * u1[0];
    p4[1] = p2[1] + dist * u1[1];
    p5[0] = p2[0] - dist * u1[0];
    p5[1] = p2[1] - dist * u1[1];
    return;
  }
//
//  Now compute the unit normal vectors to each line.
//  We choose the sign so that the unit normal to line 1 has
//  a positive dot product with line 2.
//
  u1[0] = p1[1] - p2[1];
  u1[1] = p2[0] - p1[0];
  temp = r8vec_length ( DIM_NUM, u1 );
  u1[0] = u1[0] / temp;
  u1[1] = u1[1] / temp;

  temp1 = u1[0] * ( p3[0] - p2[0] )
        + u1[1] * ( p3[1] - p2[1] );

  if ( temp1 < 0.0 )
  {
    u1[0] = -u1[0];
    u1[1] = -u1[1];
  }

  u2[0] = p3[1] - p2[1];
  u2[1] = p2[0] - p3[0];
  temp = r8vec_length ( DIM_NUM, u2 );
  u2[0] = u2[0] / temp;
  u2[1] = u2[1] / temp;

  temp2 = u2[0] * ( p1[0] - p2[0] )
        + u2[1] * ( p1[1] - p2[1] );

  if ( temp2 < 0.0 )
  {
    u2[0] = -u2[0];
    u2[1] = -u2[1];
  }
//
//  Try to catch the case where we can't determine the
//  sign of U1, because both U1 and -U1 are perpendicular
//  to (P3-P2)...and similarly for U2 and (P1-P2).
//
  if ( temp1 == 0.0 || temp2 == 0.0 )
  {
    if ( u1[0] * u2[0] + u1[1] * u2[1] < 0.0 )
    {
      u1[0] = -u1[0];
      u2[0] = -u2[0];
    }
  }
//
//  Try to catch a line turning back on itself, evidenced by
//    Cos(theta) = (P3-P2) dot (P2-P1) / ( norm(P3-P2) * norm(P2-P1) )
//  being -1, or very close to -1.
//
  temp = ( p3[0] - p2[0] ) * ( p2[0] - p1[0] )
       + ( p3[1] - p2[1] ) * ( p2[1] - p1[1] );

  temp1 = sqrt ( pow ( p3[0] - p2[0], 2 ) + pow ( p3[1] - p2[1], 2 ) );
  temp2 = sqrt ( pow ( p2[0] - p1[0], 2 ) + pow ( p2[1] - p1[1], 2 ) );

  temp = temp / ( temp1 * temp2 );

  if ( temp < -0.99 )
  {
    temp = sqrt ( pow ( p2[0] - p1[0], 2 ) + pow ( p2[1] - p1[1], 2 ) );

    p4[0] = p2[0] + dist * ( p2[0] - p1[0] ) / temp + dist * u1[0];
    p4[1] = p2[1] + dist * ( p2[1] - p1[1] ) / temp + dist * u1[1];
    p5[0] = p2[0] + dist * ( p2[0] - p1[0] ) / temp - dist * u1[0];
    p5[1] = p2[1] + dist * ( p2[1] - p1[1] ) / temp - dist * u1[1];
    return;
  }
//
//  Compute the "average" unit normal vector.
//
//  The average of the unit normals could be zero, but only when
//  the second line has the same direction and opposite sense
//  of the first, and we've already checked for that case.
//
//  Well, check again!  This problem "bit" me in the case where
//  P1 = P2, which I now treat specially just to guarantee I
//  avoid this problem!
//
  if ( u1[0] * u2[0] + u1[1] * u2[1] < 0.0 )
  {
    u2[0] = -u2[0];
    u2[1] = -u2[1];
  }

  u[0] = 0.5 * ( u1[0] + u2[0] );
  u[1] = 0.5 * ( u1[1] + u2[1] );

  temp = r8vec_length ( DIM_NUM, u );
  u[0] = u[0] / temp;
  u[1] = u[1] / temp;
//
//  You must go DIST/STHETA units along this unit normal to
//  result in a distance DIST from line1 (and line2).
//
  stheta = u[0] * u1[0] + u[1] * u1[1];

  p4[0] = p2[0] + dist * u[0] / stheta;
  p4[1] = p2[1] + dist * u[1] / stheta;
  p5[0] = p2[0] - dist * u[0] / stheta;
  p5[1] = p2[1] - dist * u[1] / stheta;

  return;
# undef DIM_NUM
}
//****************************************************************************80

bool angle_contains_ray_2d ( double p1[2], double p2[2], double p3[2], 
  double p[2] )

//****************************************************************************80
//
//  Purpose:
//
//    ANGLE_CONTAINS_RAY_2D determines if an angle contains a ray, in 2D.
//
//  Discussion:
//
//    The angle is defined by the sequence of points P1, P2, P3.
//
//    The ray is defined by the sequence of points P2, P.
//
//  Modified:
//
//    30 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], P3[2], the coordinates of points on
//    the angle.
//
//    Input, double P[2], the end point of the ray to be checked.
//    The ray is assumed to have an origin at P2.
//
//    Output, bool ANGLE_CONTAINS_RAY_2D, is true if the ray is inside
//    the angle or on its boundary, and false otherwise.
//
{
  double a1;
  double a2;

  a1 = angle_deg_2d ( p1, p2, p );

  a2 = angle_deg_2d ( p1, p2, p3 );

  if ( a1 <= a2 )
  {
    return true;
  } 
  else
  {
    return false;
  }

}
//****************************************************************************80

double angle_deg_2d ( double p1[2], double p2[2], double p3[2] )

//****************************************************************************80
//
//  Purpose:
//
//    ANGLE_DEG_2D returns the angle in degrees swept out between two rays in 2D.
//
//  Discussion:
//
//    Except for the zero angle case, it should be true that
//
//    ANGLE_DEG_2D(P1,P2,P3) + ANGLE_DEG_2D(P3,P2,P1) = 360.0
//
//  Modified:
//
//    30 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], P3[2], define the rays
//    P1 - P2 and P3 - P2 which define the angle.
//
//    Output, double ANGLE_DEG_2D, the angle swept out by the rays, measured
//    in degrees.  0 <= ANGLE_DEG_2D < 360.  If either ray has zero length,
//    then ANGLE_DEG_2D is set to 0.
//
{
# define DIM_NUM 2

  double angle_rad;
  double p[DIM_NUM];
  double pi = 3.141592653589793;
  double value;

  p[0] = ( p1[0] - p2[0] ) * ( p3[0] - p2[0] ) 
       + ( p1[1] - p2[1] ) * ( p3[1] - p2[1] );

  p[1] = ( p1[0] - p2[0] ) * ( p3[1] - p2[1] ) 
       - ( p1[1] - p2[1] ) * ( p3[0] - p2[0] );

  if ( p[0] == 0.0 && p[1] == 0.0 )
  {
    value = 0.0;
  }
  else
  {
    angle_rad = atan2 ( p[1], p[0] );

    if ( angle_rad < 0.0 )
    {
      angle_rad = angle_rad + 2.0 * pi;
    }

    value = radians_to_degrees ( angle_rad );

  }

  return value;
# undef DIM_NUM
}
//****************************************************************************80

double *angle_half_2d ( double p1[2], double p2[2], double p3[2] )

//****************************************************************************80
//
//  Purpose:
//
//    ANGLE_HALF_2D finds half an angle in 2D.
//
//  Discussion:
//
//    The original angle is defined by the sequence of points P1, P2 and P3.
//
//    The point P4 is calculated so that:
//
//      (P1,P2,P4) = (P1,P2,P3) / 2
//
//        P1
//        /
//       /   P4
//      /  .
//     / .
//    P2--------->P3
//
//  Modified:
//
//    22 May 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], P3[2], points defining the angle.
//
//    Input, double ANGLE_HALF_2D[2], a point P4 defining the half angle.
//    The vector P4 - P2 will have unit norm.
//
{
# define DIM_NUM 2

  int i;
  double norm;
  double *p4;

  p4 = new double[DIM_NUM];

  norm = sqrt ( ( p1[0] - p2[0] ) * ( p1[0] - p2[0] )
              + ( p1[1] - p2[1] ) * ( p1[1] - p2[1] ) );

  for ( i = 0; i < DIM_NUM; i++ )
  {
    p4[i] = ( p1[i] - p2[i] ) / norm;
  }

  norm = sqrt ( ( p3[0] - p2[0] ) * ( p3[0] - p2[0] )
              + ( p3[1] - p2[1] ) * ( p3[1] - p2[1] ) );

  for ( i = 0; i < DIM_NUM; i++ )
  {
    p4[i] = p4[i] + ( p3[i] - p2[i] ) / norm;
  }

  for ( i = 0; i < DIM_NUM; i++ )
  {
    p4[i] = 0.5 * p4[i];
  }

  norm = r8vec_length ( DIM_NUM, p4 );
  
  for ( i = 0; i < DIM_NUM; i++ )
  {
    p4[i] = p2[i] + p4[i] / norm;
  }

  return p4;
# undef DIM_NUM
}
//****************************************************************************80

double angle_rad_2d ( double p1[2], double p2[2], double p3[2] )

//****************************************************************************80
//
//  Purpose:
//
//    ANGLE_RAD_2D returns the angle in radians swept out between two rays in 2D.
//
//  Discussion:
//
//      ANGLE_RAD_2D ( P1, P2, P3 ) + ANGLE_RAD_2D ( P3, P2, P1 ) = 2 * PI
//
//        P1
//        /
//       /    
//      /     
//     /  
//    P2--------->P3
//
//  Modified:
//
//    24 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], P3[2], define the rays
//    P1 - P2 and P3 - P2 which define the angle.
//
//    Output, double ANGLE_RAD_3D, the angle between the two rays,
//    in radians.  This value will always be between 0 and 2*PI.  If either ray has
//    zero length, then the angle is returned as zero.
//
{
# define DIM_NUM 2

  double p[DIM_NUM];
  double pi = 3.141592653589793;
  double value;

  p[0] = ( p3[0] - p2[0] ) * ( p1[0] - p2[0] ) 
       + ( p3[1] - p2[1] ) * ( p1[1] - p2[1] );


  p[1] = ( p3[0] - p2[0] ) * ( p1[1] - p2[1] ) 
       - ( p3[1] - p2[1] ) * ( p1[0] - p2[0] );

  if ( p[0] == 0.0 && p[1] == 0.0 )
  {
    value = 0.0;
    return value;
  }

  value = atan2 ( p[1], p[0] );

  if ( value < 0.0 )
  {
    value = value + 2.0 * pi;
  }

  return value;
# undef DIM_NUM
}
//****************************************************************************80

double angle_rad_3d ( double p1[3], double p2[3], double p3[3] )

//****************************************************************************80
//
//  Purpose:
//
//    ANGLE_RAD_3D returns the angle between two vectors in 3D.
//
//  Discussion:
//
//    The routine always computes the SMALLER of the two angles between
//    two vectors.  Thus, if the vectors make an (exterior) angle of 200 
//    degrees, the (interior) angle of 160 is reported.
//
//    X dot Y = Norm(X) * Norm(Y) * Cos ( Angle(X,Y) )
//
//  Modified:
//
//    20 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[3], P2[3], P3[3], points defining an angle.
//    The rays are P1 - P2 and P3 - P2.
//
//    Output, double ANGLE_RAD_3D, the angle between the two vectors, in radians.
//    This value will always be between 0 and PI.  If either vector has
//    zero length, then the angle is returned as zero.
//
{
# define DIM_NUM 3

  double dot;
  int i;
  double v1norm;
  double v2norm;
  double value;
  
  v1norm = 0.0;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    v1norm = v1norm + pow ( p1[i] - p2[i], 2 );
  }
  v1norm = sqrt ( v1norm );

  if ( v1norm == 0.0 )
  {
    value = 0.0;
    return value;
  }

  v2norm = 0.0;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    v2norm = v2norm + pow ( p3[i] - p2[i], 2 );
  }
  v2norm = sqrt ( v2norm );

  if ( v2norm == 0.0 )
  {
    value = 0.0;
    return value;
  }

  dot = 0.0;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    dot = dot + ( p1[i] - p2[i] ) * ( p3[i] - p2[i] );
  }

  value = arc_cosine ( dot / ( v1norm * v2norm ) );

  return value;
# undef DIM_NUM
}
//****************************************************************************80

double angle_rad_nd ( int dim_num, double vec1[], double vec2[] )

//****************************************************************************80
//
//  Purpose:
//
//    ANGLE_RAD_ND returns the angle between two vectors in ND.
//
//  Discussion:
//
//    ANGLE_RAD_ND always computes the SMALLER of the two angles between
//    two vectors.  Thus, if the vectors make an (exterior) angle of 
//    1.5 PI radians, then the (interior) angle of 0.5 radians is returned.
//
//    X dot Y = Norm(X) * Norm(Y) * Cos( Angle(X,Y) )
//
//  Modified:
//
//    19 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, double VEC1[DIM_NUM], VEC2[DIM_NUM], the two vectors to be considered.
//
//    Output, double ANGLE_RAD_ND, the angle between the vectors, in radians.
//    This value will always be between 0 and PI.
//
{
  double dot;
  double v1norm;
  double v2norm;
  double value;

  dot = r8vec_dot ( dim_num, vec1, vec2 );

  v1norm = r8vec_length ( dim_num, vec1 );
  v2norm = r8vec_length ( dim_num, vec2 );

  if ( v1norm == 0.0 || v2norm == 0.0 )
  {
    value = 0.0;
  }
  else
  {
    value = acos ( dot / ( v1norm * v2norm ) );
  }

  return value;
}
//****************************************************************************80

double angle_turn_2d ( double p1[2], double p2[2], double p3[2] )

//****************************************************************************80
//
//  Purpose:
//
//    ANGLE_TURN_2D computes a turning angle in 2D.
//
//  Discussion:
//
//    This routine is most useful when considering the vertices of a
//    polygonal shape.  We wish to distinguish between angles that "turn
//    in" to the shape, (between 0 and 180 degrees) and angles that
//    "turn out" (between 180 and 360 degrees), as we traverse the boundary.
//
//    If we compute the interior angle and subtract 180 degrees, we get the
//    supplementary angle, which has the nice property that it is
//    negative for "in" angles and positive for "out" angles, and is zero if
//    the three points actually lie along a line.
//
//    Assuming P1, P2 and P3 define an angle, the TURN can be
//    defined to be either:
//
//    * the supplementary angle to the angle formed by P1=P2=P3, or
//
//    * the angle between the vector ( P3-P2) and the vector -(P1-P2),
//      where -(P1-P2) can be understood as the vector that continues
//      through P2 from the direction P1.
//
//    The turning will be zero if P1, P2 and P3 lie along a straight line.
//
//    It will be a positive angle if the turn from the previous direction
//    is counter clockwise, and negative if it is clockwise.
//
//    The turn is given in radians, and will lie between -PI and PI.
//
//  Modified:
//
//    14 March 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], P3[2], the points that form
//    the angle.
//
//    Output, double ANGLE_TURN_2D, the turn angle, between -PI and PI.
//
{
# define DIM_NUM 2

  double p[DIM_NUM];
  double pi = 3.141592653589793;
  double turn;

  p[0] = ( p3[0] - p2[0] ) * ( p1[0] - p2[0] ) 
       + ( p3[1] - p2[1] ) * ( p1[1] - p2[1] );

  p[1] = ( p3[0] - p2[0] ) * ( p1[1] - p2[1] ) 
       - ( p3[1] - p2[1] ) * ( p1[0] - p2[0] );

  if ( p[0] == 0.0 && p[1] == 0.0 )
  {
    turn = 0.0;
  }
  else
  {
    turn = pi - atan4 ( p[1], p[0] );
  }

  return turn;
# undef DIM_NUM
}
//****************************************************************************80

double anglei_deg_2d ( double p1[2], double p2[2], double p3[2] )

//****************************************************************************80
//
//  Purpose:
//
//    ANGLEI_DEG_2D returns the interior angle in degrees between two rays in 2D.
//
//  Modified:
//
//    30 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[3], P2[3], P3[3], points defining an angle.
//    The rays are P1 - P2 and P3 - P2.
//
//    Output, double ANGLEI_DEG_2D, the angle swept out by the rays, measured
//    in degrees.  This value satisfies 0 <= ANGLEI_DEG_2D < 180.0.  If either
//    ray is of zero length, then ANGLEI_deg_2D is returned as 0.
//
{
# define DIM_NUM 2

  double p[DIM_NUM];
  double pi = 3.141592653589793;
  double value;

  p[0] = ( p1[0] - p2[0] ) * ( p3[0] - p2[0] )
       + ( p1[1] - p2[1] ) * ( p3[1] - p2[1] );

  p[1] = ( p1[0] - p2[0] ) * ( p3[1] - p2[1] ) 
       - ( p1[1] - p2[1] ) * ( p3[0] - p2[0] );

  if ( p[0] == 0.0 && p[1] == 0.0 )
  {
    value = 0.0;
  }
  else
  {
    value = atan2 ( p[1], p[0] );

    if ( value < 0.0 )
    {
      value = value + 2.0 * pi;
    }

    value = radians_to_degrees ( value );

    if ( 180.0 < value )
    {
      value = 360.0 - value;
    }

  }
  return value;
# undef DIM_NUM
}
//****************************************************************************80

double anglei_rad_2d ( double p1[2], double p2[2], double p3[2] )

//****************************************************************************80
//
//  Purpose:
//
//    ANGLEI_RAD_2D returns the interior angle in radians between two rays in 2D.
//
//  Modified:
//
//    30 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[3], P2[3], P3[3], points defining an angle.
//    The rays are P1 - P2 and P3 - P2.
//
//    Output, double ANGLEI_RAD_2D, the angle swept out by the rays, measured
//    in degrees.  This value satisfies 0 <= ANGLEI_RAD_2D < PI.  If either
//    ray is of zero length, then ANGLEI_RAD_2D is returned as 0.
//
{
# define DIM_NUM 2

  double p[DIM_NUM];
  double pi = 3.141592653589793;
  double value;

  p[0] = ( p1[0] - p2[0] ) * ( p3[0] - p2[0] )
       + ( p1[1] - p2[1] ) * ( p3[1] - p2[1] );

  p[1] = ( p1[0] - p2[0] ) * ( p3[1] - p2[1] ) 
       - ( p1[1] - p2[1] ) * ( p3[0] - p2[0] );

  if ( p[0] == 0.0 && p[1] == 0.0 )
  {
    value = 0.0;
  }
  else
  {
    value = atan2 ( p[1], p[0] );

    if ( value < 0.0 )
    {
      value = value + 2.0 * pi;
    }

    if ( pi < value )
    {
      value = 2.0 * pi - value;
    }

  }
  return value;
# undef DIM_NUM
}
//****************************************************************************80

double annulus_area_2d ( double r1, double r2 )

//****************************************************************************80
//
//  Purpose:
//
//    ANNULUS_AREA_2D computes the area of a circular annulus in 2D.
//
//  Discussion:
//
//    A circular annulus with center (XC,YC), inner radius R1 and
//    outer radius R2, is the set of points (X,Y) so that
//
//      R1**2 <= (X-XC)**2 + (Y-YC)**2 <= R2**2
//
//  Modified:
//
//    02 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R1, R2, the inner and outer radii.
//
//    Output, double ANNULUS_AREA_2D, the area.
//
{
  double area;
  double pi = 3.141592653589793;

  area = pi * ( r2 + r1 ) * ( r2 - r1 );

  return area;
}
//****************************************************************************80

double annulus_sector_area_2d ( double r1, double r2, double theta1, 
  double theta2 )

//****************************************************************************80
//
//  Purpose:
//
//    ANNULUS_SECTOR_AREA_2D computes the area of an annular sector in 2D.
//
//  Discussion:
//
//    An annular sector with center PC, inner radius R1 and
//    outer radius R2, and angles THETA1, THETA2, is the set of points 
//    P so that
//
//      R1**2 <= (P(1)-PC(1))**2 + (P(2)-PC(2))**2 <= R2**2
//
//    and
//
//      THETA1 <= THETA ( P - PC ) <= THETA2
//
//  Modified:
//
//    02 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R1, R2, the inner and outer radii.
//
//    Input, double THETA1, THETA2, the angles.
//
//    Output, double ANNULUS_SECTOR_AREA_2D, the area.
//
{
  double area;

  area = 0.5 * ( theta2 - theta1 ) * ( r2 + r1 ) * ( r2 - r1 );

  return area;
}
//****************************************************************************80

double *annulus_sector_centroid_2d ( double pc[2], double r1, double r2, 
  double theta1, double theta2 )

//****************************************************************************80
//
//  Purpose:
//
//    ANNULUS_SECTOR_CENTROID_2D computes the centroid of an annular sector in 2D.
//
//  Discussion:
//
//    An annular sector with center PC, inner radius R1 and
//    outer radius R2, and angles THETA1, THETA2, is the set of points 
//    P so that
//
//      R1**2 <= (P(1)-PC(1))**2 + (P(2)-PC(2))**2 <= R2**2
//
//    and
//
//      THETA1 <= THETA ( P - PC ) <= THETA2
//
//    Thanks to Ed Segall for pointing out a mistake in the computation
//    of the angle THETA assciated with the centroid.
//
//  Modified:
//
//    02 December 2005
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John Harris, Horst Stocker,
//    Handbook of Mathematics and Computational Science,
//    Springer, 1998, QA40.S76
//
//  Parameters:
//
//    Input, double PC[2], the center.
//
//    Input, double R1, R2, the inner and outer radii.
//
//    Input, double THETA1, THETA2, the angles.
//
//    Output, double ANNULUS_SECTOR_CENTROID_2D[2], the centroid.
//
{
  double *centroid;
  double r;
  double theta;

  theta = theta2 - theta1;

  r = 4.0 * sin ( theta / 2.0 ) / ( 3.0 * theta )
    * ( r1 * r1 + r1 * r2 + r2 * r2 ) / ( r1 + r2 );

  centroid = new double[2];

  centroid[0] = pc[0] + r * cos ( theta1 + theta / 2.0 );
  centroid[1] = pc[1] + r * sin ( theta1 + theta / 2.0 );

  return centroid;
}
//****************************************************************************80

double arc_cosine ( double c )

//****************************************************************************80
//
//  Purpose:
//
//    ARC_COSINE computes the arc cosine function, with argument truncation.
//
//  Discussion:
//
//    If you call your system ACOS routine with an input argument that is
//    outside the range [-1.0, 1.0 ], you may get an unpleasant surprise.
//    This routine truncates arguments outside the range.
//
//  Modified:
//
//    13 June 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double C, the argument, the cosine of an angle.
//
//    Output, double ARC_COSINE, an angle whose cosine is C.
//
{
  double angle;
  double pi = 3.141592653589793;

  if ( c <= -1.0 )
  {
    angle = pi;
  } 
  else if ( 1.0 <= c )
  {
    angle = 0.0;
  }
  else
  {
    angle = acos ( c );
  }
  return angle;
}
//****************************************************************************80

double arc_sine ( double s )

//****************************************************************************80
//
//  Purpose:
//
//    ARC_SINE computes the arc sine function, with argument truncation.
//
//  Discussion:
//
//    If you call your system ASIN routine with an input argument that is
//    outside the range [-1.0, 1.0 ], you may get an unpleasant surprise.
//    This routine truncates arguments outside the range.
//
//  Modified:
//
//    13 June 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double S, the argument, the sine of an angle.
//
//    Output, double ARC_SINE, an angle whose sine is S.
//
{
  double angle;
  double pi = 3.141592653589793;

  if ( s <= -1.0 )
  {
    angle = - pi / 2.0;
  } 
  else if ( 1.0 <= s )
  {
    angle = pi / 2.0;
  }
  else
  {
    angle = asin ( s );
  }
  return angle;
}
//****************************************************************************80

double atan4 ( double y, double x )

//****************************************************************************80
//
//  Purpose:
//
//    ATAN4 computes the inverse tangent of the ratio Y / X.
//
//  Discussion:
//
//    ATAN4 returns an angle whose tangent is ( Y / X ), a job which
//    the built in functions ATAN and ATAN2 already do.
//
//    However:
//
//    * ATAN4 always returns a positive angle, between 0 and 2 PI,
//      while ATAN and ATAN2 return angles in the interval [-PI/2,+PI/2]
//      and [-PI,+PI] respectively;
//
//    * ATAN4 accounts for the signs of X and Y, (as does ATAN2).  The ATAN
//     function by contrast always returns an angle in the first or fourth
//     quadrants.
//
//  Modified:
//
//    13 June 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double Y, X, two quantities which represent the tangent of
//    an angle.  If Y is not zero, then the tangent is (Y/X).
//
//    Output, double ATAN4, an angle between 0 and 2 * PI, whose tangent is
//    (Y/X), and which lies in the appropriate quadrant so that the signs
//    of its cosine and sine match those of X and Y.
//
{
  double pi = 3.141592653589793;
//
//  Special cases:
//
  if ( x == 0.0 )
  {
    if ( 0.0 < y )
    {
      return ( pi / 2.0 );
    } 
    else if ( y < 0.0 )
    {
      return ( 3.0 * pi / 2.0 );
    } 
    else if ( y == 0.0 )
    {
      return ( 0.0 );
    }
  } 
  else if ( y == 0.0 )
  {
    if ( 0.0 < x )
    {
      return 0.0;
    } 
    else if ( x < 0.0 )
    {
      return pi;
    }
  }
//
//  We assume that ATAN2 is reliable when both arguments are positive.
//
  if        ( 0.0 < x && 0.0 < y )
  {
    return                  atan2 (  y,  x );
  }
  else if ( x < 0.0 && 0.0 < y )
  {
    return (           pi - atan2 (  y, -x ) );
  }
  else if ( x < 0.0 && y < 0.0 )
  {
    return (           pi + atan2 ( -y, -x ) );
  }
  else if ( 0.0 < x && y < 0.0 )
  {
    return ( 2.0 * pi - atan2 ( -y,  x ) );
  }

  return 0.0;
}
//****************************************************************************80

double *ball_unit_sample_2d ( int *seed )

//****************************************************************************80
//
//  Purpose:
//
//    BALL_UNIT_SAMPLE_2D picks a random point in the unit ball in 2D.
//
//  Discussion:
//
//    The unit ball is the set of points such that
//
//      X * X + Y * Y <= 1.
//
//  Modified:
//
//    25 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input/output, int *SEED, a seed for the random number generator.
//
//    Output, double BALL_UNIT_SAMPLE_2D[2], a random point in the unit ball.
//
{
# define DIM_NUM 2

  double pi = 3.141592653589793;
  double r;
  double theta;
  double *x;

  r = r8_uniform_01 ( seed );
  r = sqrt ( r );

  theta = r8_uniform_01 ( seed );
  theta = 2.0 * pi * theta;

  x = new double[DIM_NUM];

  x[0] = r * cos ( theta );
  x[1] = r * sin ( theta );

  return x;
# undef DIM_NUM
}
//****************************************************************************80

double *ball_unit_sample_3d ( int *seed )

//****************************************************************************80
//
//  Purpose:
//
//    BALL_UNIT_SAMPLE_3D picks a random point in the unit ball in 3D.
//
//  Modified:
//
//    25 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input/output, int *SEED, a seed for the random number generator.
//
//    Output, double BALL_UNIT_SAMPLE_3D[3], the sample point.
//
{
# define DIM_NUM 3

  double phi;
  double pi = 3.141592653589793;
  double r;
  double theta;
  double vdot;
  double *x;
//
//  Pick a uniformly random VDOT, which must be between -1 and 1.
//  This represents the dot product of the random vector with the Z unit vector.
//
//  This works because the surface area of the sphere between
//  Z and Z + dZ is independent of Z.  So choosing Z uniformly chooses
//  a patch of area uniformly.
//
  vdot = r8_uniform_01 ( seed );
  vdot = 2.0 * vdot - 1.0;

  phi = arc_cosine ( vdot );
//
//  Pick a uniformly random rotation between 0 and 2 Pi around the
//  axis of the Z vector.
//
  theta = r8_uniform_01 ( seed );
  theta = 2.0 * pi * theta;
//
//  Pick a random radius R.
//
  r = r8_uniform_01 ( seed );
  r = pow ( ( double ) r, (double ) ( 1.0 / 3.0 ) );

  x = new double[DIM_NUM];

  x[0] = r * cos ( theta ) * sin ( phi );
  x[1] = r * sin ( theta ) * sin ( phi );
  x[2] = r * cos ( phi );

  return x;
# undef DIM_NUM
}
//****************************************************************************80

double *ball_unit_sample_nd ( int dim_num, int *seed )

//****************************************************************************80
//
//  Purpose:
//
//    BALL_UNIT_SAMPLE_ND picks a random point in the unit ball in ND.
//
//  Discussion:
//
//    DIM_NUM-1 random Givens rotations are applied to the point ( 1, 0, 0, ..., 0 ).
//
//    The I-th Givens rotation is in the plane of coordinate axes I and I+1,
//    and has the form:
//
//     [ cos ( theta )  - sin ( theta ) ] * x(i)      = x'(i)
//     [ sin ( theta )    cos ( theta ) ]   x(i+1)      x'(i+1)
//
//    Finally, a scaling is applied to set the point at a distance R
//    from the center, in a way that results in a uniform distribution.
//
//  Modified:
//
//    25 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the space.
//
//    Input/output, int *SEED, a seed for the random number generator.
//
//    Output, double BALL_UNIT_SAMPLE_ND[DIM_NUM], the random point.
//
{
  int i;
  double r;
  double random_cosine;
  double random_sign;
  double random_sine;
  double *x;
  double xi;

  x = new double[dim_num];

  x[0] = 1.0;
  for ( i = 1; i < dim_num; i++ )
  {
    x[i] = 0.0;
  }

  for ( i = 0; i < dim_num-1; i++ )
  {
    random_cosine = r8_uniform_01 ( seed );
    random_cosine = 2.0 * random_cosine - 1.0;
    random_sign = r8_uniform_01 ( seed );
    random_sign = ( double ) 
      ( 2 * ( int ) ( 2.0 * random_sign - 1.0 ) );
    random_sine = random_sign 
      * sqrt ( 1.0 - random_cosine * random_cosine );
    xi = x[i];
    x[i  ] = random_cosine * xi;
    x[i+1] = random_sine   * xi;
  }

  r = r8_uniform_01 ( seed );

  r = pow ( ( double ) r, 1.0 / ( double ) dim_num );

  for ( i = 0; i < dim_num; i++ )
  {
    x[i] = r * x[i];
  }

  return x;
}
//****************************************************************************80

double *basis_map_3d ( double u[3*3], double v[3*3] )

//****************************************************************************80
//
//  Purpose:
//
//    BASIS_MAP_3D computes the matrix which maps one basis to another.
//
//  Discussion:
//
//    As long as the vectors U1, U2 and U3 are linearly independent,
//    a matrix A will be computed that maps U1 to V1, U2 to V2, and
//    U3 to V3.
//
//    Depending on the values of the vectors, A may represent a
//    rotation, reflection, dilation, project, or a combination of these
//    basic linear transformations.
//
//  Modified:
//
//    29 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double U[3*3], the matrix whose columns are three "domain" or "preimage"
//    vectors, which should be linearly independent.
//
//    Input, double V[3*3], the matrix whose columns are the three "range" or 
//    "image" vectors.
//
//    Output, double BASIS_MAP_3D[3*3], a matrix with the property that A * U1 = V1,
//    A * U2 = V2 and A * U3 = V3.
//
{
  double *a;
  double *c;
  int i;
  int j;
  int k;
//
//  Compute C = the inverse of U.
//
  c = r8mat_inverse_3d ( u );

  if ( c == NULL )
  {
    return NULL;
  }
//
//  A = V * inverse ( U ).
//
  a = new double[3*3];

  for ( j = 0; j < 3; j++ )
  {
    for ( i = 0; i < 3; i++ )
    {
      a[i+j*3] = 0.0;
      for ( k = 0; k < 3; k++ )
      {
        a[i+j*3] = a[i+j*3] + v[i+k*3] * c[k+j*3];
      }
    }
  }

  delete [] c;

  return a;
}
//****************************************************************************80

bool box_01_contains_point_2d ( double p[2] )

//****************************************************************************80
//
//  Purpose:
//
//    BOX_01_CONTAINS_POINT_2D reports if a point is contained in the unit box in 2D.
//
//  Discussion:
//
//    A unit box in 2D is a rectangle with sides aligned on coordinate
//    axes.  It can be described as the set of points P satisfying:
//
//      0 <= P(1:2) <= 1.
//
//  Modified:
//
//    16 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P[2], the coordinates of the point.
//
//    Output, bool BOX_CONTAINS_POINT_2D, is TRUE if the box contains
//    the point.
//
{
  if ( 0.0 <= p[0] && p[0] <= 1.0 && 0.0 <= p[1] && p[1] <= 1.0 )
  {
    return true;
  }
  else
  {
    return false;
  }
}
//****************************************************************************80

bool box_01_contains_point_nd ( int dim_num, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    BOX_01_CONTAINS_POINT_ND reports if a point is contained in the unit box in ND.
//
//  Discussion:
//
//    A unit box is assumed to be a rectangle with sides aligned on coordinate
//    axes.  It can be described as the set of points P satisfying:
//
//      0.0 <= P(1:DIM_NUM) <= 1.0
//
//  Modified:
//
//    16 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, double P[DIM_NUM], the coordinates of the point.
//
//    Output, bool BOX_CONTAINS_POINT_ND, is TRUE if the box contains
//    the point.
//
{
  int i;

  for ( i = 0; i < dim_num; i++ )
  {
    if ( p[i] < 0.0 || 1.0 < p[i] )
    {
      return false;
    }
  }

  return true;
}
//****************************************************************************80

bool box_contains_point_2d ( double p1[2], double p2[2], double p[2] )

//****************************************************************************80
//
//  Purpose:
//
//    BOX_CONTAINS_POINT_2D reports if a point is contained in a box in 2D.
//
//  Discussion:
//
//    A box in 2D is a rectangle with sides aligned on coordinate
//    axes.  It can be described by its low and high corners, P1 and P2
//    as the set of points P satisfying:
//
//      P1(1:2) <= P(1:2) <= P2(1:2).
//
//  Modified:
//
//    16 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], the minimum and maximum X and Y
//    values, which define the box.
//
//    Input, double P[2], the coordinates of the point.
//
//    Output, bool BOX_CONTAINS_POINT_2D, is TRUE if the box contains
//    the point.
//
{
  if ( p1[0] <= p[0] && p[0] <= p2[0] && p1[1] <= p[1] && p[1] <= p2[1] )
  {
    return true;
  }
  else
  {
    return false;
  }
}
//****************************************************************************80

bool box_contains_point_nd ( int dim_num, double p1[], double p2[], double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    BOX_CONTAINS_POINT_ND reports if a point is contained in a box in ND.
//
//  Discussion:
//
//    A box in ND is a rectangle with sides aligned on coordinate
//    axes.  It can be described by its low and high corners, P1 and P2
//    as the set of points P satisfying:
//
//      P1(1:DIM_NUM) <= P(1:DIM_NUM) <= P2(1:DIM_NUM).
//
//  Modified:
//
//    16 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, double P1[DIM_NUM], P2[DIM_NUM], the minimum and maximum X and Y
//    values, which define the box.
//
//    Input, double P[DIM_NUM], the coordinates of the point.
//
//    Output, bool BOX_CONTAINS_POINT_ND, is TRUE if the box contains
//    the point.
//
{
  int i;

  for ( i = 0; i < dim_num; i++ )
  {
    if ( p[i] < p1[i] || p2[i] < p[i] )
    {
      return false;
    }
  }

  return true;
}
//****************************************************************************80

void box_ray_int_2d ( double p1[2], double p2[2], double pa[2], 
  double pb[2], double pint[2] )

//****************************************************************************80
//
//  Purpose:
//
//    BOX_RAY_INT_2D: intersection ( box, ray ) in 2D.
//
//  Discussion:
//
//    A box in 2D is a rectangle with sides aligned on coordinate
//    axes.  It can be described by its low and high corners, P1 and P2
//    as the set of points P satisfying:
//
//      P1(1:2) <= P(1:2) <= P2(1:2).
//
//    The origin of the ray is assumed to be inside the box.  This
//    guarantees that the ray will intersect the box in exactly one point.
//
//  Modified:
//
//    30 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], the lower left corner of the box.
//
//    Input, double P2[2], the upper right corner of the box.
//
//    Input, double PA[2], the origin of the ray, which should be
//    inside the box.
//
//    Input, double PB[2], a second point on the ray.
//
//    Output, double PINT[2], the point on the box intersected by the ray.
//
{
# define DIM_NUM 2

  bool inside;
  int ival;
  double pc[DIM_NUM];
  double pd[DIM_NUM];
  int side;

  for ( side = 1; side <= 4; side++ )
  {
    if ( side == 1 )
    {
      pc[0] = p1[0];
      pc[1] = p1[1];
      pd[0] = p2[0];
      pd[1] = p1[1];
    }
    else if ( side == 2 )
    {
      pc[0] = p2[0];
      pc[1] = p1[1];
      pd[0] = p2[0];
      pd[1] = p2[1];
    }
    else if ( side == 3 )
    {
      pc[0] = p2[0];
      pc[1] = p2[1];
      pd[0] = p1[0];
      pd[1] = p2[1];
    }
    else if ( side == 4 )
    {
      pc[0] = p1[0];
      pc[1] = p2[1];
      pd[0] = p1[0];
      pd[1] = p1[1];
    }

    inside = angle_contains_ray_2d ( pc, pa, pd, pb );

    if ( inside )
    {
      break;
    }

    if ( side == 4 )
    {
      cout << "\n";
      cout << "BOX_RAY_INT_2D - Fatal error!\n";
      cout << "  No intersection could be found.\n";
      exit ( 1 );
    }

  }

  lines_exp_int_2d ( pa, pb, pc, pd, &ival, pint );

  return;
# undef DIM_NUM
}
//****************************************************************************80

int box_segment_clip_2d ( double p1[2], double p2[2], double pa[2], 
  double pb[2] )

//****************************************************************************80
//
//  Purpose:
//
//    BOX_SEGMENT_CLIP_2D uses a box to clip a line segment in 2D.
//
//  Discussion:
//
//    A box in 2D is a rectangle with sides aligned on coordinate
//    axes.  It can be described by its low and high corners, P1 and P2
//    as the set of points P satisfying:
//
//      P1(1:2) <= P(1:2) <= P2(1:2).
//
//    A line segment is the finite portion of a line that lies between
//    two points.
//
//  Modified:
//
//    16 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], the minimum and maximum X and Y
//    values, which define the box.
//
//    Input/output, double PA[2], PB[2]; on input, the endpoints 
//    of a line segment.  On output, the endpoints of the portion of the
//    line segment that lies inside the box.  However, if no part of the
//    initial line segment lies inside the box, the output value is the
//    same as the input value.
//
//    Output, int BOX_SEGMENT_CLIP_LINE_2D:
//    -1, no part of the line segment is within the box.
//     0, no clipping was necessary. 
//     1, P1 was clipped.
//     2, P2 was clipped.
//     3, P1 and P2 were clipped.
//
{
  bool clip_a;
  bool clip_b;
  int ival;
  double x;
  double y;

  clip_a = false;
  clip_b = false;
//
//  Require that XMIN <= X.
//
  if ( pa[0] < p1[0] && pb[0] < p1[0] )
  {
    ival = -1;
    return ival;
  }

  if ( pa[0] < p1[0] && p1[0] <= pb[0] )
  {
    x = p1[0];
    y = pa[1] + ( pb[1] - pa[1] ) * ( x - pa[0] ) / ( pb[0] - pa[0] );
    pa[0] = x;
    pa[1] = y;
    clip_a = true;
  }
  else if ( p1[0] <= pa[0] && pb[0] < p1[0] )
  {
    x = p1[0];
    y = pa[1] + ( pb[1] - pa[1] ) * ( x - pa[0] ) / ( pb[0] - pa[0] );
    pb[0] = x;
    pb[1] = y;
    clip_b = true;
  }
//
//  Require that X <= XMAX.
//
  if ( p2[0] < pa[0] && p2[0] < pb[0] )
  {
    ival = -1;
    return ival;
  }

  if ( p2[0] < pa[0] && pb[0] <= p2[0] )
  {
    x = p2[0];
    y = pa[1] + ( pb[1] - pa[1] ) * ( x - pa[0] ) / ( pb[0] - pa[0] );
    pa[0] = x;
    pa[1] = y;
    clip_a = true;
  }
  else if ( pa[0] <= p2[0] && p2[0] < pb[0] )
  {
    x = p2[0];
    y = pa[1] + ( pb[1] - pa[1] ) * ( x - pa[0] ) / ( pb[0] - pa[0] );
    pb[0] = x;
    pb[1] = y;
    clip_b = true;
  }
//
//  Require that YMIN <= Y.
//
  if ( pa[1] < p1[1] && pb[1] < p1[1] )
  {
    ival = -1;
    return ival;
  }

  if ( pa[1] < p1[1] && p1[1] <= pb[1] )
  {
    y = p1[1];
    x = pa[0] + ( pb[0] - pa[0] ) * ( y - pa[1] ) / ( pb[1] - pa[1] );

    pa[0] = x;
    pa[1] = y;
    clip_a = true;
  }
  else if ( p1[1] <= pa[1] && pb[1] < p1[1] )
  {
    y = p1[1];
    x = pa[0] + ( pb[0] - pa[0] ) * ( y - pa[1] ) / ( pb[1] - pa[1] );
    pb[0] = x;
    pb[1] = y;
    clip_b = true;
  }
//
//  Require that Y <= YMAX.
//
  if ( p2[1] < pa[1] && p2[1] < pb[1] )
  {
    ival = -1;
    return ival;
  }

  if ( p2[1] < pa[1] && pb[1] <= p2[1] )
  {
    y = p2[1];
    x = pa[0] + ( pb[0] - pa[0] ) * ( y - pa[1] ) / ( pb[1] - pa[1] );
    pa[0] = x;
    pa[1] = y;
    clip_a = true;
  }
  else if ( pa[1] <= p2[1] && p2[1] < pb[1] )
  {
    y = p2[1];
    x = pa[0] + ( pb[0] - pa[0] ) * ( y - pa[1] ) / ( pb[1] - pa[1] );
    pb[0] = x;
    pb[1] = y;
    clip_b = true;
  }

  ival = 0;

  if ( clip_a )
  {
    ival = ival + 1;
  }

  if ( clip_b)
  {
    ival = ival + 2;
  }

  return ival;
}
//****************************************************************************80

char ch_cap ( char c )

//****************************************************************************80
//
//  Purpose:
//
//    CH_CAP capitalizes a single character.
//
//  Discussion:
//
//    This routine should be equivalent to the library "toupper" function.
//
//  Modified:
//
//    19 July 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char C, the character to capitalize.
//
//    Output, char CH_CAP, the capitalized character.
//
{
  if ( 97 <= c && c <= 122 ) 
  {
    c = c - 32;
  }   

  return c;
}
//****************************************************************************80

bool ch_eqi ( char c1, char c2 )

//****************************************************************************80
//
//  Purpose:
//
//    CH_EQI is true if two characters are equal, disregarding case.
//
//  Modified:
//
//    13 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char C1, char C2, the characters to compare.
//
//    Output, bool CH_EQI, is true if the two characters are equal,
//    disregarding case.
//
{
  if ( 97 <= c1 && c1 <= 122 ) 
  {
    c1 = c1 - 32;
  } 
  if ( 97 <= c2 && c2 <= 122 ) 
  {
    c2 = c2 - 32;
  }     

  return ( c1 == c2 );
}
//****************************************************************************80

void circle_arc_point_near_2d ( double r, double pc[2], double theta1, 
  double theta2, double p[2], double pn[2], double *dist )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_ARC_POINT_NEAR_2D : nearest point on a circular arc.
//
//  Discussion:
//
//    A circular arc is defined by the portion of a circle (R,PC)
//    between two angles (THETA1,THETA2).
//
//    A point P on a circular arc satisfies
//
//      ( X - PC(1) ) * ( X - PC(1) ) + ( Y - PC(2) ) * ( Y - PC(2) ) = R * R
//
//    and
//
//      Theta1 <= Theta <= Theta2
//
//  Modified:
//
//    28 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the center of the circle.
//
//    Input, double THETA1, THETA2, the angles defining the arc,
//    in radians.  Normally, THETA1 < THETA2.
//
//    Input, double P[2], the point to be checked.
//
//    Output, double PN[2], a point on the circular arc which is
//    nearest to the point.
//
//    Output, double *DIST, the distance to the nearest point.
//
{
# define DIM_NUM 2

  int i;
  double pi = 3.141592653589793;
  double r2;
  double theta;
//
//  Special case, the zero circle.
//
  if ( r == 0.0 )
  {
    r8vec_copy ( DIM_NUM, pc, pn );

    *dist = 0.0;
    for ( i = 0; i < DIM_NUM; i++ )
    {
      *dist = *dist + pow ( p[i] - pn[i], 2 );
    }
    *dist = sqrt ( *dist );
    return;
  }
//
//  Determine the angle made by the point.
//
  theta = atan4 ( p[1] - pc[1], p[0] - pc[0] );
//
//  If the angle is between THETA1 and THETA2, then you can
//  simply project the point onto the arc.
//
  if ( r8_modp ( theta  - theta1,  2.0 * pi ) <= 
       r8_modp ( theta2 - theta1,  2.0 * pi ) )
  {
    r2 = 0.0;
    for ( i = 0; i < DIM_NUM; i++ )
    {
      r2 = r2 + pow ( p[i] - pc[i], 2 );
    }
    r2 = sqrt ( r2 );
    for ( i = 0; i < DIM_NUM; i++ )
    {
      pn[i] = pc[i] + ( p[i] - pc[i] ) * r / r2;
    }
  }
//
//  Otherwise, if the angle is less than the negative of the
//  average of THETA1 and THETA2, it's on the side of the arc
//  where the endpoint associated with THETA2 is closest.
//
  else if ( r8_modp ( theta - 0.5 * ( theta1 + theta2 ), 2.0 * pi ) 
    <= pi )
  {
    pn[0] = pc[0] + r * cos ( theta2 );
    pn[1] = pc[1] + r * sin ( theta2 );
  }
//
//  Otherwise, the endpoint associated with THETA1 is closest.
//
  else
  {
    pn[0] = pc[0] + r * cos ( theta1 );
    pn[1] = pc[1] + r * sin ( theta1 );
  }

  *dist = 0.0;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    *dist = *dist + pow ( p[i] - pn[i], 2 );
  }
  *dist = sqrt ( *dist );

  return;
# undef DIM_NUM
}
//****************************************************************************80

double circle_area_2d ( double r )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_AREA_2D computes the area of a circle in 2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Output, double CIRCLE_AREA_2D, the area of the circle.
//
{
  double area;
  double pi = 3.141592653589793;

  area = pi * r * r;

  return area;
}
//****************************************************************************80

void circle_dia2imp_2d ( double p1[2], double p2[2], double *r, double pc[2] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_DIA2IMP_2D converts a diameter to an implicit circle in 2D.
//
//  Discussion:
//
//    The diameter form of a circle is:
//
//      P1 and P2 are endpoints of a diameter.
//
//    The implicit form of a circle in 2D is:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//  Modified:
//
//    21 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], are the X and Y coordinates
//    of two points which form a diameter of the circle.
//
//    Output, double *R, the computed radius of the circle.
//
//    Output, double PC[2], the computed center of the circle.
//
{
  *r = 0.5 * sqrt ( pow ( p1[0] - p2[0], 2 )
                  + pow ( p1[1] - p2[1], 2 ) );

  pc[0] = 0.5 * ( p1[0] + p2[0] );
  pc[1] = 0.5 * ( p1[1] + p2[1] );

  return;
}
//****************************************************************************80

int circle_exp_contains_point_2d ( double p1[2], double p2[2], double p3[2], 
  double p[2] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_EXP_CONTAINS_POINT_2D determines if an explicit circle contains a point in 2D.
//
//  Discussion:
//
//    The explicit form of a circle in 2D is:
//
//      The circle passing through points P1, P2 and P3.
//
//  Modified:
//
//    28 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], P3[2], the coordinates of three
//    points that lie on a circle.
//
//    Input, double P[2], the coordinates of a point, whose position
//    relative to the circle is desired.
//
//    Output, int CIRCLE_EXP_CONTAINS_POINT_2D:
//   -1, the three points are distinct and noncolinear,
//      and the point lies inside the circle.
//    0, the three points are distinct and noncolinear,
//      and the point lies on the circle.
//    1, the three points are distinct and noncolinear,
//      and the point lies outside the circle.
//    2, the three points are distinct and colinear,
//      and the point lies on the line.
//    3, the three points are distinct and colinear,
//      and the point does not lie on the line.
//    4, two points are distinct, and the point lies on the line.
//    5, two points are distinct, and the point does not lie on the line.
//    6, all three points are equal, and the point is equal to them,
//    7, all three points are equal, and the point is not equal to them.
//
{
  double a[4*4];
  double det;
  int inside;
//
//  P1 = P2?
//
  if ( r8vec_eq ( 2, p1, p2 ) )
  {
    if ( r8vec_eq ( 2, p1, p3 ) )
    {
      if ( r8vec_eq ( 2, p1, p ) )
      {
        inside = 6;
      }
      else
      {
        inside = 7;
      }
    }
    else
    {

      det = ( p1[0] - p3[0] ) * ( p[1]  - p3[1] ) 
          - ( p[0]  - p3[0] ) * ( p1[1] - p3[1] );

      if ( det == 0.0 )
      {
        inside = 4;
      }
      else
      {
        inside = 5;
      }
    }
    return inside;
  }
//
//  P1 does not equal P2.  Does P1 = P3?
//
  if ( r8vec_eq ( 2, p1, p3 ) )
  {
    det = ( p1[0] - p2[0] ) * ( p[1] - p2[1] )
        - ( p[0] - p2[0] ) * ( p1[1] - p2[1] );

    if ( det == 0.0 )
    {
      inside = 4;
    }
    else
    {
      inside = 5;
    }
    return inside;
  }
//
//  The points are distinct.  Are they colinear?
//
  det = ( p1[0] - p2[0] ) * ( p3[1] - p2[1] ) 
      - ( p3[0] - p2[0] ) * ( p1[1] - p2[1] );

  if ( det == 0.0 )
  {
    det = ( p1[0] - p2[0] ) * ( p[1] - p2[1] ) 
        - ( p[0] - p2[0] ) * ( p1[1] - p2[1] );

    if ( det == 0.0 )
    {
      inside = 2;
    }
    else
    {
      inside = 3;
    }

    return inside;

  }
//
//  The points are distinct and non-colinear.
//
//  Compute the determinant
//
  a[0+0*4] = p1[0];
  a[1+0*4] = p2[0];
  a[2+0*4] = p3[0];
  a[3+0*4] = p[0];

  a[0+1*4] = p1[1];
  a[1+1*4] = p2[1];
  a[2+1*4] = p3[1];
  a[3+1*4] = p[1];

  a[0+2*4] = p1[0] * p1[0] + p1[1] * p1[1];
  a[1+2*4] = p2[0] * p2[0] + p2[1] * p2[1];
  a[2+2*4] = p3[0] * p3[0] + p3[1] * p3[1];
  a[3+2*4] = p[0]  * p[0]  + p[1]  * p[1];

  a[0+3*4] = 1.0;
  a[1+3*4] = 1.0;
  a[2+3*4] = 1.0;
  a[3+3*4] = 1.0;

  det = r8mat_det_4d ( a );

  if ( det < 0.0 )
  {
    inside = 1;
  }
  else if ( det == 0.0 )
  {
    inside = 0;
  }
  else
  {
    inside = -1;
  }

  return inside;
}
//****************************************************************************80

void circle_exp2imp_2d ( double p1[2], double p2[2], double p3[2], double *r, 
  double pc[2] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_EXP2IMP_2D converts a circle from explicit to implicit form in 2D.
//
//  Discussion:
//
//    The explicit form of a circle in 2D is:
//
//      The circle passing through points P1, P2 and P3.
//
//    Points P on an implicit circle in 2D satisfy the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//    Any three points define a circle, as long as they don't lie on a straight
//    line.  (If the points do lie on a straight line, we could stretch the
//    definition of a circle to allow an infinite radius and a center at
//    some infinite point.)
//
//    Instead of the formulas used here, you can use the linear system
//    approach in the routine TRIANGLE_OUTCIRCLE_2D.
//
//    The diameter of the circle can be found by solving a 2 by 2 linear system.
//    This is because the vectors P2 - P1 and P3 - P1 are secants of the circle,
//    and each forms a right triangle with the diameter.  Hence, the dot product
//    of P2 - P1 with the diameter is equal to the square of the length
//    of P2 - P1, and similarly for P3 - P1.  These two equations determine the
//    diameter vector originating at P1.
//
//    If all three points are equal, return a circle of radius 0 and
//    the obvious center.
//
//    If two points are equal, return a circle of radius half the distance
//    between the two distinct points, and center their average.
//
//  Modified:
//
//    10 March 2006
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Joseph ORourke,
//    Computational Geometry,
//    Second Edition,
//    Cambridge, 1998,
//    ISBN: 0521649765,
//    LC: QA448.D38.
//
//  Parameters:
//
//    Input, double P1[2], P2[2], P3[2], are the coordinates
//    of three points that lie on the circle.  These points should be
//    distinct, and not collinear.
//
//    Output, double *R, the radius of the circle.  Normally, R will be positive.
//    R will be (meaningfully) zero if all three points are
//    equal.  If two points are equal, R is returned as the distance between
//    two nonequal points.  R is returned as -1 in the unlikely event that
//    the points are numerically collinear; philosophically speaking, R
//    should actually be "infinity" in this case.
//
//    Output, double PC[2], the center of the circle.
//
{
# define DIM_NUM 2

  double a;
  double b;
  double c;
  double d;
  double e;
  double f;
  double g;
//
//  If all three points are equal, then the
//  circle of radius 0 and center P1 passes through the points.
//
  if ( r8vec_eq ( DIM_NUM, p1, p2 ) && r8vec_eq ( DIM_NUM, p1, p3 ) )
  {
    *r = 0.0;
    r8vec_copy ( DIM_NUM, p1, pc );
    return;
  }
//
//  If exactly two points are equal, then the circle is defined as
//  having the obvious radius and center.
//
  if ( r8vec_eq ( DIM_NUM, p1, p2 ) )
  {
    *r = 0.5 * sqrt ( ( p1[0] - p3[0] ) * ( p1[0] - p3[0] ) 
                    + ( p1[1] - p3[1] ) * ( p1[1] - p3[1] ) );
    pc[0] = 0.5 * ( p1[0] + p3[0] );
    pc[1] = 0.5 * ( p1[1] + p3[1] );
    return;
  }
  else if ( r8vec_eq ( DIM_NUM, p1, p3 ) )
  {
    *r = 0.5 * sqrt ( ( p1[0] - p2[0] ) * ( p1[0] - p2[0] ) 
                    + ( p1[1] - p2[1] ) * ( p1[1] - p2[1] ) );
    pc[0] = 0.5 * ( p1[0] + p2[0] );
    pc[1] = 0.5 * ( p1[1] + p2[1] );
    return;
  }
  else if ( r8vec_eq ( DIM_NUM, p2, p3 ) )
  {
    *r = 0.5 * sqrt ( ( p1[0] - p2[0] ) * ( p1[0] - p2[0] ) 
                    + ( p1[1] - p2[1] ) * ( p1[1] - p2[1] ) );
    pc[0] = 0.5 * ( p1[0] + p2[0] );
    pc[1] = 0.5 * ( p1[1] + p2[1] );
    return;
  }

  a = p2[0] - p1[0];
  b = p2[1] - p1[1];
  c = p3[0] - p1[0];
  d = p3[1] - p1[1];

  e = a * ( p1[0] + p2[0] ) + b * ( p1[1] + p2[1] );
  f = c * ( p1[0] + p3[0] ) + d * ( p1[1] + p3[1] );
//
//  Our formula is:
//
//    G = a * ( d - b ) - b * ( c - a )
//
//  but we get slightly better results using the original data.
//
  g = a * ( p3[1] - p2[1] ) - b * ( p3[0] - p2[0] );
//
//  We check for collinearity.  A more useful check would compare the
//  absolute value of G to a small quantity.
//
  if ( g == 0.0 )
  {
    pc[0] = 0.0;
    pc[1] = 0.0;
    *r = -1.0;
    return;
  }
//
//  The center is halfway along the diameter vector from P1.
//
  pc[0] = 0.5 * ( d * e - b * f ) / g;
  pc[1] = 0.5 * ( a * f - c * e ) / g;
//
//  Knowing the center, the radius is now easy to compute.
//
  *r = sqrt ( ( p1[0] - pc[0] ) * ( p1[0] - pc[0] ) 
            + ( p1[1] - pc[1] ) * ( p1[1] - pc[1] ) );

  return;
# undef DIM_NUM
}
//****************************************************************************80

bool circle_imp_contains_point_2d ( double r, double pc[2], double p[2] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP_CONTAINS_POINT_2D determines if an implicit circle contains a point in 2D.
//
//  Discussion:
//
//    An implicit circle in 2D satisfies the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//  Modified:
//
//    28 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the coordinates of the center of the circle.
//
//    Input, double P[2], the point to be checked.
//
//    Output, bool CIRCLE_IMP_CONTAINS_POINT_2D, is true if the point 
//    is inside or on the circle, false otherwise.
//
{
  if ( pow ( p[0] - pc[0], 2 ) + pow ( p[1] - pc[1], 2 ) <= r * r )
  {
    return true;
  }
  else
  {
    return false;
  }
}
//****************************************************************************80

void circle_imp_line_par_int_2d ( double r, double pc[2], double x0, double y0, 
  double f, double g, int *int_num, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP_LINE_PAR_INT_2D: ( implicit circle, parametric line ) intersection in 2D.
//
//  Discussion:
//
//    An implicit circle in 2D satisfies the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//    The parametric form of a line in 2D is:
//
//      X = X0 + F * T
//      Y = Y0 + G * T
//
//    For normalization, we choose F*F+G*G = 1 and 0 <= F.
//
//  Modified:
//
//    29 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the coordinates of the center of the circle.
//
//    Input, double F, G, X0, Y0, the parametric parameters of the line.
//
//    Output, int *INT_NUM, the number of intersecting points found.
//    INT_NUM will be 0, 1 or 2.
//
//    Output, double P[2*2], contains the X and Y coordinates of
//    the intersecting points.
//
{
  double root;
  double t;

  root = r * r * ( f * f + g * g )
    - ( f * ( pc[1] - y0 ) - g * ( pc[0] - x0 ) ) 
    * ( f * ( pc[1] - y0 ) - g * ( pc[0] - x0 ) );

  if ( root < 0.0 )
  {
    *int_num = 0;
  }
  else if ( root == 0.0 )
  {

    *int_num = 1;

    t = ( f * ( pc[0] - x0 ) + g * ( pc[1] - y0 ) ) / ( f * f + g * g );
    p[0+0*2] = x0 + f * t;
    p[1+0*2] = y0 + g * t;

  }
  else if ( 0.0 < root )
  {

    *int_num = 2;

    t = ( ( f * ( pc[0] - x0 ) + g * ( pc[1] - y0 ) ) - sqrt ( root ) )
      / ( f * f + g * g );

    p[0+0*2] = x0 + f * t;
    p[1+0*2] = y0 + g * t;

    t = ( ( f * ( pc[0] - x0 ) + g * ( pc[1] - y0 ) ) + sqrt ( root ) )
      / ( f * f + g * g );

    p[0+1*2] = x0 + f * t;
    p[1+1*2] = y0 + g * t;
  }

  return;
}
//****************************************************************************80

double circle_imp_point_dist_2d ( double r, double pc[2], double p[2]  )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP_POINT_DIST_2D: distance ( implicit circle, point ) in 2D.
//
//  Discussion:
//
//    The distance is zero if the point is on the circle.
//
//    An implicit circle in 2D satisfies the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//  Modified:
//
//    29 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the coordinates of the center of the circle.
//
//    Input, double P[2], the point to be checked.
//
//    Output, double CIRCLE_IMP_POINT_DIST_2D, the distance of the point 
//    to the circle.
//
{
  double value;

  value = sqrt ( r8_abs ( pow ( p[0] - pc[0], 2 ) + pow ( p[1] - pc[1], 2 )
    - r * r ) );

  return value;
}
//****************************************************************************80

double circle_imp_point_dist_signed_2d ( double r, double pc[2], double p[2] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP_POINT_DIST_SIGNED_2D: signed distance ( implicit circle, point ) in 2D.
//
//  Discussion:
//
//    The signed distance is zero if the point is on the circle.
//    The signed distance is negative if the point is inside the circle.
//
//    An implicit circle in 2D satisfies the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//  Modified:
//
//    29 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the coordinates of the center of the circle.
//
//    Input, double P[2], the point to be checked.
//
//    Output, double CIRCLE_IMP_POINT_DIST_SIGNED_2D, the signed distance 
//    of the point to the circle.  If the point is inside the circle, 
//    the signed distance is negative.
//
{
  double t;
  double value;

  t = pow ( p[0] - pc[0], 2 ) + pow ( p[1] - pc[1], 2 ) - r * r;

  value = r8_sign ( t ) * sqrt ( r8_abs ( t ) );

  return value;
}
//****************************************************************************80

double circle_imp_point_near_2d ( double r, double pc[2], double p[2], double pn[2] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP_POINT_NEAR_2D: nearest ( implicit circle, point ) in 2D.
//
//  Discussion:
//
//    This routine finds the distance from a point to an implicitly
//    defined circle, and returns the point on the circle that is
//    nearest to the given point.
//
//    If the given point is the center of the circle, than any point
//    on the circle is "the" nearest.
//
//    An implicit circle in 2D satisfies the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//  Modified:
//
//    08 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the center of the circle.
//
//    Input, double P[2], the point to be checked.
//
//    Output, double PN[2], the nearest point on the circle.
//
//    Output, double CIRCLE_IMP_POINT_NEAR_2D, the distance of the point to the circle.
//
{
# define DIM_NUM 2

  double dist;
  int i;
  double r2;

  if ( r8vec_eq ( DIM_NUM, p, pc ) )
  {
    dist = r;
    for ( i = 0; i < DIM_NUM; i++ )
    {
      pn[i] = pc[i] + r / sqrt ( ( double ) ( DIM_NUM ) );
    }
    return dist;
  }

  r2 = 0.0;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    r2 = r2 + pow ( p[i] - pc[i], 2 );
  }
  r2 = sqrt ( r2 );

  dist = r8_abs (  r2 - r );

  for ( i = 0; i < DIM_NUM; i++ )
  {
    pn[i] = pc[i] + r * ( p[i] - pc[i] ) / r2;
  }

  return dist;
# undef DIM_NUM
}
//****************************************************************************80

double *circle_imp_points_2d ( double r, double pc[2], int n )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP_POINTS_2D returns N equally spaced points on an implicit circle in 2D.
//
//  Discussion:
//
//    The first point is always ( PC[0] + R, PC[1] ), and subsequent points
//    proceed counter clockwise around the circle.
//
//    An implicit circle in 2D satisfies the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//  Modified:
//
//    10 March 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the coordinates of the center of the circle.
//
//    Input, int N, the number of points desired.  N must be at least 1.
//
//    Output, double CIRCLE_IMP_POINTS_2D[2*N], points on the circle.
//
{
  double angle;
  int j;
  double *p;
  double pi = 3.141592653589793;

  p = new double[2*n];

  for ( j = 0; j < n; j++ )
  {
    angle = ( 2.0 * pi * ( double ) j ) / ( double ) n ;
    p[0+j*2] = pc[0] + r * cos ( angle );
    p[1+j*2] = pc[1] + r * sin ( angle );
  }

  return p;
}
//****************************************************************************80

double *circle_imp_points_3d ( double r, double pc[3], double nc[3], int n )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP_POINTS_3D returns points on an implicit circle in 3D.
//
//  Discussion:
//
//    Points P on an implicit circle in 3D satisfy the equations:
//
//      ( P(1) - PC(1) )**2
//    + ( P(2) - PC(2) )**2
//    + ( P(3) - PC(3) )**2 = R**2
//
//    and
//
//      ( P(1) - PC(1) ) * NC(1)
//    + ( P(2) - PC(2) ) * NC(2)
//    + ( P(3) - PC(3) ) * NC(3) = 0
//
//  Modified:
//
//    10 March 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[3], the center of the circle.
//
//    Input, double NC[3], a nonzero vector that is normal to
//    the plane of the circle.  It is customary, but not necessary,
//    that this vector have unit norm.
//
//    Input, int N, the number of points desired.
//    N must be at least 1.
//
//    Output, double CIRCLE_IMP_POINTS_3D[3*N], the coordinates of points
//    on the circle.
//
{
# define DIM_NUM 3

  int i;
  int j;
  double *n1;
  double *n2;
  double *p;
  double pi = 3.141592653589793;
  double theta;
//
//  Get two unit vectors N1 and N2 which are orthogonal to each other,
//  and to NC.
//
  n1 = new double[DIM_NUM];
  n2 = new double[DIM_NUM];

  plane_normal_basis_3d ( pc, nc, n1, n2 );
//
//  Rotate R units away from PC in the plane of N1 and N2.
//
  p = new double[DIM_NUM*n];

  for ( j = 0; j < n; j++ )
  {
    theta = ( 2.0 * pi * ( double ) ( j ) ) / ( double ) ( n );

    for ( i = 0; i < DIM_NUM; i++ )
    {
      p[i+j*DIM_NUM] = pc[i] + r * ( cos ( theta ) * n1[i] 
                                   + sin ( theta ) * n2[i] );
    }
  }

  delete [] n1;
  delete [] n2;

  return p;
# undef DIM_NUM
}
//****************************************************************************80

void circle_imp_points_arc_2d ( double r, double pc[2], double theta1, 
  double theta2, int n, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP_POINTS_ARC_2D returns N points on an arc of an implicit circle in 2D.
//
//  Discussion:
//
//    The first point is ( PC[0] + R * COS ( THETA1 ), PC[1] + R * SIN ( THETA1 ) );
//    The last point is  ( PC[0] + R * COS ( THETA2 ), PC[1] + R * SIN ( THETA2 ) );
//    and the intermediate points are evenly spaced in angle between these,
//    and in counter clockwise order.
//
//    An implicit circle in 2D satisfies the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//  Modified:
//
//    29 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the coordinates of the center of the circle.
//
//    Input, double THETA1, THETA2, the angular coordinates of the first
//    and last points to be drawn, in radians.
//
//    Input, int N, the number of points desired.  N must be at least 1.
//
//    Output, double P[2*N], the coordinates of points on the circle.
//
{
  int i;
  double pi = 3.141592653589793;
  double theta;
  double theta3;
//
//  THETA3 is the smallest angle, no less than THETA1, which
//  coincides with THETA2.
//
  theta3 = theta1 + r8_modp ( theta2 - theta1, 2.0 * pi );

  for ( i = 0; i < n; i++ )
  {
    if ( 1 < n )
    {
      theta = ( ( double ) ( n - i - 1 ) * theta1 
              + ( double ) (     i     ) * theta3 ) 
              / ( double ) ( n     - 1 );
    }
    else
    {
      theta = 0.5 * ( theta1 + theta3 );
    }

    p[0+i*2] = pc[0] + r * cos ( theta );
    p[1+i*2] = pc[1] + r * sin ( theta );
  }

  return;
}
//****************************************************************************80

void circle_imp_print_2d ( double r, double pc[2], char *title )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP_PRINT_2D prints an implicit circle in 2D.
//
//  Discussion:
//
//    An implicit circle in 2D satisfies the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//  Modified:
//
//    26 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the center of the circle.
//
//    Input, char *TITLE, an optional title.
//
{
# define DIM_NUM 2

  if ( 0 < s_len_trim ( title ) )
  {
    cout << "\n";
    cout << title << "\n";
  }

  cout << "\n";
  cout << "  Radius = " << r << "\n";
  cout << "  Center = (" << pc[0] << ",  " << pc[1] << ")\n";

  return;
# undef DIM_NUM
}
//****************************************************************************80

void circle_imp_print_3d ( double r, double pc[3], double nc[3], char *title )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP_PRINT_2D prints an implicit circle in 3D.
//
//  Discussion:
//
//    Points P on an implicit circle in 3D satisfy the equations:
//
//      ( P(1) - PC(1) )**2
//    + ( P(2) - PC(2) )**2
//    + ( P(3) - PC(3) )**2 = R**2
//
//    and
//
//      ( P(1) - PC(1) ) * NC(1)
//    + ( P(2) - PC(2) ) * NC(2)
//    + ( P(3) - PC(3) ) * NC(3) = 0
//
//  Modified:
//
//    10 March 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[3], the center of the circle.
//
//    Input, double NC[3], the normal vector to the circle.
//
//    Input, character *TITLE, an optional title.
//
{
# define DIM_NUM 3

  if ( 0 < s_len_trim ( title ) )
  {
    cout << "\n";
    cout << title << "\n";
  }

  cout << "\n";
  cout << "  Radius = " << r << "\n";
  cout << "  Center = (" << pc[0] 
                 << ", " << pc[1] 
                 << ", " << pc[2] << ")\n";
  cout << "  Normal = (" << nc[0] 
                 << ", " << nc[1] 
                 << ", " << nc[2] << ")\n";

  return;
# undef DIM_NUM
}
//****************************************************************************80

void circle_imp2exp_2d ( double r, double pc[2], double p1[2], double p2[2], 
  double p3[2] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_IMP2EXP_2D converts a circle from implicit to explicit form in 2D.
//
//  Discussion:
//
//    Points P on an implicit circle in 2D satisfy the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//    The explicit form of a circle in 2D is:
//
//      The circle passing through points P1, P2 and P3.
//
//  Modified:
//
//    17 May 2007
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Joseph ORourke,
//    Computational Geometry,
//    Second Edition,
//    Cambridge, 1998,
//    ISBN: 0521649765,
//    LC: QA448.D38.
//
//  Parameters:
//
//    Input, double R, PC[2], the radius and center of the circle.
//
//    Output, double P1[2], P2[2], P3[2], three points on the circle.
//
{
  double pi = 3.141592653589793;
  double theta;

  theta = 0.0;
  p1[0] = pc[0] + r * cos ( theta );
  p1[1] = pc[1] + r * sin ( theta );

  theta = 2.0 * pi / 3.0;
  p2[0] = pc[0] + r * cos ( theta );
  p2[1] = pc[1] + r * sin ( theta );

  theta = 4.0 * pi / 3.0;
  p3[0] = pc[0] + r * cos ( theta );
  p3[1] = pc[1] + r * sin ( theta );

  return;
}
//****************************************************************************80

double *circle_llr2imp_2d ( double p1[], double p2[], double q1[], double q2[], 
  double r )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_LLR2IMP_2D converts a circle from LLR to implicit form in 2D.
//
//  Discussion:
//
//    The LLR form of a circle in 2D is:
//
//      The circle of radius R tangent to the lines L1 and L2.
//
//    The implicit form of a circle in 2D is:
//
//      ( P(1) - PC(1) )**2 + ( P(2) - PC(2) )**2 = R**2
//
//    Let S be the scaled distance of a point on L1 from P1 to P2,
//    and let N1 be a unit normal vector to L1.  Then a point P that is
//    R units from L1 satisfies:
//
//      P = P1 + s * ( P2 - P1 ) + R * N1.
//
//    Let t be the scaled distance of a point on L2 from Q1 to Q2,
//    and let N2 be a unit normal vector to L2.  Then a point Q that is
//    R units from L2 satisfies:
//
//      Q = Q1 + t * ( Q2 - Q1 ) + R * N2.
//
//    For the center of the circle, then, we have P = Q, that is
//
//      ( P2 - P1 ) * s + ( Q2 - Q1 ) * t = - P1 - Q1 - R * ( N1 + N2 )
//
//    This is a linear system for ( s and t ) from which we can compute
//    the points of tangency, and the center.
//
//    Note that we have four choices for the circle based on the use
//    of plus or minus N1 and plus or minus N2.
//
//  Modified:
//
//    16 November 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], two points on line 1.
//
//    Input, double Q1[2], Q2[2], two points on line 2.
//
//    Input, double R, the radius of the circle.
//
//    Output, double CIRCLE_LLR2IMP_2D[2*4], the centers of the circles.
//
{
  double *a;
  double *b;
  double det;
  double *n1;
  double *n2;
  double *pc;
  double *x;

  a = new double[2*2];
  b = new double[2];
  pc = new double[2*4];
//
//  Compute the normals N1 and N2.
//
  n1 = line_exp_normal_2d ( p1, p2 );

  n2 = line_exp_normal_2d ( q1, q2 );
//
//  Set the linear system.
//
  a[0+0*2] =   p2[0] - p1[0];
  a[1+0*2] =   p2[1] - p1[1];
  a[0+1*2] = - q2[0] + q1[0];
  a[1+1*2] = - q2[1] + q1[1];
//
//  Solve the 4 linear systems, using every combination of
//  signs on the normal vectors.
//
  b[0] = - p1[0] + q1[0] + r * n1[0] + r * n2[0];
  b[1] = - p1[1] + q1[1] + r * n1[1] + r * n2[1];

  x = r8mat_solve_2d ( a, b, &det );

  pc[0+2*0] = p1[0] + ( p2[0] - p1[0] ) * x[0] - r * n1[0];
  pc[1+2*0] = p1[1] + ( p2[1] - p1[1] ) * x[0] - r * n1[1];

  delete [] x;

  b[0] = - p1[0] + q1[0] + r * n1[0] - r * n2[0];
  b[1] = - p1[1] + q1[1] + r * n1[1] - r * n2[1];

  x = r8mat_solve_2d ( a, b, &det );

  pc[0+2*1] = p1[0] + ( p2[0] - p1[0] ) * x[0] - r * n1[0];
  pc[1+2*1] = p1[1] + ( p2[1] - p1[1] ) * x[0] - r * n1[1];

  delete [] x;

  b[0] = - p1[0] + q1[0] - r * n1[0] + r * n2[0];
  b[1] = - p1[1] + q1[1] - r * n1[1] + r * n2[1];

  x = r8mat_solve_2d ( a, b, &det );

  pc[0+2*2] = p1[0] + ( p2[0] - p1[0] ) * x[0] + r * n1[0];
  pc[1+2*2] = p1[1] + ( p2[1] - p1[1] ) * x[0] + r * n1[1];

  delete [] x;

  b[0] = - p1[0] + q1[0] - r * n1[0] - r * n2[0];
  b[1] = - p1[1] + q1[1] - r * n1[1] - r * n2[1];

  x = r8mat_solve_2d ( a, b, &det );

  pc[0+2*3] = p1[0] + ( p2[0] - p1[0] ) * x[0] + r * n1[0];
  pc[1+2*3] = p1[1] + ( p2[1] - p1[1] ) * x[0] + r * n1[1];

  delete [] a;
  delete [] b;
  delete [] n1;
  delete [] n2;
  delete [] x;

  return pc;
}
//****************************************************************************80

double circle_lune_area_2d ( double r, double pc[2], double theta1, double theta2 )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_LUNE_AREA_2D returns the area of a circular lune in 2D.
//
//  Discussion:
//
//    A lune is formed by drawing a circular arc, and joining its endpoints.
//
//  Modified:
//
//    10 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the center of the circle.
//
//    Input, double THETA1, THETA2, the angles defining the arc,
//    in radians.  Normally, THETA1 < THETA2.
//
//    Output, double CIRCLE_LUNE_AREA_2D, the area of the lune.
//
{
  double area;
  double area_sector;
  double area_triangle;

  area_sector = circle_sector_area_2d ( r, pc, theta1, theta2 );
  area_triangle = circle_triangle_area_2d ( r, pc, theta1, theta2 );
  area = area_sector - area_triangle;

  return area;
}
//****************************************************************************80

double *circle_lune_centroid_2d ( double r, double pc[2], double theta1,
  double theta2 )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_LUNE_CENTROID_2D returns the centroid of a circular lune in 2D.
//
//  Discussion:
//
//    A lune is formed by drawing a circular arc, and joining its endpoints.
//
//  Modified:
//
//    30 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Adrian Bowyer, John Woodwark,
//    A Programmer's Geometry,
//    Butterworths, 1983.
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the coordinates of the center of the circle.
//
//    Input, double THETA1, THETA2, the angles of the first and last points
//    on the circular arc.
//
//    Output, double CIRCLE_LUNE_CENTROID_2D[2], the coordinates of the centroid 
//    of the lune.
//
{
# define DIM_NUM 2

  double *centroid;
  double d;
  double theta;

  theta = theta2 - theta1;

  if ( theta == 0.0 )
  {
    d = r;
  }
  else
  {
    d = 4.0 * r * pow ( ( sin ( 0.5 * theta ) ), 3 ) /
      ( 3.0 * ( theta - sin ( theta ) ) );
  }

  centroid = new double[DIM_NUM];

  centroid[0] = pc[0] + d * cos ( theta );
  centroid[1] = pc[1] + d * sin ( theta );

  return centroid;
# undef DIM_NUM
}
//****************************************************************************80

void circle_pppr2imp_3d ( double p1[], double p2[], double p3[], double r, 
  double pc[], double normal[] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_PPR2IMP_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.
//
//  Modified:
//
//    12 November 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[3], P2[3], two points on the circle.
//
//    Input, double P3[3], a third point.
//
//    Input, double R, the radius of the circle.
//
//    Output, double PC[3*2], the centers of the two circles.
//
//    Output, double NORMAL[3], the normal to the circles.
//
{
# define DIM_NUM 3

  double dist;
  double dot;
  double h;
  int i;
  int j;
  double length;
  double *v;
//
//  Compute the distance from P1 to P2.
//
  dist = r8vec_distance ( DIM_NUM, p1, p2 );
//
//  If R is smaller than DIST, we don't have a circle.
//
  if ( 2.0 * r < dist )
  {
    for ( j = 0; j < 2; j++ )
    {
      for ( i = 0; i < DIM_NUM; i++ )
      {
        pc[i+j*DIM_NUM] = 0.5 * ( p1[i] + p2[i] );
      }
    }
  }
//
//  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 = new double[DIM_NUM];

  for ( i = 0; i < DIM_NUM; i++ )
  {
    v[i] = p3[i] - p1[i];
  }
  dot = 0.0;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    dot = dot + v[i] * ( p2[i] - p1[i] );
  }
  dot = dot / dist;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    v[i] = v[i] - dot * ( p2[i] - p1[i] ) / dist;
  }

  length = r8vec_length ( DIM_NUM, v );
  for ( i = 0; i < DIM_NUM; i++ )
  {
    v[i] = v[i] / length;
  }
//
//  We can go with or against the given normal direction.
//
  for ( i = 0; i < DIM_NUM; i++ )
  {
    pc[i+0*DIM_NUM] = 0.5 * ( p2[i] + p1[i] ) + h * v[i];
    pc[i+1*DIM_NUM] = 0.5 * ( p2[i] + p1[i] ) - h * v[i];
  }
  delete [] v;

  plane_exp_normal_3d ( p1, p2, p3, normal );

  return;
# undef DIM_NUM
}
//****************************************************************************80

double *circle_ppr2imp_2d ( double p1[], double p2[], double r )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_PPR2IMP_2D converts a circle from PPR to implicit form in 2D.
//
//  Discussion:
//
//    The PPR form of a circle in 2D is:
//
//      The circle of radius R passing through points P1 and P2.
//
//    The implicit form of a circle in 2D is:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = R * R
//
//    There may be zero, one, or two circles that satisfy the
//    requirements of the PPR form.
//
//    If there is no such circle, then the two "solutions" are set to
//    the midpoint of P1 and P2.
//
//    If there is one circle, then the two solutions will be set equal
//    to the midpoint of P1 and P2.
//
//    If there are two distinct circles, then (PC[0],PC[1]) is the first center,
//    and (PC[2],PC[3]) is the second.
//
//    This calculation is equivalent to finding the intersections of
//    circles of radius R at points P1 and P2.
//
//  Modified:
//
//    11 November 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[2], P2[2], two points on the circle.
//
//    Input, double R, the radius of the circle.
//
//    Output, double PC[2*2], the centers of the two circles.
//
{
# define DIM_NUM 2

  double dist;
  double h;
  int i;
  int j;
  double *pc;

  pc = new double[DIM_NUM*2];
//
//  Compute the distance from P1 to P2.
//
  dist = sqrt ( pow ( p2[0] - p1[0], 2 ) + pow ( p2[1] - p1[1], 2 ) );
//
//  If R is smaller than DIST, we don't have a circle.
//
  if ( 2.0 * r < dist )
  {
    for ( j = 0; j < 2; j++ )
    {
      for ( i = 0; i < DIM_NUM; i++ )
      {
        pc[i+j*DIM_NUM] = 0.5 * ( p1[i] + p2[i] );
      }
    }
  }
//
//  H is the distance from the midpoint of (P1,P2) to the center.
//
  h = sqrt ( ( r + 0.5 * dist ) * ( r - 0.5 * dist ) );
//
//  The center is found by going midway between P1 and P2, and then
//  H units in the unit perpendicular direction.
//
//  We actually have two choices for the normal direction.
//
  pc[0+0*DIM_NUM] = 0.5 * ( p2[0] + p1[0] ) + h * ( p2[1] - p1[1] ) / dist;
  pc[1+0*DIM_NUM] = 0.5 * ( p2[1] + p1[1] ) - h * ( p2[0] - p1[0] ) / dist;

  pc[0+1*DIM_NUM] = 0.5 * ( p2[0] + p1[0] ) - h * ( p2[1] - p1[1] ) / dist;
  pc[1+1*DIM_NUM] = 0.5 * ( p2[1] + p1[1] ) + h * ( p2[0] - p1[0] ) / dist;

  return pc;
# undef DIM_NUM
}
//****************************************************************************80

double circle_sector_area_2d ( double r, double pc[2], double theta1,
  double theta2 )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_SECTOR_AREA_2D computes the area of a circular sector in 2D.
//
//  Discussion:
//
//    A circular sector is formed by a circular arc, and the two straight line 
//    segments that join its ends to the center of the circle.
//
//    A circular sector is defined by
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//    and
//
//      Theta1 <= Theta <= Theta2
//
//  Modified:
//
//    10 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the coordinates of the center of the circle.
//
//    Input, double THETA1, THETA2, the angles of the first and last points
//    on the circular arc.
//
//    Output, double CIRCLE_SECTOR_AREA_2D, the area of the circle.
//
{
  double area;

  area = 0.5 * r * r * ( theta2 - theta1 );

  return area;
}
//****************************************************************************80

double *circle_sector_centroid_2d ( double r, double pc[2], double theta1, 
  double theta2 )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_SECTOR_CENTROID_2D returns the centroid of a circular sector in 2D.
//
//  Discussion:
//
//    A circular sector is formed by a circular arc, and the two straight line
//    segments that join its ends to the center of the circle.
//
//    A circular sector is defined by
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//    and
//
//      Theta1 <= Theta <= Theta2
//
//  Modified:
//
//    30 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Adrian Bowyer, John Woodwark,
//    A Programmer's Geometry,
//    Butterworths, 1983.
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the coordinates of the center of the circle.
//
//    Input, double THETA1, THETA2, the angles of the first and last points
//    on the circular arc.
//
//    Output, double CIRCLE_SECTOR_CENTROID_2D[2], the coordinates 
//    of the centroid of the sector.
//
{
# define DIM_NUM 2

  double *centroid;
  double d;
  double theta;

  theta = theta2 - theta1;

  if ( theta == 0.0 )
  {
    d = 2.0 * r / 3.0;
  }
  else
  {
    d = 4.0 * r * sin ( 0.5 * theta ) / ( 3.0 * theta );
  }

  centroid = new double[DIM_NUM];

  centroid[0] = pc[0] + d * cos ( theta );
  centroid[1] = pc[1] + d * sin ( theta );

  return centroid;
# undef DIM_NUM
}
//****************************************************************************80

bool circle_sector_contains_point_2d ( double r, double pc[2], double theta1, 
  double theta2, double p[2] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_SECTOR_CONTAINS_POINT_2D : is a point inside a circular sector?
//
//  Discussion:
//
//    A circular sector is formed by a circular arc, and the two straight line 
//    segments that join its ends to the center of the circle.
//
//    A circular sector is defined by
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//    and
//
//      Theta1 <= Theta <= Theta2
//
//  Modified:
//
//    22 October 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the center of the circle.
//
//    Input, double THETA1, THETA2, the angles defining the arc,
//    in radians.  Normally, THETA1 < THETA2.
//
//    Input, double P[2], the point to be checked.
//
//    Output, logical CIRCLE_SECTOR_CONTAINS_POINT_2D, is TRUE if the point is 
//    inside or on the circular sector.
//
{
# define DIM_NUM 2

  bool inside;
  double pi = 3.141592653589793;
  double theta;

  inside = false;
//
//  Is the point inside the (full) circle?
//
  if ( pow ( p[0] - pc[0], 2 ) + pow ( p[1] - pc[1], 2 ) <= r * r )
  {
//
//  Is the point's angle within the arc's range?
//  Try to force the angles to lie between 0 and 2 * PI.
//
    theta = atan4 ( p[1] - pc[1], p[0] - pc[0] );

    if ( r8_modp ( theta - theta1, 2.0 * pi ) <= 
         r8_modp ( theta2 - theta1, 2.0 * pi ) )
    {
      inside = true;
    }
  }

  return inside;
# undef DIM_NUM
}
//****************************************************************************80

void circle_sector_print_2d ( double r, double pc[2], double theta1, double theta2 )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_SECTOR_PRINT_2D prints a circular sector in 2D.
//
//  Discussion:
//
//    A circular sector is formed by a circular arc, and the two straight line 
//    segments that join its ends to the center of the circle.
//
//    A circular sector is defined by
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//    and
//
//      Theta1 <= Theta <= Theta2
//
//  Modified:
//
//    08 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the center of the circle.
//
//    Input, double THETA1, THETA2, the angles defining the arc,
//    in radians.  Normally, THETA1 < THETA2.
//
{
# define DIM_NUM 2

  cout << "\n";
  cout << "  Circular sector definition:\n";
  cout << "\n";
  cout << "    Radius = " << setw(12) << r << "\n";
  cout << "    Center = " << setw(12) << pc[0] 
       << "  "            << setw(12) << pc[1] << "\n";
  cout << "    Theta  = " << setw(12) << theta1
       << "  "            << setw(12) << theta2 << "\n";

  return;
# undef DIM_NUM
}
//****************************************************************************80

double circle_triangle_area_2d ( double r, double pc[2], double theta1, 
  double theta2 )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_TRIANGLE_AREA_2D returns the area of a circle triangle in 2D.
//
//  Discussion:
//
//    A circle triangle is formed by drawing a circular arc, and considering
//    the triangle formed by the endpoints of the arc plus the center of
//    the circle.
//
//    Note that for angles greater than PI, the triangle will actually
//    have NEGATIVE area.
//
//  Modified:
//
//    10 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R, the radius of the circle.
//
//    Input, double PC[2], the center of the circle.
//
//    Input, double THETA1, THETA2, the angles defining the arc,
//    in radians.  Normally, THETA1 < THETA2.
//
//    Output, double AREA, the (signed) area of the triangle.
//
{
  double area;

  area = 0.5 * r * r * sin ( theta2 - theta1 );

  return area;
}
//****************************************************************************80

void circle_triple_angles_2d ( double r1, double r2, double r3, double *angle1, 
  double *angle2, double *angle3 )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLE_TRIPLE_ANGLE_2D returns an angle formed by three circles in 2D.
//
//  Discussion:
//
//    A circle triple is a set of three tangent circles.  We assume
//    that no circle is contained in another.
//
//    We consider the triangle formed by joining the centers of the circles.
//
//  Modified:
//
//    02 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Kenneth Stephenson,
//    Circle Packing, The Theory of Discrete Analytic Functions,
//    Cambridge, 2005.
//
//  Parameters:
//
//    Input, double R1, R2, R3, the radii of the circles.
//
//    Input, double *ANGLE1, *ANGLE2, *ANGLE3, the angles
//    in the triangle.
//
{
  *angle1 = arc_cosine ( 
    pow ( r1 + r2, 2 ) + pow ( r1 + r3, 2 ) - pow ( r2 + r3, 2 ) ) / 
    ( 2.0 * ( r1 + r2 ) * ( r1 + r3 ) );

  *angle2 = arc_cosine ( 
    pow ( r2 + r3, 2 ) + pow ( r2 + r1, 2 ) - pow ( r3 + r1, 2 ) ) / 
    ( 2.0 * ( r2 + r3 ) * ( r2 + r1 ) );

  *angle3 = arc_cosine (
    pow ( r3 + r1, 2 ) + pow ( r3 + r2, 2 ) - pow ( r1 + r2, 2 ) ) /
    ( 2.0 * ( r3 + r1 ) * ( r3 + r2 ) );

  return;
}
//****************************************************************************80

void circles_imp_int_2d ( double r1, double pc1[2], double r2, double pc2[2],
  int *int_num, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    CIRCLES_IMP_INT_2D: finds the intersection of two implicit circles in 2D.
//
//  Discussion:
//
//    Two circles can intersect in 0, 1, 2 or infinitely many points.
//
//    The 0 and 2 intersection cases are numerically robust; the 1 and
//    infinite intersection cases are numerically fragile.  The routine
//    uses a tolerance to try to detect the 1 and infinite cases.
//
//    An implicit circle in 2D satisfies the equation:
//
//      pow ( P[0] - PC[0], 2 ) + pow ( P[1] - PC[1], 2 ) = pow ( R, 2 )
//
//    Thanks to Mario Pintaric for pointing out, on 13 March 2006,
//    a place where (R1-R2) had been mistakenly written as (R1-R1).
//
//  Modified:
//
//    13 March 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double R1, the radius of the first circle.
//
//    Input, double PC1[2], the coordinates of the center of the first circle.
//
//    Input, double R2, the radius of the second circle.
//
//    Input, double PC2[2], the coordinates of the center of the second circle.
//
//    Output, int *INT_NUM, the number of intersecting points found.
//    INT_NUM will be 0, 1, 2 or 3.  3 indicates that there are an infinite
//    number of intersection points.
//
//    Output, double P[2*2], if INT_NUM is 1 or 2, the
//    coordinates of the intersecting points.
//
{
  double distsq;
  double root;
  double sc1;
  double sc2;
  double t1;
  double t2;
  double tol;

  tol = r8_epsilon ( );

  p[0+0*2] = 0.0;
  p[1+0*2] = 0.0;
  p[0+1*2] = 0.0;
  p[1+1*2] = 0.0;
//
//  Take care of the case in which the circles have the same center.
//
  t1 = ( r8_abs ( pc1[0] - pc2[0] ) + r8_abs ( pc1[1] - pc2[1] ) ) / 2.0;
  t2 = ( r8_abs ( pc1[0] ) + r8_abs ( pc2[0] ) 
       + r8_abs ( pc1[1] ) + r8_abs ( pc2[1] ) + 1.0 ) / 5.0;

  if ( t1 <= tol * t2 )
  {
    t1 = r8_abs ( r1 - r2 );
    t2 = ( r8_abs ( r1 ) + r8_abs ( r2 ) + 1.0 ) / 3.0;

    if ( t1 <= tol * t2 )
    {
      *int_num = 3;
    }
    else
    {
      *int_num = 0;
    }
    return;
  }

  distsq = ( pc1[0] - pc2[0] ) * ( pc1[0] - pc2[0] ) 
         + ( pc1[1] - pc2[1] ) * ( pc1[1] - pc2[1] );

  root = 2.0 * ( r1 * r1 + r2 * r2 ) * distsq - distsq * distsq
    - ( r1 - r2 ) * ( r1 - r2 ) * ( r1 + r2 ) * ( r1 + r2 );

  if ( root < -tol )
  {
    *int_num = 0;
    return;
  }

  sc1 = ( distsq - ( r2 * r2 - r1 * r1 ) ) / distsq;

  if ( root < tol )
  {
    *int_num = 1;
    p[0+0*2] = pc1[0] + 0.5 * sc1 * ( pc2[0] - pc1[0] );
    p[1+0*2] = pc1[1] + 0.5 * sc1 * ( pc2[1] - pc1[1] );
    return;
  }

  sc2 = sqrt ( root ) / distsq;

  *int_num = 2;

  p[0+0*2] = pc1[0] + 0.5 * sc1 * ( pc2[0] - pc1[0] ) 
                    - 0.5 * sc2 * ( pc2[1] - pc1[1] );
  p[1+0*2] = pc1[1] + 0.5 * sc1 * ( pc2[1] - pc1[1] ) 
                    + 0.5 * sc2 * ( pc2[0] - pc1[0] );

  p[0+1*2] = pc1[0] + 0.5 * sc1 * ( pc2[0] - pc1[0] ) 
                    + 0.5 * sc2 * ( pc2[1] - pc1[1] );
  p[1+1*2] = pc1[1] + 0.5 * sc1 * ( pc2[1] - pc1[1] ) 
                    - 0.5 * sc2 * ( pc2[0] - pc1[0] );

  return;
}
//****************************************************************************80

double cone_area_3d ( double h, double r )

//****************************************************************************80
//
//  Purpose:
//
//    CONE_AREA_3D computes the surface area of a right circular cone in 3D.
//
//  Modified:
//
//    22 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double H, R, the height of the cone, and the radius of the
//    circle that forms the base of the cone.
//
//    Output, double CONE_AREA_3D, the surface area of the cone.
//
{
  double area;
  double pi = 3.141592653589793;

  area = pi * r * sqrt ( h * h + r * r );

  return area;
}
//****************************************************************************80

double *cone_centroid_3d ( double r, double pc[3], double pt[3] )

//****************************************************************************80
//
//  Purpose:
//
//    CONE_CENTROID_3D returns the centroid of a cone in 3D.
//
//  Modified:
//
//    02 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Adrian Bowyer, John Woodwark,
//    A Programmer's Geometry,
//    Butterworths, 1983.
//
//  Parameters:
//
//    Input, double R, the radius of the circle at the base of the cone.
//
//    Input, double PC[3], the coordinates of the center of the circle.
//
//    Input, double PT[3], the coordinates of the tip of the cone.
//
//    Output, double CONE_CENTROID_3D[3], the coordinates of the centroid of the cone.
//
{
  double *centroid;
  int dim_num = 3;
  int i;

  centroid = new double[dim_num];

  for ( i = 0; i < dim_num; i++ )
  {
    centroid[i] = 0.75 * pc[i] + 0.25 * pt[i];
  }

  return centroid;
}
//****************************************************************************80

double cone_volume_3d ( double h, double r )

//****************************************************************************80
//
//  Purpose:
//
//    CONE_VOLUME_3D computes the volume of a right circular cone in 3D.
//
//  Modified:
//
//    22 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double H, R, the height of the cone, and the radius of the
//    circle that forms the base of the cone.
//
//    Output, double CONE_VOLUME_3D, the volume of the cone.
//
{
  double pi = 3.141592653589793;
  double volume;

  volume = pi * r * r * h / 3.0;

  return volume;
}
//****************************************************************************80

void conv3d ( char axis, double theta, int n, double cor3[], double cor2[] )

//****************************************************************************80
//
//  Purpose:
//
//    CONV3D converts 3D data to a 2D projection.
//
//  Discussion:
//
//    A "presentation angle" THETA is used to project the 3D point
//    (X3D, Y3D, Z3D) to the 2D projection (XVAL,YVAL).
//
//    If COR = 'X':
//
//      X2D = Y3D - sin ( THETA ) * X3D
//      Y2D = Z3D - sin ( THETA ) * X3D
//
//  Modified:
//
//    02 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char AXIS, the coordinate to be projected.
//    AXIS should be 'X', 'Y', or 'Z'.
//
//    Input, double THETA, the presentation angle in degrees.
//
//    Input, int N, the number of values to be projected.
//
//    Input, double COR3[3*N], the point coordinates.
//
//    Output, double COR2[2*N], the 2D projections.
//
{
  int j;
  double stheta;

  stheta = sin ( degrees_to_radians ( theta ) );

  if ( axis == 'X' || axis == 'x' )
  {
    for ( j = 0; j < n; j++ )
    {
      cor2[0+j*2] = cor3[2+j*2] - stheta * cor3[0+j*2];
      cor2[1+j*2] = cor3[3+j*2] - stheta * cor3[0+j*2];
    }
  }
  else if ( axis == 'Y' || axis == 'y' )
  {
    for ( j = 0; j < n; j++ )
    {
      cor2[0+j*2] = cor3[0+j*2] - stheta * cor3[1+j*2];
      cor2[1+j*2] = cor3[2+j*2] - stheta * cor3[1+j*2];
    }
  }
  else if ( axis == 'Z' || axis == 'z' )
  {
    for ( j = 0; j < n; j++ )
    {
      cor2[0+j*2] = cor3[0+j*2] - stheta * cor3[2+j*2];
      cor2[1+j*2] = cor3[1+j*2] - stheta * cor3[2+j*2];
    }
  }
  else
  {
    cout << "\n";
    cout << "CONV3D - Fatal error!\n";
    cout << "  Illegal coordinate index = " << axis << "\n";
    exit ( 1 );
  }

  return;
}
//****************************************************************************80

double cos_deg ( double angle )

//****************************************************************************80
//
//  Purpose:
//
//    COS_DEG returns the cosine of an angle given in degrees.
//
//  Modified:
//
//    22 May 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ANGLE, the angle, in degrees.
//
//    Output, double COS_DEG, the cosine of the angle.
//
{
# define DEGREES_TO_RADIANS ( 3.141592653589793 / 180.0 )

  double angle_rad;
  double value;

  angle_rad = DEGREES_TO_RADIANS * angle;

  value = cos ( angle_rad );

  return value;
# undef DEGREES_TO_RADIANS
}
//****************************************************************************80

double cot_deg ( double angle )

//****************************************************************************80
//
//  Purpose:
//
//    COT_DEG returns the cotangent of an angle given in degrees.
//
//  Modified:
//
//    22 May 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ANGLE, the angle, in degrees.
//
//    Output, double COT_DEG, the cotangent of the angle.
//
{
# define DEGREES_TO_RADIANS ( 3.141592653589793 / 180.0 )

  double angle_rad;
  double value;

  angle_rad = DEGREES_TO_RADIANS * angle;

  value = cos ( angle_rad ) / sin ( angle_rad );

  return value;
# undef DEGREES_TO_RADIANS
}
//****************************************************************************80

double cot_rad ( double angle )

//****************************************************************************80
//
//  Purpose:
//
//    COT_RAD returns the cotangent of an angle.
//
//  Modified:
//
//    21 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ANGLE, the angle, in radians.
//
//    Output, double COT_RAD, the cotangent of the angle.
//
{
  double value;

  value = cos ( angle ) / sin ( angle );

  return value;
}
//****************************************************************************80

double csc_deg ( double angle )

//****************************************************************************80
//
//  Purpose:
//
//    CSC_DEG returns the cosecant of an angle given in degrees.
//
//  Modified:
//
//    22 May 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ANGLE, the angle, in degrees.
//
//    Output, double CSC_DEG, the cosecant of the angle.
//
{
# define DEGREES_TO_RADIANS ( 3.141592653589793 / 180.0 )

  double angle_rad;
  double value;

  angle_rad = DEGREES_TO_RADIANS * angle;

  value = 1.0 / sin ( angle_rad );

  return value;
# undef DEGREES_TO_RADIANS
}
//****************************************************************************80

void cube_shape_3d ( int point_num, int face_num, int face_order_max, 
  double point_coord[], int face_order[], int face_point[] )

//****************************************************************************80
//
//  Purpose:
//
//    CUBE_SHAPE_3D describes a cube in 3D.
//
//  Discussion:
//
//    The vertices lie on the unit sphere.
//
//    The dual of the octahedron is the cube.
//
//  Modified:
//
//    10 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, int FACE_NUM, the number of faces.
//
//    Input, int FACE_ORDER_MAX, the maximum number of vertices 
//    per face.
//
//    Output, double POINT_COORD[3*POINT_NUM], the point coordinates.
//
//    Output, int FACE_ORDER[FACE_NUM], the number of vertices per face.
//
//    Output, int FACE_POINT[FACE_ORDER_MAX*FACE_NUM]; FACE_POINT(I,J)
//    contains the index of the I-th point in the J-th face.  The
//    points are listed in the counter clockwise direction defined
//    by the outward normal at the face.
//
{
# define DIM_NUM 3

  double a = sqrt ( 1.0 / 3.0 );

  static int face_order_save[6] = {
    4, 4, 4, 4, 4, 4 };
  static int face_point_save[4*6] = {
     1, 4, 3, 2, 
     1, 2, 6, 5, 
     2, 3, 7, 6, 
     3, 4, 8, 7, 
     1, 5, 8, 4, 
     5, 6, 7, 8 }; 
  static double point_coord_save[DIM_NUM*8] = {
     -a, -a, -a, 
      a, -a, -a, 
      a,  a, -a, 
     -a,  a, -a, 
     -a, -a,  a, 
      a, -a,  a, 
      a,  a,  a, 
     -a,  a,  a };

  i4vec_copy ( face_num, face_order_save, face_order );
  i4vec_copy ( face_order_max*face_num, face_point_save, face_point );
  r8vec_copy ( DIM_NUM*point_num, point_coord_save, point_coord );

  return;
# undef DIM_NUM
}
//****************************************************************************80

void cube_size_3d ( int *point_num, int *edge_num, int *face_num, 
  int *face_order_max )

//****************************************************************************80
//
//  Purpose:
//
//    CUBE_SIZE_3D gives "sizes" for a cube in 3D.
//
//  Modified:
//
//    22 July 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, int *POINT_NUM, the number of points.
//
//    Output, int *EDGE_NUM, the number of edges.
//
//    Output, int *FACE_NUM, the number of faces.
//
//    Output, int *FACE_ORDER_MAX, the maximum order of any face.
//
{
  *point_num = 8;
  *edge_num = 12;
  *face_num = 6;
  *face_order_max = 4;

  return;
}
//****************************************************************************80

double cylinder_point_dist_3d ( double p1[3], double p2[3], double r, 
  double p[3] )

//****************************************************************************80
//
//  Purpose:
//
//    CYLINDER_POINT_DIST_3D determines the distance from a cylinder to a point in 3D.
//
//  Discussion:
//
//    We are computing the distance to the SURFACE of the cylinder.
//
//    The surface of a (right) (finite) cylinder in 3D is defined by an axis,
//    which is the line segment from point P1 to P2, and a radius R.  The points 
//    on the surface of the cylinder are:
//    * points at a distance R from the line through P1 and P2, and whose nearest
//      point on the line through P1 and P2 is strictly between P1 and P2, 
//    PLUS
//    * points at a distance less than or equal to R from the line through P1
//      and P2, whose nearest point on the line through P1 and P2 is either 
//      P1 or P2.
//
//  Modified:
//
//    23 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[3], P2[3], the first and last points
//    on the axis line of the cylinder.
//
//    Input, double R, the radius of the cylinder.
//
//    Input, double P[3], the point.
//
//    Output, double CYLINDER_POINT_DIST_3D, the distance from the point 
//    to the cylinder.
//
{
# define DIM_NUM 3

  double axis[DIM_NUM];
  double axis_length = 0.0;
  double distance = 0.0;
  int i;
  double off_axis_component = 0.0;
  double p_dot_axis = 0.0;
  double p_length = 0.0;
  double v1[DIM_NUM];

  for ( i = 0; i < DIM_NUM; i++ )
  {
    axis[i] = p2[i] - p1[i];
  }

  axis_length = r8vec_length ( DIM_NUM, axis );

  if ( axis_length == 0.0 )
  {
    distance = -r8_huge ( );
    return distance;
  }

  for ( i = 0; i < DIM_NUM; i++ )
  {
    axis[i] = axis[i] / axis_length;
  }

  p_dot_axis = 0.0;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    p_dot_axis = p_dot_axis + ( p[i] - p1[i] ) * axis[i];
  }
//
//  Case 1: Below bottom cap.
//
  if ( p_dot_axis <= 0.0 )
  {
    distance = disk_point_dist_3d ( p1, r, axis, p );
  }
//
//  Case 2: between cylinder planes.
//
  else if ( p_dot_axis <= axis_length )
  {
    for ( i = 0; i < DIM_NUM; i++ )
    {
      v1[i] = p[i] - p1[i];
    }
    p_length = r8vec_length ( DIM_NUM, v1 );
    off_axis_component = sqrt ( pow ( p_length, 2 ) - pow ( p_dot_axis, 2 ) );

    distance = r8_abs ( off_axis_component - r );

    if ( off_axis_component < r )
    {
      distance = r8_min ( distance, axis_length - p_dot_axis );
      distance = r8_min ( distance, p_dot_axis );
    }
  }
//
//  Case 3: Above the top cap.
//  
  else if ( axis_length < p_dot_axis )
  {
    distance = disk_point_dist_3d ( p2, r, axis, p );
  }

  return distance;
# undef DIM_NUM
}
//****************************************************************************80

double cylinder_point_dist_signed_3d ( double p1[3], double p2[3], double r, 
  double p[3] )

//****************************************************************************80
//
//  Purpose:
//
//    CYLINDER_POINT_DIST_SIGNED_3D: signed distance from cylinder to point in 3D.
//
//  Discussion:
//
//    We are computing the signed distance to the SURFACE of the cylinder.
//
//    The surface of a (right) (finite) cylinder in 3D is defined by an axis,
//    which is the line segment from point P1 to P2, and a radius R.  The points 
//    on the surface of the cylinder are:
//    * points at a distance R from the line through P1 and P2, and whose nearest
//      point on the line through P1 and P2 is strictly between P1 and P2, 
//    PLUS
//    * points at a distance less than or equal to R from the line through P1
//      and P2, whose nearest point on the line through P1 and P2 is either 
//      P1 or P2.
//
//    Points inside the surface have a negative distance.
//    Points on the surface have a zero distance.
//    Points outside the surface have a positive distance.
//
//  Modified:
//
//    26 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[3], P2[3], the first and last points
//    on the axis line of the cylinder.
//
//    Input, double R, the radius of the cylinder.
//
//    Input, double P[3], the point.
//
//    Output, double CYLINDER_POINT_DIST_SIGNED_3D, the signed distance from the point 
//    to the cylinder.
//
{
# define DIM_NUM 3

  double axis[DIM_NUM];
  double axis_length = 0.0;
  double distance = 0.0;
  int i;
  double off_axis_component = 0.0;
  double p_dot_axis = 0.0;
  double p_length = 0.0;
  double v1[DIM_NUM];

  for ( i = 0; i < DIM_NUM; i++ )
  {
    axis[i] = p2[i] - p1[i];
  }

  axis_length = r8vec_length ( DIM_NUM, axis );

  if ( axis_length == 0.0 )
  {
    distance = -r8_huge ( );
    return distance;
  }

  for ( i = 0; i < DIM_NUM; i++ )
  {
    axis[i] = axis[i] / axis_length;
  }

  p_dot_axis = 0.0;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    p_dot_axis = p_dot_axis + ( p[i] - p1[i] ) * axis[i];
  }
//
//  Case 1: Below bottom cap.
//
  if ( p_dot_axis <= 0.0 )
  {
    distance = disk_point_dist_3d ( p1, r, axis, p );
  }
//
//  Case 2: between cylinder planes.
//
  else if ( p_dot_axis <= axis_length )
  {
    for ( i = 0; i < DIM_NUM; i++ )
    {
      v1[i] = p[i] - p1[i];
    }
    p_length = r8vec_length ( DIM_NUM, v1 );
    off_axis_component = sqrt ( pow ( p_length, 2 ) - pow ( p_dot_axis, 2 ) );

    distance = off_axis_component - r;

    if ( distance < 0.0 )
    {
      distance = r8_max ( distance, p_dot_axis - axis_length);
      distance = r8_max ( distance, -p_dot_axis );
    }
  }
//
//  Case 3: Above the top cap.
//  
  else if ( axis_length < p_dot_axis )
  {
    distance = disk_point_dist_3d ( p2, r, axis, p );
  }

  return distance;
# undef DIM_NUM
}
//****************************************************************************80

bool cylinder_point_inside_3d ( double p1[3], double p2[3], double r, double p[3] )

//****************************************************************************80
//
//  Purpose:
//
//    CYLINDER_POINT_INSIDE_3D determines if a cylinder contains a point in 3D.
//
//  Discussion:
//
//    The surface and interior of a (right) (finite) cylinder in 3D is defined 
//    by an axis, which is the line segment from point P1 to P2, and a 
//    radius R.  The points contained in the volume include:
//    * points at a distance less than or equal to R from the line through P1
//      and P2, whose nearest point on the line through P1 and P2 is, in fact,
//      P1, P2, or any point between them.
//
//  Modified:
//
//    23 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[3], P2[3], the first and last points
//    on the axis line of the cylinder.
//
//    Input, double R, the radius of the cylinder.
//
//    Input, double P[3], the point.
//
//    Output, bool CYLINDER_POINT_INSIDE_3D, is TRUE if the point is inside the cylinder.
//
{
# define DIM_NUM 3

  double axis[DIM_NUM];
  double axis_length;
  int i;
  bool inside;
  double off_axis_component;
  double p_dot_axis;
  double p_length;
  double v1[DIM_NUM];

  for ( i = 0; i < DIM_NUM; i++ )
  {
    axis[i] = p2[i] - p1[i];
  }

  axis_length = r8vec_length ( DIM_NUM, axis );

  if ( axis_length == 0.0 )
  {
    inside = false;
    return inside;
  }

  for ( i = 0; i < DIM_NUM; i++ )
  {
    axis[i] = axis[i] / axis_length;
  }

  p_dot_axis = 0.0;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    p_dot_axis = p_dot_axis + ( p[i] - p1[i] ) * axis[i];
  }
//
//  If the point lies below or above the "caps" of the cylinder, we're done.
//
  if ( p_dot_axis < 0.0 || axis_length < p_dot_axis )
  {
    inside = false;
  }
//
//  Otherwise, determine the distance from P to the axis.
//
  else
  {
    for ( i = 0; i < DIM_NUM; i++ )
    {
      v1[i] = p[i] - p1[i];
    }
    p_length = r8vec_length ( DIM_NUM, v1 );

    off_axis_component = sqrt ( pow ( p_length, 2 ) - pow ( p_dot_axis, 2 ) );

    if ( off_axis_component <= r )
    {
      inside = true;
    }
    else
    {
      inside = false;
    }
  }

  return inside;
# undef DIM_NUM
}
//****************************************************************************80

double *cylinder_point_near_3d ( double p1[3], double p2[3], double r, double p[3] )

//****************************************************************************80
//
//  Purpose:
//
//    CYLINDER_POINT_NEAR_3D determines the nearest point on a cylinder to a point in 3D.
//
//  Discussion:
//
//    We are computing the nearest point on the SURFACE of the cylinder.
//
//    The surface of a (right) (finite) cylinder in 3D is defined by an axis,
//    which is the line segment from point P1 to P2, and a radius R.  The points 
//    on the surface of the cylinder are:
//    * points at a distance R from the line through P1 and P2, and whose nearest
//      point on the line through P1 and P2 is strictly between P1 and P2, 
//    PLUS
//    * points at a distance less than or equal to R from the line through P1
//      and P2, whose nearest point on the line through P1 and P2 is either 
//      P1 or P2.
//
//  Modified:
//
//    19 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[3], P2[3], the first and last points
//    on the axis line of the cylinder.
//
//    Input, double R, the radius of the cylinder.
//
//    Input, double P[3], the point.
//
//    Output, double CYLINDER_POINT_NEAR_3D[3], the nearest point on the cylinder.
//
{
# define DIM_NUM 3

  double axial_component;
  double axis[DIM_NUM];
  double axis_length;
  double distance;
  int i;
  double *normal;
  double off_axis[DIM_NUM];
  double off_axis_component;
  double *pn;

  pn = new double[DIM_NUM];

  for ( i = 0; i < DIM_NUM; i++ )
  {
    axis[i] = p2[i] - p1[i];
  }
  axis_length = r8vec_length ( DIM_NUM, axis );
  for ( i = 0; i < DIM_NUM; i++ )
  {
    axis[i] = axis[i] / axis_length;
  }

  axial_component = 0.0;
  for ( i = 0; i < DIM_NUM; i++ )
  {
    axial_component = axial_component + ( p[i] - p1[i] ) * axis[i];
  }
  for ( i = 0; i < DIM_NUM; i++ )
  {
    off_axis[i] = p[i] - p1[i] - axial_component * axis[i];
  }

  off_axis_component = r8vec_length ( DIM_NUM, off_axis );
//
//  Case 1: Below bottom cap.
//
  if ( axial_component <= 0.0 )
  {
    if ( off_axis_component <= r )
    {
      for ( i = 0; i < DIM_NUM; i++ )
      {
        pn[i] = p1[i] + off_axis[i];
      }
    }
    else
    {
      for ( i = 0; i < DIM_NUM; i++ )
      {
        pn[i] = p1[i] + ( r / off_axis_component ) * off_axis[i];
      }
    }
  }
//
//  Case 2: between cylinder planes.
//
  else if ( axial_component <= axis_length )
  {
    if ( off_axis_component == 0.0 )
    {
      normal = r8vec_any_normal ( DIM_NUM, axis );
      for ( i = 0; i < DIM_NUM; i++ )
      {
        pn[i] = p[i] + r * normal[i];
      }
      delete [] normal;
    }
    else
    {
      distance = r8_abs ( off_axis_component - r );

      for ( i = 0; i < DIM_NUM; i++ )
      {
        pn[i] = p1[i] + axial_component * axis[i] 
              + ( r / off_axis_component ) * off_axis[i];
      }
      if ( off_axis_component < r )
      {
        if ( axis_length - axial_component < distance )
        {
          distance = axis_length - axial_component;
          for ( i = 0; i < DIM_NUM; i++ )
          {
            pn[i] = p2[i] + off_axis[i];
          }
        }
        if ( axial_component < distance )
        {
          distance = axial_component;
          for ( i = 0; i < DIM_NUM; i++ )
          {
            pn[i] = p1[i] + off_axis[i];
          }
        }
      }
    }
  }
//
//  Case 3: Above the top cap.
//  
  else if ( axis_length < axial_component )
  {
    if ( off_axis_component <= r )
    {
      for ( i = 0; i < DIM_NUM; i++ )
      {
        pn[i] = p2[i] + off_axis[i];
      }
    }
    else
    {
      for ( i = 0; i < DIM_NUM; i++ )
      {
        pn[i] = p2[i] + ( r / off_axis_component ) * off_axis[i];
      }
    }
  }

  return pn;
# undef DIM_NUM
}
//****************************************************************************80

double *cylinder_sample_3d ( double p1[3], double p2[3], double r, int n, 
  int *seed )

//****************************************************************************80
//
//  Purpose:
//
//    CYLINDER_SAMPLE_3D samples a cylinder in 3D.
//
//  Discussion:
//
//    We are sampling the interior of a right finite cylinder in 3D.
//
//    The interior of a (right) (finite) cylinder in 3D is defined by an axis,
//    which is the line segment from point P1 to P2, and a radius R.  The points
//    on or inside the cylinder are:
//    * points whose distance from the line through P1 and P2 is less than
//      or equal to R, and whose nearest point on the line through P1 and P2
//      lies (nonstrictly) between P1 and P2.
//
//  Modified:
//
//    27 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[3], P2[3], the first and last points
//    on the axis line of the cylinder.
//
//    Input, double R, the radius of the cylinder.
//
//    Input, int N, the number of sample points to compute.
//
//    Input/output, int *SEED, the random number seed.
//
//    Input, double CYLINDER_SAMPLE_3D[3*N], the sample points.
//
{
# define DIM_NUM 3

  double axis[DIM_NUM];
  double axis_length;
  int i;
  int j;
  double *p;
  double pi = 3.141592653589793;
  double radius;
  double theta;
  double v2[DIM_NUM];
  double v3[DIM_NUM];
  double z;
//
//  Compute the axis vector.
//
  for ( i = 0; i < DIM_NUM; i++ )
  {
    axis[i] = p2[i] - p1[i];
  }
  axis_length = r8vec_length ( DIM_NUM, axis );
  for ( i = 0; i < DIM_NUM; i++ )
  {
    axis[i] = axis[i] / axis_length;
  }
//
//  Compute vectors V2 and V3 that form an orthogonal triple with AXIS.
//
  plane_normal_basis_3d ( p1, axis, v2, v3 );
//
//  Assemble the randomized information.
//
  p = new double[DIM_NUM*n];

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < DIM_NUM; i++ )
    {
      radius = r * sqrt ( r8_uniform_01 ( seed ) );
      theta = 2.0 * pi * r8_uniform_01 ( seed );
      z = axis_length * r8_uniform_01 ( seed );

      p[i+j*DIM_NUM] = p1[i] 
        + z                      * axis[i] 
        + radius * cos ( theta ) * v2[i]   
        + radius * sin ( theta ) * v3[i];
    }
  }

  return p;
# undef DIM_NUM
}
//****************************************************************************80

double cylinder_volume_3d ( double p1[3], double p2[3], double r )

//****************************************************************************80
//
//  Purpose:
//
//    CYLINDER_VOLUME_3D determines the volume of a cylinder in 3D.
//
//  Discussion:
//
//    A (right) (finite) cylinder in 3D is the set of points
//    contained on or inside a circle of radius R, whose center
//    lies along the line segment from point P1 to P2, and whose
//    plane is perpendicular to that line segment.
//
//  Modified:
//
//    11 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double P1[3], P2[3], the first and last points
//    on the axis line of the cylinder.
//
//    Input, double R, the radius of the cylinder.
//
//    Output, double CYLINDER_VOLUME_3D, the volume of the cylinder.
//
{
# define DIM_NUM 3

  double h;
  double pi = 3.141592653589793;
  double volume;

  h = r8vec_distance ( DIM_NUM, p1, p2 );

  volume = pi * r * r * h;

  return volume;
# undef DIM_NUM
}
//****************************************************************************80

double degrees_to_radians ( double angle )

//****************************************************************************80
//
//  Purpose:
//
//    DEGREES_TO_RADIANS converts an angle from degrees to radians.
//
//  Modified:
//
//    27 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ANGLE, an angle in degrees.
//
//    Output, double DEGREES_TO_RADIANS, the equivalent angle
//    in radians.
//
{
  double pi = 3.141592653589793;
  double value;

  value = ( angle / 180.0 ) * pi;

  return value;
}
//****************************************************************************80

double dge_det ( int n, double a[], int pivot[] )

//****************************************************************************80
//
//  Purpose:
//
//    DGE_DET computes the determinant of a matrix factored by SGE_FA.
//
//  Discussion:
//
//    The doubly dimensioned array A is treated as a one dimensional vector,
//    stored by COLUMNS:  
//
//      A(0,0), A(1,0), A(2,0), ..., A(N-1,0) // A(1,0), A(1,1), ... A(N-1,1)
//
//    Entry A(I,J) is stored as A[I+J*N]
//
//  Modified:
//
//    04 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Jack Dongarra, James Bunch, Cleve Moler, Pete Stewart,
//    LINPACK User's Guide,
//    SIAM, 1979
//
//  Parameters:
//
//    Input, int N, the order of the matrix.
//    N must be positive.
//
//    Input, double A[N*N], the LU factors computed by DGE_FA.
//
//    Input, int PIVOT[N], as computed by DGE_FA.
//
//    Output, double DGE_DET, the determinant of the matrix.
//
{
  double det;
  int i;

  det = 1.0;

  for ( i = 0; i < n; i++ )
  {
    det = det * a[i+i*n];
    if ( pivot[i] != i+1 )
    {
      det = -det;
    }
  }

  return det;
}
//****************************************************************************80

int dge_fa ( int n, double a[], int pivot[] )

//****************************************************************************80
//
//  Purpose:
//
//    DGE_FA factors a general matrix.
//
//  Discussion:
//
//    DGE_FA is a simplified version of the LINPACK routine SGEFA.
//
//    The doubly dimensioned array A is treated as a one dimensional vector,
//    stored by COLUMNS:  
//
//      A(0,0), A(1,0), A(2,0), ..., A(N-1,0) // A(1,0), A(1,1), ... A(N-1,1)
//
//    Entry A(I,J) is stored as A[I+J*N]
//
//  Modified:
//
//    05 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Jack Dongarra, James Bunch, Cleve Moler, Pete Stewart,
//    LINPACK User's Guide,
//    SIAM, 1979
//
//  Parameters:
//
//    Input, int N, the order of the matrix.
//    N must be positive.
//
//    Input/output, double A[N*N], the matrix to be factored.
//    On output, A contains an upper triangular matrix and the multipliers
//    which were used to obtain it.  The factorization can be written
//    A = L * U, where L is a product of permutation and unit lower
//    triangular matrices and U is upper triangular.
//
//    Output, int PIVOT[N], a vector of pivot indices.
//
//    Output, int DGE_FA, singularity flag.
//    0, no singularity detected.
//    nonzero, the factorization failed on the DGE_FA-th step.
//
{
  int i;
  int ii;
  int info;
  int j;
  int k;
  int l;
  double t;

  info = 0;

  for ( k = 1; k <= n-1; k++ )
  {
//
//  Find L, the index of the pivot row.
//
    l = k;
    for ( i = k+1; i <= n; i++ )
    {
      if ( r8_abs ( a[l-1+(k-1)*n] ) < r8_abs ( a[i-1+(k-1)*n] ) )
      {
        l = i;
      }
    }

    pivot[k-1] = l;
//
//  If the pivot index is zero, the algorithm has failed.
//
    if ( a[l-1+(k-1)*n] == 0.0 )
    {
      info = k;
      cout << "\n";
      cout << "DGE_FA - Fatal error!\n";
      cout << "  Zero pivot on step " << info << "\n";
      return info;
    }
//
//  Interchange rows L and K if necessary.
//
    if ( l != k )
    {
      t              = a[l-1+(k-1)*n];
      a[l-1+(k-1)*n] = a[k-1+(k-1)*n];
      a[k-1+(k-1)*n] = t;
    }
//
//  Normalize the values that lie below the pivot entry A(K,K).
//
    for ( j = k+1; j <= n; j++ )
    {
      a[j-1+(k-1)*n] = -a[j-1+(k-1)*n] / a[k-1+(k-1)*n];
    }
//
//  Row elimination with column indexing.
//
    for ( j = k+1; j <= n; j++ )
    {
      if ( l != k )
      {
        t              = a[l-1+(j-1)*n];
        a[l-1+(j-1)*n] = a[k-1+(j-1)*n];
        a[k-1+(j-1)*n] = t;
      }

      for ( ii = k; ii < n; ii++ )
      {
        a[ii+(j-1)*n] = a[ii+(j-1)*n] + a[ii+(k-1)*n] * a[k-1+(j-1)*n];
      }
    }
  }

  pivot[n-1] = n;

  if ( a[n-1+(n-1)*n] == 0.0 )
  {
    info = n;
    cout << "\n";
    cout << "DGE_FA - Fatal error!\n";
    cout << "  Zero pivot on step " << info << "\n";
  }

  return info;
}
//****************************************************************************80

void dge_sl ( int n, double a[], int pivot[], double b[], int job )

//****************************************************************************80
//
//  Purpose:
//
//    DGE_SL solves a system factored by SGE_FA.
//
//  Discussion:
//
//    DGE_SL is a simplified version of the LINPACK routine SGESL.
//
//    The doubly dimensioned array A is treated as a one dimensional vector,
//    stored by COLUMNS:  
//
//      A(0,0), A(1,0), A(2,0), ..., A(N-1,0) // A(1,0), A(1,1), ... A(N-1,1)
//
//    Entry A(I,J) is stored as A[I+J*N]
//
//  Modified:
//
//    06 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the matrix.
//    N must be positive.
//
//    Input, double A[N*N], the LU factors from DGE_FA.
//
//    Input, int PIVOT[N], the pivot vector from DGE_FA.
//
//    Input/output, double B[N].
//    On input, the right hand side vector.
//    On output, the solution vector.
//
//    Input, int JOB, specifies the operation.
//    0, solve A * x = b.
//    nonzero, solve A' * x = b.
//
{
  int i;
  int k;
  int l;
  double t;
//
//  Solve A * x = b.
//
  if ( job == 0 )
  {
//
//  Solve PL * Y = B.
//
    for ( k = 1; k <= n-1; k++ )
    {
      l = pivot[k-1];

      if ( l != k )
      {
        t      = b[l-1];
        b[l-1] = b[k-1];
        b[k-1] = t;
      }

      for ( i = k+1; i <= n; i++ )
      {
        b[i-1] = b[i-1] + a[i-1+(k-1)*n] * b[k-1];
      }
    }
//
//  Solve U * X = Y.
//
    for ( k = n; 1 <= k; k-- )
    {
      b[k-1] = b[k-1] / a[k-1+(k-1)*n];
      for ( i = 1; i <= k-1; i++ )
      {
        b[i-1] = b[i-1] - a[i-1+(k-1)*n] * b[k-1];
      }
    }
//
//  Solve A' * X = B.
//
  }
  else
  {
//
//  Solve U' * Y = B.
//
    for ( k = 1; k <= n; k++ )
    {
      t = 0.0;
      for ( i = 1; i <= k-1; i++ )
      {
        t = t + b[i-1] * a[i-1+(k-1)*n];
      }
      b[k-1] = ( b[k-1] - t ) / a[k-1+(k-1)*n];
    }
//
//  Solve ( PL )' * X = Y.
//
    for ( k = n-1; 1 <= k; k-- )
    {
      t = 0.0;
      for ( i = k+1; i <= n; i++ )
      {
        t = t + b[i-1] * a[i-1+(k-1)*n];
      }
      b[k-1] = b[k-1] + t;

      l = pivot[k-1];

      if ( l != k )
      {
        t      = b[l-1];
        b[l-1] = b[k-1];
        b[k-1] = t;
      }
    }
  }

  return;
}
//****************************************************************************80

double *direction_pert_3d ( double sigma, double vbase[3], int *seed )

//****************************************************************************80
//
//  Purpose:
//
//    DIRECTION_PERT_3D randomly perturbs a direction vector in 3D.
//
//  Modified:
//
//    04 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double SIGMA, determines the strength of the perturbation.
//    SIGMA <= 0 results in a completely random direction.
//    1 <= SIGMA results in VBASE.
//    0 < SIGMA < 1 results in a perturbation from VBASE, which is
//    large when SIGMA is near 0, and small when SIGMA is near 1.
//
//    Input, double VBASE[3], the base direction vector, which should have
//    unit norm.
//
//    Input/output, int *SEED, a seed for the random number generator.
//
//    Output, double DIRECTION_PERT_3D[3], the perturbed vector, which will have unit norm.
//
{
# define DIM_NUM 3

  double dphi;
  double p[DIM_NUM];
  double phi;
  double pi = 3.141592653589793;
  double psi;
  double theta;
  double vdot;
  double *vran;
  double x;
//
//  0 <= SIGMA, just use the base vector.
//
  vran = new double[DIM_NUM];

  if ( 1.0 <= sigma )
  {
    r8vec_copy ( DIM_NUM, vbase, vran );
  }
  else if ( sigma <= 0.0 )
  {
    vdot = r8_uniform_01 ( seed );
    vdot = 2.0 * vdot - 1.0;
    phi = acos ( vdot );
    theta = r8_uniform_01 ( seed );
    theta = 2.0 * pi * theta;

    vran[0] = cos ( theta ) * sin ( phi );
    vran[1] = sin ( theta ) * sin ( phi );
    vran[2] = cos ( phi );
  }
  else
  {
    phi = acos ( vbase