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

using namespace std;

# include "geometry.H"

int main ( void );
void test0005 ( void );
void test001 ( void );
void test002 ( void );
void test0023 ( void );
void test0025 ( void );
void test003 ( void );
void test0032 ( void );
void test0035 ( void );
void test004 ( void );
void test0045 ( void );
void test005 ( void );
void test006 ( void );
void test007 ( void );
void test0075 ( void );
void test008 ( void );
void test0085 ( void );
void test0087 ( void );
void test009 ( void );

void test010 ( void );
void test011 ( void );
void test012 ( void );
void test0125 ( void );
void test0126 ( void );
void test0127 ( void );
void test013 ( void );
void test014 ( void );
void test015 ( void );
void test0155 ( void );
void test0156 ( void );
void test016 ( void );
void test0165 ( void );
void test017 ( void );
void test018 ( void );
void test0183 ( void );
void test0185 ( void );
void test019 ( void );

void test020 ( void );
void test0201 ( void );
void test02015 ( void );
void test0202 ( void );
void test0203 ( void );
void test02035 ( void );
void test0204 ( void );
void test0205 ( void );
void test021 ( void );
void test022 ( void );
void test023 ( void );
void test0232 ( void );
void test0234 ( void );
void test0235 ( void );
void test0236 ( void );
void test0238 ( void );
void test024 ( void );
void test0243 ( void );
void test0245 ( void );
void test025 ( void );
void test0255 ( void );
void test026 ( void );
void test027 ( void );
void test028 ( void );
void test029 ( void );

void test030 ( void );
void test031 ( void );
void test0315 ( void );
void test032 ( void );
void test0321 ( void );
void test0322 ( void );
void test0323 ( void );
void test0325 ( void );
void test0327 ( void );
void test033 ( void );
void test0335 ( void );
void test0336 ( void );
void test0337 ( void );
void test034 ( void );
void test0345 ( void );
void test0346 ( void );
void test035 ( void );
void test038 ( void );
void test0385 ( void );
void test03855 ( void );
void test0386 ( void );
void test039 ( void );

void test040 ( void );
void test041 ( void );
void test0415 ( void );
void test0416 ( void );
void test0418 ( void );
void test042 ( void );
void test043 ( void );
void test044 ( void );
void test045 ( void );
void test046 ( void );
void test047 ( void );
void test0475 ( void );
void test048 ( void );
void test0485 ( void );
void test049 ( void );
void test0493 ( void );
void test0495 ( void );

void test050 ( void );
void test051 ( void );
void test052 ( void );
void test053 ( void );
void test054 ( void );
void test055 ( void );
void test056 ( void );
void test057 ( void );
void test058 ( void );
void test059 ( void );

void test060 ( void );
void test061 ( void );
void test0615 ( void );
void test062 ( void );
void test063 ( void );
void test064 ( void );
void test065 ( void );
void test066 ( void );
void test068 ( void );
void test0685 ( void );

void test0755 ( void );
void test0757 ( void );
void test076 ( void );
void test0765 ( void );
void test078 ( void );
void test0782 ( void );
void test0784 ( void );
void test0786 ( void );
void test079 ( void );

void test080 ( void );
void test0801 ( void );
void test0803 ( void );
void test0805 ( void );
void test0807 ( void );
void test081 ( void );
void test082 ( void );
void test0825 ( void );
void test083 ( void );
void test084 ( void );
void test0844 ( void );
void test0845 ( void );
void test0846 ( void );
void test085 ( void );

void test170 ( void );
void test171 ( void );
void test1715 ( void );
void test171 ( void );
void test172 ( void );
void test173 ( void );
void test174 ( void );
void test1745 ( void );
void test1746 ( void );
void test175 ( void );
void test176 ( void );
void test177 ( void );
void test178 ( void );
void test1787 ( void );
void test036 ( void );
void test0365 ( void );
void test0366 ( void );
void test0367 ( void );
void test0368 ( void );
void test037 ( void );
void test179 ( void );

void test180 ( void );
void test1805 ( void );
void test181 ( void );
void test182 ( void );
void test183 ( void );
void test1835 ( void );
void test1837 ( void );
void test1838 ( void );
void test1839 ( void );
void test184 ( void );
void test185 ( void );
void test186 ( void );
void test187 ( void );
void test188 ( void );
void test189 ( void );
void test1893 ( void );
void test1895 ( void );

void test190 ( void );
void test191 ( void );
void test192 ( void );
void test193 ( void );
void test194 ( void );
void test195 ( void );
void test1955 ( void );
void test196 ( void );
void test197 ( void );
void test198 ( void );
void test199 ( void );

void test200 ( void );
void test201 ( void );
void test202 ( void );
void test203 ( void );
void test2031 ( void );
void test2032 ( void );
void test20321 ( void );
void test20322 ( void );
void test20323 ( void );
void test203232 ( void );
void test203233 ( void );
void test203234 ( void );
void test203235 ( void );
void test20324 ( void );
void test20325 ( void );
void test2033 ( void );
void test2035 ( void );
void test2036 ( void );
void test204 ( void );
void test205 ( void );
void test206 ( void );
void test20605 ( void );
void test2061 ( void );
void test2062 ( void );
void test209 ( void );
void test20655 ( void );
void test2066 ( void );
void test2094 ( void );
void test2101 ( void );
void test2067 ( void );
void test21015 ( void );
void test2068 ( void );
void test2069 ( void );
void test207 ( void );
void test2075 ( void );
void test208 ( void );

void test2102 ( void );
void test2070 ( void );
void test20701 ( void );
void test211 ( void );
void test2103 ( void );
void test2071 ( void );
void test20715 ( void );
void test2095 ( void );
void test2072 ( void );
void test2115 ( void );
void test212 ( void );
void test213 ( void );
void test219 ( void );

void test220 ( void );
void test221 ( void );
void test222 ( void );
void test2225 ( void );
void test223 ( void );
void test224 ( void );
void test2245 ( void );
void test225 ( void );
void test226 ( void );
void test227 ( void );

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

int main ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    MAIN is the main program for GEOMETRY_PRB.
//
//  Discussion:
//
//    GEOMETRY_PRB tests routines from the GEOMETRY library.
//
//  Modified:
//
//   29 October 2007
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );

  cout << "\n";
  cout << "GEOMETRY_PRB\n";
  cout << "  C++ version\n";
  cout << "\n";
  cout << "  Test the routines in the GEOMETRY library.\n";

  test0005 ( );
  test001 ( );
  test002 ( );
  test0023 ( );
  test0025 ( );
  test003 ( );
  test0032 ( );
  test0035 ( );
  test004 ( );
  test0045 ( );
  test005 ( );
  test006 ( );
  test007 ( );
  test0075 ( );
  test008 ( );
  test0085 ( );
  test0087 ( );
  test009 ( );

  test010 ( );
  test011 ( );
  test012 ( );
  test0125 ( );
  test0126 ( );
  test0127 ( );
  test013 ( );
  test014 ( );
  test015 ( );
  test0155 ( );
  test0156 ( );
  test016 ( );
  test0165 ( );
  test017 ( );
  test018 ( );
  test0183 ( );
  test0185 ( );
  test019 ( );

  test020 ( );
  test0201 ( );
  test02015 ( );
  test0202 ( );
  test0203 ( );
  test02035 ( );
  test0204 ( );
  test0205 ( );
  test021 ( );
  test022 ( );
  test023 ( );
  test0232 ( );
  test0234 ( );
  test0235 ( );
  test0236 ( );
  test0238 ( );
  test024 ( );
  test0243 ( );
  test0245 ( );
  test025 ( );
  test0255 ( );
  test026 ( );
  test027 ( );
  test028 ( );
  test029 ( );

  test030 ( );
  test031 ( );
  test0315 ( );
  test032 ( );
  test0321 ( );
  test0322 ( );
  test0323 ( );
  test0325 ( );
  test0327 ( );
  test033 ( );
  test0335 ( );
  test0336 ( );
  test0337 ( );
  test034 ( );
  test0345 ( );
  test0346 ( );
  test035 ( );
  test038 ( );
  test0385 ( );
  test03855 ( );
  test0386 ( );
  test039 ( );

  test040 ( );
  test041 ( );
  test0415 ( );
  test0416 ( );
  test0418 ( );
  test042 ( );
  test043 ( );
  test044 ( );
  test045 ( );
  test046 ( );
  test047 ( );
  test0475 ( );
  test048 ( );
  test0485 ( );
  test049 ( );
  test0493 ( );
  test0495 ( );

  test050 ( );
  test051 ( );
  test052 ( );
  test053 ( );
  test054 ( );
  test055 ( );
  test056 ( );
  test057 ( );
  test058 ( );
  test059 ( );

  test060 ( );
  test061 ( );
  test0615 ( );
  test062 ( );
  test063 ( );
  test064 ( );
  test065 ( );
  test066 ( );
  test068 ( );
  test0685 ( );

  test0755 ( );
  test0757 ( );
  test076 ( );
  test0765 ( );
  test078 ( );
  test0782 ( );
  test0784 ( );
  test0786 ( );
  test079 ( );

  test080 ( );
  test0803 ( );
  test0805 ( );
  test0807 ( );
  test081 ( );
  test082 ( );
  test0825 ( );
  test083 ( );
  test084 ( );
  test0844 ( );
  test0845 ( );
  test0846 ( );
  test085 ( );

  test170 ( );
  test171 ( );
  test1715 ( );
  test172 ( );
  test173 ( );
  test174 ( );
  test1745 ( );
  test1746 ( );
  test175 ( );
  test176 ( );
  test177 ( );
  test178 ( );
  test1787 ( );
  test036 ( );
  test0365 ( );
  test0366 ( );
  test0367 ( );
  test0368 ( );
  test037 ( );
  test179 ( );

  test180 ( );
  test1805 ( );
  test181 ( );
  test182 ( );
  test183 ( );
  test1835 ( );
  test1837 ( );
  test1838 ( );
  test1839 ( );
  test184 ( );
  test185 ( );
  test186 ( );
  test187 ( );
  test188 ( );
  test189 ( );
  test1893 ( );
  test1895 ( );

  test190 ( );
  test191 ( );
  test192 ( );
  test193 ( );
  test194 ( );
  test195 ( );
  test1955 ( );
  test196 ( );
  test197 ( );
  test198 ( );
  test199 ( );

  test200 ( );
  test201 ( );
  test202 ( );
  test203 ( );
  test2031 ( );
  test2032 ( );
  test20321 ( );
  test20322 ( );
  test20323 ( );
  test203232 ( );
  test203233 ( );
  test203234 ( );
  test203235 ( );
  test20324 ( );
  test20325 ( );
  test2033 ( );
  test2035 ( );
  test2036 ( );
  test204 ( );
  test205 ( );
  test206 ( );
  test20605 ( );
  test2061 ( );
  test2062 ( );
  test209 ( );
  test20655 ( );
  test2066 ( );
  test2094 ( );
  test2101 ( );
  test2067 ( );
  test21015 ( );
  test2068 ( );
  test2069 ( );
  test207 ( );
  test2075 ( );
  test208 ( );

  test2102 ( );
  test2070 ( );
  test20701 ( );
  test211 ( );
  test2103 ( );
  test2071 ( );
  test20715 ( );
  test2095 ( );
  test2072 ( );
  test2115 ( );
  test212 ( );
  test213 ( );
  test219 ( );

  test220 ( );
  test221 ( );
  test222 ( );
  test2225 ( );
  test223 ( );
  test224 ( );
  test2245 ( );
  test225 ( );
  test226 ( );
  test227 ( );

  cout << "\n";
  cout << "GEOMETRY_PRB\n";
  cout << "  Normal end of execution.\n";

  cout << "\n";
  timestamp ( );

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

void test0005 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0005 tests ANGLE_BOX_2D.
//
//  Modified:
//
//    19 July 2006
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double dist;
  double p1[DIM_NUM];
  double p2[DIM_NUM];
  double p3[DIM_NUM];
  double p4[DIM_NUM];
  double p5[DIM_NUM];

  cout << "\n";
  cout << "TEST0005\n";
  cout << "  ANGLE_BOX_2D\n";
  cout << "\n";
  cout << "  Compute P4 and P5, normal to\n";
  cout << "  line through P1 and P2, and\n";
  cout << "  line through P2 and P3,\n";
  cout << "  and DIST units from P2.\n";
//
//  These points define the lines
//    y = 0
//  and
//    y = 2x-6
//
  p1[0] = 0.0;
  p1[1] = 0.0;
  p2[0] = 3.0;
  p2[1] = 0.0;
  p3[0] = 4.0;
  p3[1] = 2.0;
  dist = 1.0;

  cout << "\n";
  cout << "  DIST " << setw(14) << dist  << "\n";
  cout << "  P1:  " << setw(14) << p1[0] << "  " << setw(14) << p1[1] << "\n";
  cout << "  P2:  " << setw(14) << p2[0] << "  " << setw(14) << p2[1] << "\n";
  cout << "  P3:  " << setw(14) << p3[0] << "  " << setw(14) << p3[1] << "\n";

  angle_box_2d ( dist, p1, p2, p3, p4, p5 );

  cout << "  P4:  " << setw(14) << p4[0] << "  " << setw(14) << p4[1] << "\n";
  cout << "  P5:  " << setw(14) << p5[0] << "  " << setw(14) << p5[1] << "\n";
//
//  These points define the lines
//    y = 0
//  and
//    y = 2x-6
//
  p1[0] = 0.0;
  p1[1] = 0.0;
  p2[0] = 3.0;
  p2[1] = 0.0;
  p3[0] = 2.0;
  p3[1] = -2.0;
  dist = 1.0;

  cout << "\n";
  cout << "  DIST " << setw(14) << dist  << "\n";
  cout << "  P1:  " << setw(14) << p1[0] << "  " << setw(14) << p1[1] << "\n";
  cout << "  P2:  " << setw(14) << p2[0] << "  " << setw(14) << p2[1] << "\n";
  cout << "  P3:  " << setw(14) << p3[0] << "  " << setw(14) << p3[1] << "\n";

  angle_box_2d ( dist, p1, p2, p3, p4, p5 );

  cout << "  P4:  " << setw(14) << p4[0] << "  " << setw(14) << p4[1] << "\n";
  cout << "  P5:  " << setw(14) << p5[0] << "  " << setw(14) << p5[1] << "\n";
//
//  By setting P1 = P2, we are asking to extend the line
//    y = 2x-6
//  from P3 to P2 through to the other side.
//
  p1[0] = 3.0;
  p1[1] = 0.0;
  p2[0] = 3.0;
  p2[1] = 0.0;
  p3[0] = 2.0;
  p3[1] = -2.0;
  dist = 1.0;

  cout << "\n";
  cout << "  DIST " << setw(14) << dist  << "\n";
  cout << "  P1:  " << setw(14) << p1[0] << "  " << setw(14) << p1[1] << "\n";
  cout << "  P2:  " << setw(14) << p2[0] << "  " << setw(14) << p2[1] << "\n";
  cout << "  P3:  " << setw(14) << p3[0] << "  " << setw(14) << p3[1] << "\n";

  angle_box_2d ( dist, p1, p2, p3, p4, p5 );

  cout << "  P4:  " << setw(14) << p4[0] << "  " << setw(14) << p4[1] << "\n";
  cout << "  P5:  " << setw(14) << p5[0] << "  " << setw(14) << p5[1] << "\n";

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

void test001 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST001 tests ANGLE_CONTAINS_RAY_2D.
//
//  Modified:
//
//    19 July 2006
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define TEST_NUM 6

  int angle;
  int angle_num = 12;
  double angle_rad;
  bool inside;
  double p[DIM_NUM];
  double p1[DIM_NUM];
  double p2[DIM_NUM];
  double p3[DIM_NUM];
  double pi = 3.141592653589793;
  int test;

  cout << "\n";
  cout << "TEST001\n";
  cout << "  ANGLE_CONTAINS_RAY_2D sees if a ray lies within an angle.\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    if ( test == 0 )
    {
      p1[0] = 1.0;
      p1[1] = 0.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] = 1.0;
      p3[1] = 1.0;
    }
    else if ( test == 1 )
    {
      p1[0] = 1.0;
      p1[1] = 0.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] = 0.0;
      p3[1] = 1.0;
    }
    else if ( test == 2 )
    {
      p1[0] =  1.0;
      p1[1] = -1.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] = 0.0;
      p3[1] = 1.0;
    }
    else if ( test == 3 )
    {
      p1[0] = 1.0;
      p1[1] = 0.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] = -1.0;
      p3[1] =  0.0;
    }
    else if ( test == 4 )
    {
      p1[0] = 1.0;
      p1[1] = 0.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] =  0.0;
      p3[1] = -1.0;
    }
    else if ( test == 5 )
    {
      p1[0] = 1.0;
      p1[1] = 0.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] =  1.0;
      p3[1] = -0.01;
    }

    r8vec_print ( DIM_NUM, p1, "  Vertex A" );
    r8vec_print ( DIM_NUM, p2, "  Vertex B" );
    r8vec_print ( DIM_NUM, p3, "  Vertex C" );

    cout << "\n";
    cout << "       X            Y       Inside?\n";
    cout << "\n";

    for ( angle = 0; angle <= angle_num; angle++ )
    {
      angle_rad = ( double ) ( angle ) * 2.0 * pi / ( double ) angle_num;

      p[0] = cos ( angle_rad );
      p[1] = sin ( angle_rad );

      inside = angle_contains_ray_2d ( p1, p2, p3, p );

      cout << "  " << setw(12) << p[0]
           << "  " << setw(12) << p[1]
           << "  " << inside << "\n";
    }

  }

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

void test002 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST002 tests ANGLE_DEG_2D and ANGLE_RAD_ND.
//
//  Modified:
//
//    19 July 2006
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  int i;
  int angle_num = 12;
  double temp1;
  double temp2;
  double temp3;
  double thetad;
  double thetar;
  double v1[DIM_NUM];
  double v2[DIM_NUM];
  double v3[DIM_NUM];

  cout << "\n";
  cout << "TEST002\n";
  cout << "  ANGLE_DEG_2D computes an angle;\n";
  cout << "  ANGLE_RAD_ND computes an angle.\n";
  cout << "\n";
  cout << "      X           Y          Theta         atan2  ANGLE_RAD_ND, ANGLE_DEG_2D\n";
  cout << "\n";

  v1[0] = 1.0;
  v1[1] = 0.0;
  v3[0] = 0.0;
  v3[1] = 0.0;

  for ( i = 0; i <= angle_num; i++ )
  {
    thetad = ( double ) ( i ) * 360.0 / ( double ) ( angle_num );
    thetar = degrees_to_radians ( thetad );

    v2[0] = cos ( thetar );
    v2[1] = sin ( thetar );

    temp1 = radians_to_degrees ( atan2 ( v2[1], v2[0] ) );

    temp2 = angle_rad_nd ( DIM_NUM, v1, v2 );

    temp3 = angle_deg_2d ( v1, v3, v2 );

    cout << "  " << setw(10) << v2[0]
         << "  " << setw(10) << v2[1]
         << "  " << setw(10) << thetad
         << "  " << setw(10) << temp1
         << "  " << setw(10) << temp2
         << "  " << setw(10) << temp3 << "\n";
  }

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

void test0023 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0023 tests ANGLE_HALF_2D;
//
//  Modified:
//
//    19 July 2006
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double angle_deg;
  double p1[DIM_NUM];
  double p2[DIM_NUM];
  double p3[DIM_NUM];
  double *p4;
  double r;

  cout << "\n";
  cout << "TEST0023\n";
  cout << "  ANGLE_HALF_2D computes the half angle between two rays;\n";
  cout << "  The angle is defined by the points (P1,P2,P3)\n";
  cout << "  or by the rays P2-->P3, P2-->P1\n";

  p2[0] = 5.0;
  p2[1] = 3.0;

  angle_deg = 75.0;
  r = 3.0;
  p1[0] = p2[0] + r * cos_deg ( angle_deg );
  p1[1] = p2[1] + r * sin_deg ( angle_deg );

  angle_deg = 15.0;
  r = 2.0;
  p3[0] = p2[0] + r * cos_deg ( angle_deg );
  p3[1] = p2[1] + r * sin_deg ( angle_deg );

  r8vec_print ( DIM_NUM, p1, "  Point P1:" );
  r8vec_print ( DIM_NUM, p2, "  Point P2:" );
  r8vec_print ( DIM_NUM, p3, "  Point P3:" );

  p4 = angle_half_2d ( p1, p2, p3 );

  r8vec_print ( DIM_NUM, p4,
    "  End point of unit ray from P2, defining half angle, P4:" );

  angle_deg = 45.0;
  r = 1.0;
  p4[0] = p2[0] + r * cos_deg ( angle_deg );
  p4[1] = p2[1] + r * sin_deg ( angle_deg );

  r8vec_print ( DIM_NUM, p4, "  Expected value of P4:" );

  delete [] p4;

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

void test0025 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0025 tests ANGLE_RAD_2D;
//
//  Modified:
//
//    19 July 2006
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define TEST_NUM 6

  double angle;
  double p1[DIM_NUM];
  double p2[DIM_NUM];
  double p3[DIM_NUM];
  int test;

  cout << "\n";
  cout << "TEST0025\n";
  cout << "  ANGLE_RAD_2D computes the angle between two rays;\n";
  cout << "\n";

  for ( test = 1; test <= TEST_NUM; test++ )
  {
    if ( test == 1 )
    {
      p1[0] = 1.0;
      p1[1] = 0.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] = 1.0;
      p3[1] = 1.0;
    }
    else if ( test == 2 )
    {
      p1[0] = 1.0;
      p1[1] = 0.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] = 0.0;
      p3[1] = 1.0;
    }
    else if ( test == 3 )
    {
      p1[0] = 1.0;
      p1[1] = -1.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] = 0.0;
      p3[1] = 1.0;
    }
    else if ( test == 4 )
    {
      p1[0] = 1.0;
      p1[1] = 0.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] = -1.0;
      p3[1] =  0.0;
    }
    else if ( test == 5 )
    {
      p1[0] = 1.0;
      p1[1] = 0.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] =  0.0;
      p3[1] = -1.0;
    }
    else if ( test == 6 )
    {
      p1[0] = 1.0;
      p1[1] = 0.0;

      p2[0] = 0.0;
      p2[1] = 0.0;

      p3[0] =  1.0;
      p3[1] = -0.01;
    }

    r8vec_print ( DIM_NUM, p1, "  Vertex A" );
    r8vec_print ( DIM_NUM, p2, "  Vertex B" );
    r8vec_print ( DIM_NUM, p3, "  Vertex C" );

    angle = angle_rad_2d ( p1, p2, p3 );

    cout << "\n";
    cout << "  Angle = " << angle << "\n";
    cout << "\n";
  }
  return;
# undef DIM_NUM
# undef TEST_NUM
}
//****************************************************************************80********80

void test003 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST003 tests ANGLE_RAD_3D;
//
//  Modified:
//
//    19 July 2006
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3
# define TEST_NUM 3

  int i;
  double p1[DIM_NUM];
  double p1_test[DIM_NUM*TEST_NUM] = {
    1.0, 0.0, 0.0,
    1.0, 2.0, 3.0,
    0.0, 0.0, 1.0 };
  double p2[DIM_NUM] = { 0.0, 0.0, 0.0 };
  double p3[DIM_NUM] = { 0.0, 0.0, 1.0 };
  double temp1;
  double temp2;
  int test;

  cout << "\n";
  cout << "TEST003\n";
  cout << "  ANGLE_RAD_3D computes an angle;\n";
  cout << "\n";
  cout << "         X           Y           Z   ANGLE_RAD_3D  (Degrees)\n";
  cout << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    for ( i = 0; i < DIM_NUM; i++ )
    {
      p1[i] = p1_test[i+test*DIM_NUM];
    }
    temp1 = angle_rad_3d ( p1, p2, p3 );
    temp2 = radians_to_degrees ( temp1 );

    for ( i = 0; i < DIM_NUM; i++ )
    {
      cout << "  " << setw(10) << p1[i];
    }
    cout << "  " << setw(10) << temp1
         << "  " << setw(10) << temp2 << "\n";
  }

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

void test0032 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0032 tests ANGLE_TURN_2D.
//
//  Modified:
//
//    14 March 2006
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define TEST_NUM 13

  double p1[DIM_NUM];
  double p2[DIM_NUM] = { 0.0, 0.0 };
  double p3[DIM_NUM] = { 1.0, 0.0 };
  double pi = 3.141592653589793;
  int test;
  double theta;
  double theta_degrees;
  double turn;

  cout << "\n";
  cout << "TEST0032\n";
  cout << "  ANGLE_TURN_2D computes the turning angle\n";
  cout << "  defined by the line segments [P1,P2] and [P2,P3].\n";

  cout << "\n";
  cout << "  Our three points are:\n";
  cout << "\n";
  cout << "    P1 = (C,S)\n";
  cout << "    P2 = (0,0)\n";
  cout << "    P3 = (1,0)\n";
  cout << "\n";
  cout << "  C = cosine ( theta ), S = sine ( theta ).\n";
  cout << "\n";
  cout << "  Test  Theta    Turn\n";
  cout << "\n";

  for ( test = 1; test <= TEST_NUM; test++ )
  {
    theta = 2.0 * pi * ( double ) ( test - 1 )
      / ( double ) ( TEST_NUM - 1 );

    theta_degrees = 360.0 * ( double ) ( test - 1 )
      / ( double ) ( TEST_NUM - 1 );

    p1[0] = cos ( theta );
    p1[1] = sin ( theta );

    turn = angle_turn_2d ( p1, p2, p3 );

    cout << "  " << setw(4)  << test
         << "  " << setw(5)  << theta_degrees
         << "  " << setw(14) << turn << "\n";
  }

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

void test0035 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0035 tests ANNULUS_SECTOR_CENTROID_2D.
//
//  Modified:
//
//    02 December 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double *centroid;
  double pc[DIM_NUM] = { 5.0, 3.0 };
  double r1 = 2.0;
  double r2 = 3.0;
  double theta1;
  double theta2;

  theta1 = degrees_to_radians ( 30.0 );
  theta2 = degrees_to_radians ( 60.0 );

  cout << "\n";
  cout << "TEST0035\n";
  cout << "  ANNULUS_SECTOR_CENTROID_2D computes the centroid of a\n";
  cout << "    circular annulus.\n";
  cout << "\n";
  cout << "  The circle has center        " << pc[0]
       << "  " << pc[1] << "\n";
  cout << "  The inner radius is R1 =     " << r1     << "\n";
  cout << "  The outer radius is R2 =     " << r2     << "\n";
  cout << "  The first angle is THETA1 =  " << theta1 << "\n";
  cout << "  The second angle is THETA2 = " << theta2 << "\n";

  centroid = annulus_sector_centroid_2d ( pc, r1, r2, theta1, theta2 );

  cout << "\n";
  cout << "  Centroid: "
       << setw(12) << centroid[0]
       << setw(12) << centroid[1] << "\n";

  delete [] centroid;

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

void test004 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST004 tests ARC_COSINE;
//
//  Modified:
//
//    19 July 2006
//
//  Author:
//
//    John Burkardt
//
{
# define TEST_NUM 9

  double temp1;
  double temp2;
  int test;
  double x;
  double x_test[TEST_NUM] = {
    5.0, 1.2, 1.0, 0.9, 0.5, 0.0, -0.9, -1.0, -1.01 };

  cout << "\n";
  cout << "TEST004\n";
  cout << "  ARC_COSINE computes an angle with a given cosine;\n";
  cout << "\n";
  cout << "  X  ARC_COSINE(X) (Degrees)\n";
  cout << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    x = x_test[test];

    temp1 = arc_cosine ( x );
    temp2 = radians_to_degrees ( temp1 );

    cout << "  " << setw(12) << x
         << "  " << setw(12) << temp1
         << "  " << setw(12) << temp2 << "\n";
  }

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

void test0045 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0045 tests ARC_SINE;
//
//  Modified:
//
//    07 August 2005
//
//  Author:
//
//    John Burkardt
//
{
# define TEST_NUM 9

  double temp1;
  double temp2;
  int test;
  double x;
  double x_test[TEST_NUM] = {
    5.0, 1.2, 1.0, 0.9, 0.5, 0.0, -0.9, -1.0, -1.01 };

  cout << "\n";
  cout << "TEST0045\n";
  cout << "  ARC_SINE computes an angle with a given sine;\n";
  cout << "\n";
  cout << "  X  ARC_SINE(X)  (Degrees)\n";
  cout << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    x = x_test[test];

    temp1 = arc_sine ( x );
    temp2 = radians_to_degrees ( temp1 );

    cout << "  " << setw(12) << x
         << "  " << setw(12) << temp1
         << "  " << setw(12) << temp2 << "\n";
  }

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

void test005 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST005 tests ATAN4;
//
//  Modified:
//
//    19 July 2006
//
//  Author:
//
//    John Burkardt
//
{
# define TEST_NUM 8

  double temp1;
  double temp2;
  double temp3;
  int test;
  double x;
  double x_test[TEST_NUM] = { 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 0.0 };
  double y;
  double y_test[TEST_NUM] = { 0.0, 1.0, 2.0, 0.0, -1.0, -1.0, -1.0, -1.0 };

  cout << "\n";
  cout << "TEST005\n";
  cout << "  ATAN4 computes an angle with a given tangent.\n";
  cout << "\n";
  cout << "  X, Y, ATAN(Y/X), ATAN2(Y,X), ATAN4(Y,X)\n";
  cout << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    x = x_test[test];
    y = y_test[test];

    if ( x != 0.0 )
    {
      temp1 = atan ( y / x );
    }
    else
    {
      temp1 = r8_huge ( );
    }

    temp2 = atan2 ( y, x );
    temp3 = atan4 ( y, x );

    cout << "  " << setw(12) << x
         << "  " << setw(12) << y
         << "  " << setw(12) << temp1
         << "  " << setw(12) << temp2
         << "  " << setw(12) << temp3<< "\n";
  }

  cout << "\n";
  cout << "  Repeat, but display answers in degrees.\n";
  cout << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    x = x_test[test];
    y = y_test[test];

    if ( x != 0.0 )
    {
      temp1 = radians_to_degrees ( atan ( y / x ) );
    }
    else
    {
      temp1 = r8_huge ( );
    }

    temp2 = radians_to_degrees ( atan2 ( y, x ) );
    temp3 = radians_to_degrees ( atan4 ( y, x ) );

    cout << "  " << setw(12) << x
         << "  " << setw(12) << y
         << "  " << setw(12) << temp1
         << "  " << setw(12) << temp2
         << "  " << setw(12) << temp3<< "\n";
  }

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

void test006 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST006 tests BALL_UNIT_SAMPLE_2D.
//
//  Modified:
//
//    07 August 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double average[DIM_NUM];
  double average_r;
  double average_theta;
  int i;
  int j;
  int sample_num = 1000;
  int seed = 123456789;
  double temp;
  double theta;
  double *x;

  cout << "\n";
  cout << "TEST006\n";
  cout << "  For the unit ball in 2 dimensions (the disk):\n";
  cout << "  BALL_UNIT_SAMPLE_2D samples;\n";

  cout << "\n";
  cout << "  A few sample values:\n";
  cout << "\n";

  for ( i = 1; i <= 5; i++ )
  {
    x = ball_unit_sample_2d ( &seed );
    for ( j = 0; j < DIM_NUM; j++ )
    {
      cout<< "  " << setw(10) << x[j];
    }
    cout << "\n";
    delete [] x;
  }

  cout << "\n";
  cout << "  Number of sample points = " << sample_num << "\n";

  r8vec_zero ( DIM_NUM, average );

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_2d ( &seed );
    for ( j = 0; j < DIM_NUM; j++ )
    {
      average[j] = average[j] + x[j];
    }
    delete [] x;
  }

  for ( j = 0; j < DIM_NUM; j++ )
  {
    average[j] = average[j] / ( double ) sample_num;
  }

  cout << "\n";
  cout << "  Now average the points, which should get a value\n";
  cout << "  close to zero, and closer as N increases.\n";
  cout << "\n";
  cout << "  Average:        ";
  for ( j = 0; j < DIM_NUM; j++ )
  {
    cout << "  " << setw(10) << average[j];
  }
  cout << "\n";

  average_r = 0.0;

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_2d ( &seed );
    temp = 0.0;
    for ( j = 0; j < DIM_NUM; j++ )
    {
      temp = temp + x[j] * x[j];
    }
    average_r = average_r + sqrt ( temp );
    delete [] x;
  }

  average_r = average_r / ( double ) sample_num;

  cout << "\n";
  cout << "  Now average the distance of the points from the center,\n";
  cout << "  which should be 1/sqrt(2) = "
    << 1.0 / sqrt ( 2.0 ) << "\n";
  cout << "\n";
  cout << "  Average:        " << average_r << "\n";

  average_theta = 0.0;

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_2d ( &seed );
    theta = atan4 ( x[1], x[0] );
    average_theta = average_theta + theta;
    delete [] x;
  }

  average_theta = average_theta / ( double ) sample_num;

  cout << "\n";
  cout << "  Now average the angle THETA,\n";
  cout << "  which should be PI.\n";
  cout << "\n";
  cout << "  Average:        " << average_theta << "\n";

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

void test007 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST007 tests BALL_UNIT_SAMPLE_3D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  double average[DIM_NUM];
  double average_phi;
  double average_r;
  double average_theta;
  int i;
  int j;
  int sample_num = 1000;
  double phi;
  double r;
  int seed = 123456789;
  double theta;
  double *x;

  cout << "\n";
  cout << "TEST007\n";
  cout << "  For the unit ball in 3 dimensions:\n";
  cout << "  BALL_UNIT_SAMPLE_3D samples;\n";

  cout << "\n";
  cout << "  A few sample values:\n";
  cout << "\n";

  for ( i = 1; i <= 5; i++ )
  {
    x = ball_unit_sample_3d ( &seed );
    for ( j = 0; j < DIM_NUM; j++ )
    {
      cout << "  " << setw(8) << x[j];
    }
    cout << "\n";
    delete [] x;
  }

  cout << "\n";
  cout << "  Number of sample points = " << sample_num << "\n";

  r8vec_zero ( DIM_NUM, average );

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_3d ( &seed );
    for ( j = 0; j < DIM_NUM; j++ )
    {
      average[j] = average[j] + x[j];
    }
    delete [] x;
  }

  for ( j = 0; j < DIM_NUM; j++ )
  {
    average[j] = average[j] / ( double ) ( sample_num );
  }

  cout << "\n";
  cout << "  Now average the points, which should get a value\n";
  cout << "  close to zero, and closer as sample_num increases.\n";
  cout << "\n";
  cout << "  Average:";
  for ( i = 0; i < DIM_NUM; i++ )
  {
    cout << "  " << setw(8) << average[i];
  }
  cout << "\n";

  seed = 123456789;

  average_r = 0.0;

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_3d ( &seed );
    r = 0.0;
    for ( j = 0; j < DIM_NUM; j++ )
    {
      r = r + x[j] * x[j];
    }
    r = sqrt ( r );
    average_r = average_r + r;
    delete [] x;
  }

  average_r = average_r / ( double ) ( sample_num );

  cout << "\n";
  cout << "  Now average the distance of the points from\n";
  cout << "  the center, which should be the \n";
  cout << "  1/2^(1/n) = " <<  pow ( 0.5, 1.0 / ( double ) ( DIM_NUM ) ) << "\n";
  cout << "\n";
  cout << "  Average:        " << average_r << "\n";

  average_theta = 0.0;

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_3d ( &seed );
    theta = atan4 ( x[1], x[0] );
    average_theta = average_theta + theta;
    delete [] x;
  }

  average_theta = average_theta / ( double ) ( sample_num );

  cout << "\n";
  cout << "  Now average the angle THETA,\n";
  cout << "  which should be PI.\n";
  cout << "\n";
  cout << "  Average:        " << average_theta << "\n";

  average_phi = 0.0;

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_3d ( &seed );
    r = 0.0;
    for ( j = 0; j < DIM_NUM; j++ )
    {
      r = r + x[j] * x[j];
    }
    r = sqrt ( r );
    phi = acos ( x[2] / r );
    average_phi = average_phi + phi;
    delete [] x;
  }

  average_phi = average_phi / ( double ) ( sample_num );

  cout << "\n";
  cout << "  Now average the angle PHI,\n";
  cout << "  which should be PI/2.\n";
  cout << "\n";
  cout << "  Average:        " << average_phi << "\n";

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

void test0075 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0075 tests BALL_UNIT_SAMPLE_ND.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  double average[DIM_NUM];
  double average_phi;
  double average_r;
  double average_theta;
  int i;
  int j;
  double pi = 3.141592653589793;
  int sample_num = 1000;
  double phi;
  double r;
  int seed = 123456789;
  double theta;
  double *x;

  cout << "\n";
  cout << "TEST0075\n";
  cout << "  For the unit ball in N dimensions:\n";
  cout << "  BALL_UNIT_SAMPLE_ND samples;\n";

  cout << "\n";
  cout << "  A few sample values:\n";
  cout << "\n";

  for ( i = 1; i <= 5; i++ )
  {
    x = ball_unit_sample_nd ( DIM_NUM, &seed );
    for ( j = 0; j < DIM_NUM; j++ )
    {
      cout << "  " << setw(8) << x[j];
    }
    cout << "\n";
    delete [] x;
  }

  cout << "\n";
  cout << "  Number of sample points = " << sample_num << "\n";

  r8vec_zero ( DIM_NUM, average );

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_nd ( DIM_NUM, &seed );
    for ( j = 0; j < DIM_NUM; j++ )
    {
      average[j] = average[j] + x[j];
    }
    delete [] x;
  }

  for ( j = 0; j < DIM_NUM; j++ )
  {
    average[j] = average[j] / ( double ) ( sample_num );
  }

  cout << "\n";
  cout << "  Now average the points, which should get a value\n";
  cout << "  close to zero, and closer as N increases.\n";
  cout << "\n";
  cout << "  Average:        ";
  for ( j = 0; j < DIM_NUM; j++ )
  {
    cout << "  " << setw(8) << average[j];
  }
  cout << "\n";

  seed = 123456789;

  average_r = 0.0;

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_nd ( DIM_NUM, &seed );
    r = 0.0;
    for ( j = 0; j < DIM_NUM; j++ )
    {
      r = r + x[j] * x[j];
    }
    r = sqrt ( r );
    average_r = average_r + r;
    delete [] x;
  }

  average_r = average_r / ( double ) ( sample_num );

  cout << "\n";
  cout << "  Now average the distance of the points from\n";
  cout << "  the center, which should be the\n";
  cout << "  1/2^(1/dim_num) = " 
       << pow ( 0.5, 1.0 / ( double ) ( DIM_NUM ) ) << "\n";
  cout << "\n";
  cout << "  Average:        " << average_r << "\n";

  average_theta = 0.0;

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_nd ( DIM_NUM, &seed );
    theta = atan4 ( x[1], x[0] );
    average_theta = average_theta + theta;
    delete [] x;
  }

  average_theta = average_theta / ( double ) ( sample_num );

  cout << "\n";
  cout << "  Now average the angle THETA,\n";
  cout << "  which should be PI.\n";
  cout << "\n";
  cout << "  Average:        " << average_theta << "\n";

  average_phi = 0.0;

  for ( i = 1; i <= sample_num; i++ )
  {
    x = ball_unit_sample_nd ( DIM_NUM, &seed );
    r = 0.0;
    for ( j = 0; j < DIM_NUM; j++ )
    {
      r = r + x[j] * x[j];
    }
    r = sqrt ( r );
    phi = acos ( x[2] / r );
    average_phi = average_phi + phi;
    delete [] x;
  }

  average_phi = average_phi / ( double ) ( sample_num );

  cout << "\n";
  cout << "  Now average the angle PHI,\n";
  cout << "  which should be PI/2.\n";
  cout << "\n";
  cout << "  Average:        " << average_phi << "\n";

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

void test008 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST008 tests BASIS_MAP_3D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
  double *a = NULL;
  double c[3*3];
  int i;
  int j;
  int k;
  double u[3*3] = {
    1.0, 2.0, 3.0,
    0.0, 0.0, 1.0,
    1.0, 0.0, 2.0 };
  double v[3*3] = {
    14.0, 4.0, 4.0,
     3.0, 1.0, 0.0,
     7.0, 3.0, 2.0 };

  cout << "\n";
  cout << "TEST008\n";
  cout << "  BASIS_MAP_3D computes the linear transform A\n";
  cout << "  which maps vectors U1, U2 and U3 to vectors\n";
  cout << "  V1, V2 and V3.\n";

  r8mat_print ( 3, 3, u, "  The matrix U" );

  r8mat_print ( 3, 3, v, "  The matrix V" );

  a = basis_map_3d ( u, v );

  if ( a == NULL )
  {
    cout << "\n";
    cout << "  The matrix [ U1 | U2 | U3 ] was singular.\n";
    cout << "  No transformation was computed.\n";
    return;
  }

  r8mat_print ( 3, 3, a, "  The transformation matrix" );

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

  r8mat_print ( 3, 3, c, "  The product matrix A * [ U1 | U2 | U3 ]" );

  delete [] a;
  return;
}
//****************************************************************************80********80

void test0085 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0085 tests BOX_01_CONTAINS_POINT_2D.
//
//  Modified:
//
//    08 July 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM
# define N 46

  int i;
  int j;
  double p[2];
  double xhi =  1.2;
  double xlo = -0.3;
  double yhi =  1.4;
  double ylo = -0.1;

  cout << "\n";
  cout << "TEST0085\n";
  cout << "  BOX_01_CONTAINS_POINT_2D reports if the unit box\n";
  cout << "  contains a point.\n";
  cout << "\n";
  cout << "  We will call the function repeatedly, and draw\n";
  cout << "  a sketch of the unit square.\n";
  cout << "\n";

  for ( i = 1; i <= N; i++ )
  {
    p[1] = ( ( double ) ( N - i     ) * yhi
           + ( double ) (     i - 1 ) * ylo )
           / ( double ) ( N     - 1 );

    cout << "  ";
    for ( j = 1; j <= N; j++ )
    {
      p[0] = ( ( double ) ( N - j     ) * xlo
             + ( double ) (     j - 1 ) * xhi )
             / ( double ) ( N     - 1 );

      if ( box_01_contains_point_2d ( p ) )
      {
        cout << '*';
      }
      else
      {
        cout << '-';
      }
    }
    cout << "\n";
  }

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

void test0087 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0087 tests BOX_CONTAINS_POINT_2D.
//
//  Modified:
//
//    08 July 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define N 46

  int i;
  int j;
  double p[2];
  double p1[2] = { -0.1, 0.3 };
  double p2[2] = {  1.1, 0.9 };
  double xhi =  1.2;
  double xlo = -0.3;
  double yhi =  1.4;
  double ylo = -0.1;

  cout << "\n";
  cout << "TEST0087\n";
  cout << "  BOX_CONTAINS_POINT_2D reports if a box\n";
  cout << "  contains a point.\n";
  cout << "\n";
  cout << "  We will call the function repeatedly, and draw\n";
  cout << "  a sketch of the box.\n";
  cout << "\n";

  for ( i = 1; i <= N; i++ )
  {
    p[1] = ( ( double ) ( N - i     ) * yhi
           + ( double ) (     i - 1 ) * ylo )
           / ( double ) ( N     - 1 );

    cout << "  ";
    for ( j = 1; j <= N; j++ )
    {
      p[0] = ( ( double ) ( N - j     ) * xlo
             + ( double ) (     j - 1 ) * xhi )
             / ( double ) ( N     - 1 );

      if ( box_contains_point_2d ( p1, p2, p ) )
      {
        cout << '*';
      }
      else
      {
        cout << '-';
      }
    }
    cout << "\n";
  }

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

void test009 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST009 tests BOX_SEGMENT_CLIP_2D.
//
//  Modified:
//
//    12 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define TEST_NUM 5

  int i;
  int ival;
  double p1[DIM_NUM] = { -10.0, 10.0 };
  double p2[DIM_NUM] = {  10.0, 20.0 };
  double pa[DIM_NUM];
  double pb[DIM_NUM];
  double qa[DIM_NUM];
  double qb[DIM_NUM];
  double p1_test[DIM_NUM*TEST_NUM] = {
     1.0,  2.0,
    -3.0, 12.0,
   -20.0, 20.0,
   -20.0, 40.0,
    10.0, 40.0 };
  double p2_test[DIM_NUM*TEST_NUM] = {
     8.0, 16.0,
     5.0, 12.0,
     7.0, 20.0,
     0.0,  0.0,
    20.0, 30.0 };
  int test;

  cout << "\n";
  cout << "TEST009\n";
  cout << "  BOX_SEGMENT_CLIP_2D clips a line with respect to a box.\n";
  cout << "\n";
  cout << "  The lower left box corner is:\n";
  cout << "\n";
  for ( i = 0; i < DIM_NUM; i++ )
  {
    cout << "  " << setw(8) << p1[i];
  }
  cout << "\n";
  cout << "\n";
  cout << "  The upper right box corner is:\n";
  cout << "\n";
  for ( i = 0; i < DIM_NUM; i++ )
  {
    cout << "  " << setw(8) << p2[i];
  }
  cout << "\n";
  cout << "\n";
  cout << "  We list the points PA and PB, and then\n";
  cout << "  the clipped values.\n";
  cout << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    r8vec_copy ( DIM_NUM, p1_test+test*DIM_NUM, pa );
    r8vec_copy ( DIM_NUM, p2_test+test*DIM_NUM, pb );
    r8vec_copy ( DIM_NUM, pa, qa );
    r8vec_copy ( DIM_NUM, pb, qb );

    ival = box_segment_clip_2d ( p1, p2, qa, qb );

    cout << "\n";
    for ( i = 0; i < DIM_NUM; i++ )
    {
      cout << "  " << setw(8) << pa[i];
    }
    for ( i = 0; i < DIM_NUM; i++ )
    {
      cout << "  " << setw(8) << pb[i];
    }
    cout << "\n";
    cout << "\n";

    if ( ival == -1 )
    {
      cout << "  Line is outside the box.\n";
    }
    else if ( ival == 0 )
    {
      cout << "  Line is inside the box.\n";
    }
    else if ( ival == 1 )
    {
      for ( i = 0; i < DIM_NUM; i++ )
      {
        cout << "  " << setw(8) << qa[i];
      }
    }
    else if ( ival == 2 )
    {
      for ( i = 0; i < DIM_NUM; i++ )
      {
        cout << "  " << "        ";
      }
      for ( i = 0; i < DIM_NUM; i++ )
      {
        cout << "  " << setw(8) << qb[i];
      }
      cout << "\n";
    }
    else if ( ival == 3 )
      {
      for ( i = 0; i < DIM_NUM; i++ )
      {
        cout << "  " << setw(8) << qa[i];
      }
      for ( i = 0; i < DIM_NUM; i++ )
      {
        cout << "  " << setw(8) << qb[i];
      }
      cout << "\n";
    }

  }

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

void test010 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST010 tests BOX_RAY_INT_2D.
//
//  Modified:
//
//    08 July 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define TEST_NUM 3

  double p1[DIM_NUM] = { 0.0, 0.0 };
  double p2[DIM_NUM] = { 5.0, 3.0 };
  double pa[DIM_NUM];
  double pa_test[DIM_NUM*TEST_NUM] = {
    3.0, 1.0,
    4.0, 1.0,
    3.0, 1.0 };
  double pb[DIM_NUM];
  double pb_test[DIM_NUM*TEST_NUM] = {
    5.0, 5.0,
    3.0, 1.0,
    4.0, 2.0 };
  double pc_test[DIM_NUM*TEST_NUM] = {
    4.0, 3.0,
    0.0, 1.0,
    5.0, 3.0 };
  double pint[DIM_NUM];
  int test;

  cout << "\n";
  cout << "TEST010\n";
  cout << "  For a box with coordinate line sides in 2D,\n";
  cout << "  BOX_RAY_INT_2D computes the intersection of\n";
  cout << "    a shape and a ray whose origin is within\n";
  cout << "    the shape.\n";
  cout << "\n";
  cout << "  Lower left box corner:\n";
  cout << "\n";
  cout << "  " << setw(12) << p1[0]
       << "  " << setw(12) << p1[1] << "\n";
  cout << "\n";
  cout << "  Upper right box corner:\n";
  cout << "\n";
  cout << "  " << setw(12) << p2[0]
       << "  " << setw(12) << p2[1] << "\n";
  cout << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    r8vec_copy ( DIM_NUM, pa_test+test*DIM_NUM, pa );
    r8vec_copy ( DIM_NUM, pb_test+test*DIM_NUM, pb );

    box_ray_int_2d ( p1, p2, pa, pb, pint );

    cout << "\n";
    cout << "  Origin:       " << setw(12) << pa[0]
         << "  " << setw(12) << pa[1] << "\n";
    cout << "  Point 2:      " << setw(12) << pb[0]
         << "  " << setw(12) << pb[1] << "\n";
    cout << "  Intersection: " << setw(12) << pint[0]
         << "  " << setw(12) << pint[1] << "\n";
    cout << "  Correct:      " << setw(12) << pc_test[0+test*DIM_NUM]
         << "  " << setw(12) << pc_test[1+test*DIM_NUM] << "\n";
  }

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

void test011 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST011 tests CIRCLE_DIA2IMP_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double pc[DIM_NUM];
  double p1[DIM_NUM];
  double p2[DIM_NUM];
  double r;
  double theta;

  cout << "\n";
  cout << "TEST011\n";
  cout << "  CIRCLE_DIA2IMP_2D converts a diameter to an\n";
  cout << "  implicit circle in 2D.\n";

  theta = 2.0;

  p1[0] = 2.0 + 5.0 * cos ( theta );
  p1[1] = 2.0 + 5.0 * sin ( theta );
  p2[0] = 2.0 - 5.0 * cos ( theta );
  p2[1] = 2.0 - 5.0 * sin ( theta );

  r8vec_print ( DIM_NUM, p1, "  P1:" );
  r8vec_print ( DIM_NUM, p2, "  P2:" );

  circle_dia2imp_2d ( p1, p2, &r, pc );

  circle_imp_print_2d ( r, pc, "  The implicit circle:" );

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

void test012 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST012 tests CIRCLE_LUNE_AREA_2D, CIRCLE_SECTOR_AREA_2D, CIRCLE_TRIANGLE_AREA_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double area1;
  double area2;
  double area3;
  double pc[DIM_NUM] = { 0.0, 0.0 };
  double pi = 3.141592653589793;
  double r = 1.0;
  int test;
  int test_num = 12;
  double theta1;
  double theta2;

  cout << "\n";
  cout << "TEST012\n";
  cout << "  CIRCLE_LUNE_AREA_2D computes the area of a\n";
  cout << "    circular lune, defined by joining the endpoints\n";
  cout << "    of a circular arc.\n";
  cout << "  CIRCLE_SECTOR_AREA_2D computes the area of a\n";
  cout << "    circular sector, defined by joining the endpoints\n";
  cout << "    of a circular arc to the center.\n";
  cout << "  CIRCLE_TRIANGLE_AREA_2D computes the signed area of a\n";
  cout << "    triangle, defined by joining the endpoints\n";
  cout << "    of a circular arc and the center.\n";
  cout << "\n";
  cout << "      R            Theta1      Theta2        Sector       Triangle     Lune\n";
  cout << "\n";

  for ( test = 0; test <= test_num; test++ )
  {
    theta1 = 0.0;
    theta2 = ( double ) ( test ) * 2.0 * pi / ( double ) ( test_num );

    area1 = circle_sector_area_2d ( r, pc, theta1, theta2 );

    area2 = circle_triangle_area_2d ( r, pc, theta1, theta2 );

    area3 = circle_lune_area_2d ( r, pc, theta1, theta2 );

    cout << "  " << setw(6)  << r
         << "  " << setw(12) << theta1
         << "  " << setw(12) << theta2
         << "  " << setw(12) << area1
         << "  " << setw(12) << area2
         << "  " << setw(12) << area3 << "\n";
  }

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

void test0125 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0125 tests CIRCLE_LUNE_AREA_2D and SPHERE_CAP_VOLUME_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double area;
  double h;
  double haver_sine;
  double pc[DIM_NUM] = { 0.0, 0.0 };
  double pi = 3.141592653589793;
  double r = 1.0;
  int test;
  int test_num = 12;
  double theta1;
  double theta2;
  double volume1;
  double volume2;

  cout << "\n";
  cout << "TEST0125\n";
  cout << "  CIRCLE_LUNE_AREA_2D computes the area of a\n";
  cout << "    circular lune, defined by joining the endpoints\n";
  cout << "    of a circular arc (THETA1,THETA2).\n";
  cout << "  SPHERE_CAP_VOLUME_2D computes the volume (area) of a\n";
  cout << "    spherical cap, defined by a plane that cuts the\n";
  cout << "    sphere to a thickness of H units.\n";
  cout << "  SPHERE_CAP_VOLUME_ND does the same operation,\n";
  cout << "    but in N dimensions.\n";
  cout << "\n";
  cout << "  The two routines should get the same results\n";
  cout << "  if THETA1, THETA2 and H correspond.\n";
  cout << "\n";
  cout << "  Using a radius R = " << r << "\n";

  cout << "\n";
  cout << "        Theta1      Theta2      H           Lune        Cap        Cap\n";
  cout << "                                            area        vol_3d     vol_nd\n";
  cout << "\n";

  for ( test = 1; test <= test_num; test++ )
  {
    h = 2.0 * r * ( double ) ( test ) / ( double ) ( test_num );

    haver_sine = sqrt ( r * r - ( r - h ) * ( r - h ) );

    if ( h <= r )
    {
      theta2 = arc_sine ( haver_sine / r );
    }
    else
    {
      theta2 = ( pi - arc_sine ( haver_sine / r ) );
    }

    theta1 = -theta2;

    area = circle_lune_area_2d ( r, pc, theta1, theta2 );

    volume1 = sphere_cap_volume_2d ( r, h );

    volume2 = sphere_cap_volume_nd ( DIM_NUM, r, h );

    cout << "  " << setw(10) << theta1
         << "  " << setw(10) << theta2
         << "  " << setw(10) << h
         << "  " << setw(10) << area
         << "  " << setw(10) << volume1
         << "  " << setw(10) << volume2 << "\n";
  }

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

void test0126 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0126 tests SPHERE_CAP_VOLUME_3D and SPHERE_CAP_VOLUME_ND.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  double h;
  int test;
  int test_num = 12;
  double r = 1.0;
  double volume1;
  double volume2;

  cout << "\n";
  cout << "TEST0126\n";
  cout << "  SPHERE_CAP_VOLUME_3D computes the volume of a\n";
  cout << "    spherical cap, defined by a plane that cuts the\n";
  cout << "    sphere to a thickness of H units.\n";
  cout << "  SPHERE_CAP_VOLUME_ND does the same operation,\n";
  cout << "    but in N dimensions.\n";
  cout << "\n";
  cout << "  Using a radius R = " << r << "\n";

  cout << "\n";
  cout << "        H           Cap        Cap\n";
  cout << "               volume_3d  volume_nd\n";
  cout << "\n";

  for ( test = 0; test <= test_num; test++ )
  {
    h = 2.0 * r * ( double ) ( test ) / ( double ) ( test_num );

    volume1 = sphere_cap_volume_3d ( r, h );

    volume2 = sphere_cap_volume_nd ( DIM_NUM, r, h );

    cout << "  " << setw(12) << h
         << "  " << setw(12) << volume1
         << "  " << setw(12) << volume2 << "\n";
  }

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

void test0127 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0127 tests SPHERE_CAP_AREA_3D and SPHERE_CAP_AREA_ND.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  double area1;
  double area2;
  double h;
  int test;
  int test_num = 12;
  double r = 1.0;

  cout << "\n";
  cout << "TEST0127\n";
  cout << "  SPHERE_CAP_AREA_3D computes the volume of a\n";
  cout << "    3D spherical cap, defined by a plane that cuts the\n";
  cout << "    sphere to a thickness of H units.\n";
  cout << "  SPHERE_CAP_AREA_ND computes the volume of an\n";
  cout << "    ND spherical cap, defined by a plane that cuts the\n";
  cout << "    sphere to a thickness of H units.\n";
  cout << "\n";
  cout << "        R           H           Cap         Cap\n";
  cout << "                                area_3d     area_nd\n";
  cout << "\n";

  for ( test = 0; test <= test_num; test++ )
  {
    h = 2.0 * r * ( double ) ( test ) / ( double ) ( test_num );

    area1 = sphere_cap_area_3d ( r, h );

    area2 = sphere_cap_area_nd ( DIM_NUM, r, h );

    cout << "  " << setw(12) << r
         << "  " << setw(12) << h
         << "  " << setw(12) << area1
         << "  " << setw(12) << area2 << "\n";
  }

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

void test013 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST013 tests CIRCLE_LUNE_CENTROID_2D and CIRCLE_SECTOR_CENTROID_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double *centroid1;
  double *centroid2;
  double pc[DIM_NUM] = { 5.0, 3.0 };
  double pi = 3.141592653589793;
  double r = 2.0;
  int test;
  int test_num = 12;
  double theta1 = 0.0;
  double theta2;

  cout << "\n";
  cout << "TEST013\n";
  cout << "  CIRCLE_LUNE_CENTROID_2D computes the centroid of a\n";
  cout << "    circular lune, defined by joining the endpoints\n";
  cout << "    of a circular arc.\n";
  cout << "  CIRCLE_SECTOR_CENTROID_2D computes the centroid of a\n";
  cout << "    circular sector, defined by joining the endpoints\n";
  cout << "    of a circular arc to the center.\n";

  circle_imp_print_2d ( r, pc, "  The implicit circle:" );

  cout << "\n";
  cout << "  The first angle of our lune and sector is always 0.\n";
  cout << "\n";
  cout << "                         Lune                       Sector\n";
  cout << "  THETA2           X             Y             X             Y\n";
  cout << "\n";

  for ( test = 0; test <= test_num; test++ )
  {
    theta2 = ( double ) ( test ) * 2.0 * pi / ( double ) ( test_num );

    centroid1 = circle_lune_centroid_2d ( r, pc, theta1, theta2 );

    centroid2 = circle_sector_centroid_2d ( r, pc, theta1, theta2 );

    cout << "  " << setw(12) << theta2
         << "  " << setw(12) << centroid1[0]
         << "  " << setw(12) << centroid1[1]
         << "  " << setw(12) << centroid2[0]
         << "  " << setw(12) << centroid2[1] << "\n";

    delete [] centroid1;
    delete [] centroid2;
  }

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

void test014 ( void )

//****************************************************************************80
//
//  Purpose:
//
//    TEST014 tests CIRCLE_EXP_CONTAINS_POINT_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM

  int inside;
  double p1[2];
  double p2[2];
  double p3[2];
  double p4[2];

  cout << "\n";
  cout << "TEST014\n";
  cout << "  CIRCLE_EXP_CONTAINS_POINT_2D determines if a\n";
  cout << "    point lies inside a circle.\n";
  cout << "\n";
  cout << "  Possible return values are:\n";
  cout << "\n";
  cout << "  -1: The point is inside the circle.\n";
  cout << "   0: The point is on the circle.\n";
  cout << "   1: The point is outside the circle\n";
  cout << "   2: Colinear data, the point is on the line.\n";
  cout << "   3: Colinear data, the point is not on the line.\n";
  cout << "   4: Two equal data points, the point is on the line.\n";
  cout << "   5: Two equal data points, the point is not on the line.\n";
  cout << "   6: All data points equal, the point is equal.\n";
  cout << "   7: All data points equal, the point is not equal.\n";
//
//  This point is inside.
//
  p1[0] = 4.0;
  p1[1] = 2.0;

  p2[0] = 1.0;
  p2[1] = 5.0;

  p3[0] = -2.0;
  p3[1] = 2.0;

  p4[0] = 2.0;
  p4[1] = 3.0;

  inside = circle_exp_contains_point_2d ( p1, p2, p3, p4 );

  cout << "\n";
  cout << "  P1 = " << p1[0] << "  " << p1[1] << "\n";
  cout << "  P2 = " << p2[0] << "  " << p2[1] << "\n";
  cout << "  P3 = " << p3[0] << "  " << p3[1] << "\n";
  cout << "  P4 = " << p4[0] << "  " << p4[1] << "\n";
  cout << "  INSIDE = " << inside << "\n";
//
//  This point is actually right on the circle.
//
  p1[0] = 4.0;
  p1[1] = 2.0;

  p2[0] = 1.0;
  p2[1] = 5.0;

  p3[0] = -2.0;
  p3[1] = 2.0;

  p4[0] = 1.0;
  p4[1] = -1.0;

  inside = circle_exp_contains_point_2d ( p1, p2, p3, p4 );

  cout << "\n";
  cout << "  P1 = " << p1[0] << "  " << p1[1] << "\n";
  cout << "  P2 = " << p2[0] << "  " << p2[1] << "\n";
  cout << "  P3 = " << p3[0] << "  " << p3[1] << "\n";
  cout << "  P4 = " << p4[0] << "  " << p4[1] << "\n";
  cout << "  INSIDE = " << inside << "\n";
//
//  This point is outside.
//
  p1[0] = 4.0;
  p1[1] = 2.0;

  p2[0] = 1.0;
  p2[1] = 5.0;

  p3[0] = -2.0;
  p3[1] = 2.0;

  p4[0] = 4.0;
  p4[1] = 6.0;

  inside = circle_exp_contains_point_2d ( p1, p2, p3, p4 );

  cout << "\n";
  cout << "  P1 = " << p1[0] << "  " << p1[1] << "\n";
  cout << "  P2 = " << p2[0] << "  " << p2[1] << "\n";
  cout << "  P3 = " << p3[0] << "  " << p3[1] << "\n";
  cout << "  P4 = " << p4[0] << "  " << p4[1] << "\n";
  cout << "  INSIDE = " << inside << "\n";

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

void test015 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST015 tests CIRCLE_EXP2IMP_2D and TRIANGLE_CIRCUMCIRCLE_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define TEST_NUM 3

  double *p1;
  double *p2;
  double *p3;
  double pc[DIM_NUM];
  double r;
  double *t;
  int test;
  double test_t[DIM_NUM*3*TEST_NUM] = {
    4.0, 2.0,
    1.0, 5.0,
   -2.0, 2.0,
    4.0, 2.0,
    5.0, 4.0,
    6.0, 6.0,
    4.0, 2.0,
    1.0, 5.0,
    4.0, 2.0 };

  cout << "\n";
  cout << "TEST015\n";
  cout << "  CIRCLE_EXP2IMP_2D computes the radius and\n";
  cout << "    center of the circle through three points.\n";
  cout << "  TRIANGLE_CIRCUMCIRCLE_2D computes the radius and\n";
  cout << "    center of the circle through the vertices of\n";
  cout << "    a triangle.\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    t = test_t+test*DIM_NUM*3;

    r8mat_transpose_print ( DIM_NUM, 3, t, "  The triangle:" );

    p1 = t;
    p2 = p1 + DIM_NUM;
    p3 = p2 + DIM_NUM;

    circle_exp2imp_2d ( p1, p2, p3, &r, pc );

    circle_imp_print_2d ( r, pc, "  The implicit circle:" );

    triangle_circumcircle_2d ( t, &r, pc );

    circle_imp_print_2d ( r, pc, "  The triangle's circumcircle:" );
  }

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

void test0155 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0155 tests CIRCLE_EXP2IMP_2D.
//
//  Modified:
//
//    10 March 2006
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define TEST_NUM 13

  double curvature;
  double pc[DIM_NUM];
  double p1[DIM_NUM] = { 0.0, 0.0 };
  double p2[DIM_NUM] = { 1.0, 0.0 };
  double p3[DIM_NUM];
  double pi = 3.141592653589793;
  double r;
  double theta;
  double theta_degrees;
  int test;

  cout << "\n";
  cout << "TEST0155\n";
  cout << "  CIRCLE_EXP2IMP_2D computes the radius and\n";
  cout << "    center of the circle through three points.\n";
  cout << "\n";
  cout << "  We can use this routine to compute, for three\n";
  cout << "  points in space, the circle incident to those\n";
  cout << "  points, and hence the radius of that circle,\n";
  cout << "  and hence the \"curvature\" of those points.\n";

  cout << "\n";
  cout << "  Our three points are:\n";
  cout << "\n";
  cout << "    (0,0)\n";
  cout << "    (1,0)\n";
  cout << "    (C,S)\n";
  cout << "\n";
  cout << "  C = cosine ( theta), S = sine ( theta ).\n";
  cout << "\n";
  cout << "  Test  Theta  Curvature\n";
  cout << "\n";

  for ( test = 1; test <= TEST_NUM; test++ )
  {
    theta = 2.0 * pi * ( double ) ( test     - 1 )
                     / ( double ) ( TEST_NUM - 1 );

    theta_degrees = 360.0 * ( double ) ( test     - 1 )
                         /  ( double ) ( TEST_NUM - 1 );

    p3[0] = cos ( theta );
    p3[1] = sin ( theta );

    circle_exp2imp_2d ( p1, p2, p3, &r, pc );

    if ( 0.0 < r )
    {
      curvature = 1.0 / r;
    }
    else
    {
      curvature = 0.0;
    }
    cout << "  " << setw(4)  << test
         << "  " << setw(5)  << theta_degrees
         << "  " << setw(14) << curvature << "\n";
  }

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

void test0156 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0156 tests CIRCLE_EXP2IMP_2D and CIRCLE_IMP2EXP_2D.
//
//  Modified:
//
//    17 May 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double p1[DIM_NUM];
  double p2[DIM_NUM];
  double p3[DIM_NUM];
  double pc1[DIM_NUM];
  double pc2[DIM_NUM];
  double r1;
  double r2;

  cout << "\n";
  cout << "TEST0156\n";
  cout << "  CIRCLE_EXP2IMP_2D converts an explicit circle\n";
  cout << "    to an implicit circle.\n";
  cout << "  CIRCLE_IMP2EXP_2D converts an implicit circle\n";
  cout << "    to an explicit circle.\n";

  pc1[0] = 10.0;
  pc1[1] =  5.0;
  r1 = 3.0;

  circle_imp_print_2d ( r1, pc1, "  The implicit circle:" );

  circle_imp2exp_2d ( r1, pc1, p1, p2, p3 );

  r8vec_print ( DIM_NUM, p1, "  P1:" );
  r8vec_print ( DIM_NUM, p2, "  P2:" );
  r8vec_print ( DIM_NUM, p3, "  P3:" );

  circle_exp2imp_2d ( p1, p2, p3, &r2, pc2 );

  circle_imp_print_2d ( r2, pc2, "  The recovered implicit circle:" );

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

void test016 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST016 tests CIRCLE_IMP_POINTS_2D and POLYGON_AREA_2D.
//
//  Modified:
//
//    10 July 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  int n;
  double *p;
  double pc[DIM_NUM] = { 5.0, -2.0 };
  double pi = 3.141592653589793;
  double r = 2.0;
  double result;

  cout << "\n";
  cout << "TEST016\n";
  cout << "  CIRCLE_IMP_POINTS_2D gets points on a circle;\n";
  cout << "  POLYGON_AREA_2D finds the area of a polygon.\n";

  circle_imp_print_2d ( r, pc, "  The implicit circle:" );

  cout << "\n";
  cout << "  The area = " << pi * r * r << "\n";

  n = 8;

  p = circle_imp_points_2d ( r, pc, n );

  r8mat_transpose_print ( DIM_NUM, n, p, "  Sample results:" );

  delete [] p;

  cout << "\n";
  cout << "  For any N, the sampled points define a polygon\n";
  cout << "  whose area approximates the circle area.\n";
  cout << "\n";
  cout << "  N      Area\n";
  cout << "\n";

  for ( n = 3; n <= 24; n++ )
  {
    p = circle_imp_points_2d ( r, pc, n );
    result = polygon_area_2d ( n, p );
    cout << "  " << setw(6) << n
         << "  " << setw(12) << result << "\n";

    delete [] p;
  }

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

void test0165 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0165 tests CIRCLE_IMP_POINTS_3D.
//
//  Modified:
//
//    10 March 2006
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3
# define N 12

  double nc[DIM_NUM] = { 1.0, 1.0, 1.0 };
  double *p;
  double pc[DIM_NUM] = { 5.0, -2.0, 1.0 };
  double r = 2.0;

  cout << "\n";
  cout << "TEST0165\n";
  cout << "  CIRCLE_IMP_POINTS_3D gets points on a circle in 3D;\n";

  circle_imp_print_3d ( r, pc, nc, "  The implicit circle:" );

  p = circle_imp_points_3d ( r, pc, nc, N );

  r8mat_transpose_print ( DIM_NUM, N, p, "  Points on the circle:" );

  delete [] p;

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

void test017 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST017 tests CIRCLE_IMP_POINTS_ARC_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define N 13
# define DIM_NUM 2

  double r = 2.0;
  double theta1;
  double theta2;
  double p[DIM_NUM*N];
  double pc[DIM_NUM] = { 5.0, -2.0 };
  double pi = 3.141592653589793;

  theta1 = pi / 2.0;
  theta2 = 3.0 * pi / 2.0;

  cout << "\n";
  cout << "TEST017\n";
  cout << "  CIRCLE_IMP_POINTS_ARC_2D returns points on a\n";
  cout << "  circular arc.\n";
  cout << "\n";
  cout << "  The circle will have center " << pc[0] << ", "<< pc[1] << "\n";
  cout << "  and radius R = " << r << "\n";
  cout << "\n";
  cout << "  The arc extends from THETA1 = " << theta1 << "\n";
  cout << "  to THETA2 = " << theta2 << "\n";

  circle_imp_points_arc_2d ( r, pc, theta1, theta2, N, p );

  r8mat_transpose_print ( DIM_NUM, N, p, "  Sample results:" );

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

void test018 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST018 tests CIRCLE_IMP_POINT_DIST_2D and CIRCLES_IMP_INT_2D;
//
//  Modified:
//
//    10 July 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define TEST_NUM 5

  double d1;
  double d2;
  int int_num;
  double pint[DIM_NUM*2];
  double pc1[DIM_NUM] = { 0.0, 0.0 };
  double *pc2;
  double pc2_test[DIM_NUM*TEST_NUM] = {
    5.0,       5.0,
    7.0710678, 7.0710678,
    4.0,       0.0,
    6.0,       0.0,
    0.0,       0.0 };
  double r1 = 5.0;
  double r2;
  double r2_test[TEST_NUM] = { 0.5, 5.0, 3.0, 3.0, 5.0 };
  int test;

  cout << "\n";
  cout << "TEST018\n";
  cout << "  CIRCLE_IMP_POINT_DIST_2D finds the\n";
  cout << "    distance from a point to a circle.\n";
  cout << "  CIRCLES_IMP_INT_2D determines the intersections of\n";
  cout << "    two circles in 2D.\n";

  circle_imp_print_2d ( r1, pc1, "  The first circle:" );

  for ( test = 0; test < TEST_NUM; test++ )
  {
    r2 = r2_test[test];
    pc2 = pc2_test + test*DIM_NUM;

    circle_imp_print_2d ( r2, pc2, "  The second circle:" );

    circles_imp_int_2d ( r1, pc1, r2, pc2, &int_num, pint );

    if ( int_num == 0 )
    {
      cout << "\n";
      cout << "  The circles do not intersect.\n";
    }
    else if ( int_num == 1 )
    {
      cout << "\n";
      cout << "  The circles intersect at one point:\n";
      cout << "\n";
      cout << "        P         Dist 1  Dist 2\n";
      cout << "\n";
      d1 = circle_imp_point_dist_2d ( r1, pc1, pint+0*DIM_NUM );
      d2 = circle_imp_point_dist_2d ( r2, pc2, pint+0*DIM_NUM );
      cout << "  " << setw(8) << pint[0+0*DIM_NUM]
           << "  " << setw(8) << pint[1+0*DIM_NUM]
           << "  " << setw(8) << d1
           << "  " << setw(8) << d2 << "\n";
    }
    else if ( int_num == 2 )
    {
      cout << "\n";
      cout << "  The circles intersect at two points:\n";
      cout << "\n";
      cout << "        P         Dist 1  Dist 2\n";
      cout << "\n";
      d1 = circle_imp_point_dist_2d ( r1, pc1, pint+0*DIM_NUM );
      d2 = circle_imp_point_dist_2d ( r2, pc2, pint+0*DIM_NUM );
      cout << "  " << setw(8) << pint[0+0*DIM_NUM]
           << "  " << setw(8) << pint[1+0*DIM_NUM]
           << "  " << setw(8) << d1
           << "  " << setw(8) << d2 << "\n";
      d1 = circle_imp_point_dist_2d ( r1, pc1, pint+1*DIM_NUM );
      d2 = circle_imp_point_dist_2d ( r2, pc2, pint+1*DIM_NUM );
      cout << "  " << setw(8) << pint[0+1*DIM_NUM]
           << "  " << setw(8) << pint[1+1*DIM_NUM]
           << "  " << setw(8) << d1
           << "  " << setw(8) << d2 << "\n";
    }
    else if ( int_num == 3 )
    {
      cout << "\n";
      cout << "  The circles coincide (infinite intersection).\n";
    }
  }

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

void test0183 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0183 tests CIRCLE_LLR2IMP_2D.
//
//  Modified:
//
//    16 November 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double d1;
  double d2;
  double d3;
  double d4;
  double p_hi =  10.0;
  double p_lo = -10.0;
  double *pc;
  double *p1;
  double *p2;
  double *q1;
  double *q2;
  double r;
  double r_hi;
  double r_lo;
  int seed;
  int test;
  int test_num = 5;

  cout << "\n";
  cout << "TEST0183\n";
  cout << "  CIRCLE_LLR2IMP_3D is given:\n";
  cout << "  a line through P1 and P2,\n";
  cout << "  a line through Q1 and Q2,\n";
  cout << "  and a radius R,\n";
  cout << "  and determines the centers C of 4 circles\n";
  cout << "  of the given radius, tangent to both lines.\n";

  seed = 123456789;

  for ( test = 1; test <= test_num; test++ )
  {
    p1 = r8vec_uniform ( DIM_NUM, p_lo, p_hi, &seed );
    p2 = r8vec_uniform ( DIM_NUM, p_lo, p_hi, &seed );
    q1 = r8vec_uniform ( DIM_NUM, p_lo, p_hi, &seed );
    q2 = r8vec_uniform ( DIM_NUM, p_lo, p_hi, &seed );

    r_lo = 1.0;
    r_hi = 5.0;
    r = r8_uniform ( r_lo, r_hi, &seed );

    cout << "\n";
    cout << "  Radius R = " << r << "\n";

    cout << "  Point #P1: ( "
         << setw(12) << p1[0] << ","
         << setw(12) << p1[1] << ")\n";
    cout << "  Point #P2: ( "
         << setw(12) << p2[0] << ","
         << setw(12) << p2[1] << ")\n";
    cout << "  Point #Q1: ( "
         << setw(12) << q1[0] << ","
         << setw(12) << q1[1] << ")\n";
    cout << "  Point #Q2: ( "
         << setw(12) << q2[0] << ","
         << setw(12) << q2[1] << ")\n";

    pc = circle_llr2imp_2d ( p1, p2, q1, q2, r );

    cout << "  Center #1: ( "
         << setw(12) << pc[0+0*DIM_NUM] << ","
         << setw(12) << pc[1+0*DIM_NUM] << ")\n";
    cout << "  Center #2: ( "
         << setw(12) << pc[0+1*DIM_NUM] << ","
         << setw(12) << pc[1+1*DIM_NUM] << ")\n";
    cout << "  Center #3: ( "
         << setw(12) << pc[0+2*DIM_NUM] << ","
         << setw(12) << pc[1+2*DIM_NUM] << ")\n";
    cout << "  Center #4: ( "
         << setw(12) << pc[0+3*DIM_NUM] << ","
         << setw(12) << pc[1+3*DIM_NUM] << ")\n";
//
//  Check that the lines are the right distance from the center.
//
    d1 = line_exp_point_dist_2d ( p1, p2, pc+0*2 );
    d2 = line_exp_point_dist_2d ( p1, p2, pc+1*2 );
    d3 = line_exp_point_dist_2d ( p1, p2, pc+2*2 );
    d4 = line_exp_point_dist_2d ( p1, p2, pc+3*2 );

    cout << "  " << setw(12) << d1
         << "  " << setw(12) << d2
         << "  " << setw(12) << d3
         << "  " << setw(12) << d4 << "\n";

    d1 = line_exp_point_dist_2d ( q1, q2, pc+0*2 );
    d2 = line_exp_point_dist_2d ( q1, q2, pc+1*2 );
    d3 = line_exp_point_dist_2d ( q1, q2, pc+2*2 );
    d4 = line_exp_point_dist_2d ( q1, q2, pc+3*2 );

    cout << "  " << setw(12) << d1
         << "  " << setw(12) << d2
         << "  " << setw(12) << d3
         << "  " << setw(12) << d4 << "\n";

    delete [] p1;
    delete [] p2;
    delete [] pc;
    delete [] q1;
    delete [] q2;
  }

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

void test0185 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0185 tests CIRCLE_PPPR2IMP_3D.
//
//  Modified:
//
//    12 November 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  double d11;
  double d12;
  double d21;
  double d22;
  int i;
  double normal[DIM_NUM];
  double p_hi =  10.0;
  double p_lo = -10.0;
  double pc[DIM_NUM*2];
  double *p1;
  double *p2;
  double *p3;
  double r;
  double r_hi;
  double r_lo;
  int seed;
  int test;
  int test_num = 5;

  cout << "\n";
  cout << "TEST0185\n";
  cout << "  CIRCLE_PPPR2IMP_3D is given 3D points P1, P2, P3,\n";
  cout << "  and a radius R,\n";
  cout << "  and determines the centers C of two circles\n";
  cout << "  of the given radius, passing through P1 and P2\n";
  cout << "  and lying in the plane of P1, P2 and P3.\n";

  seed = 123456789;

  for ( test = 1; test <= test_num; test++ )
  {
    p1 = r8vec_uniform ( DIM_NUM, p_lo, p_hi, &seed );
    p2 = r8vec_uniform ( DIM_NUM, p_lo, p_hi, &seed );
    p3 = r8vec_uniform ( DIM_NUM, p_lo, p_hi, &seed );

    r_lo = r8vec_distance ( DIM_NUM, p1, p2 );
    r_hi = r_lo + 5.0;
    r = r8_uniform ( r_lo, r_hi, &seed );

    cout << "\n";
    cout << "  Radius R = " << r << "\n";

    cout << "  Point #1: ( "
         << setw(12) << p1[0] << ","
         << setw(12) << p1[1] << ")\n";
    cout << "  Point #2: ( "
         << setw(12) << p2[0] << ","
         << setw(12) << p2[1] << ")\n";
    cout << "  Point #3: ( "
         << setw(12) << p3[0] << ","
         << setw(12) << p3[1] << ")\n";

    circle_pppr2imp_3d ( p1, p2, p3, r, pc, normal );

    cout << "  Center #1: ( "
         << setw(12) << pc[0+0*DIM_NUM] << ","
         << setw(12) << pc[1+0*DIM_NUM] << ")\n";
    cout << "  Center #2: ( "
         << setw(12) << pc[0+1*DIM_NUM] << ","
         << setw(12) << pc[1+1*DIM_NUM] << ")\n";
//
//  Check that the points are the right distance from the center.
//
    d11 = r8vec_distance ( DIM_NUM, p1, pc+0*DIM_NUM );
    d21 = r8vec_distance ( DIM_NUM, p2, pc+0*DIM_NUM );
    d12 = r8vec_distance ( DIM_NUM, p1, pc+1*DIM_NUM );
    d22 = r8vec_distance ( DIM_NUM, p2, pc+1*DIM_NUM );

    cout << "  " << setw(12) << d11
         << "  " << setw(12) << d21
         << "  " << setw(12) << d12
         << "  " << setw(12) << d22 << "\n";
//
//  Check that the radial vector to the point is perpendicular to NORMAL.
//
    d11 = 0.0;
    d21 = 0.0;
    d12 = 0.0;
    d22 = 0.0;
    for ( i = 0; i < DIM_NUM; i++ )
    {
      d11 = d11 + normal[i] * ( p1[i] - pc[i+0*DIM_NUM] );
      d21 = d21 + normal[i] * ( p2[i] - pc[i+0*DIM_NUM] );
      d12 = d12 + normal[i] * ( p1[i] - pc[i+1*DIM_NUM] );
      d22 = d22 + normal[i] * ( p2[i] - pc[i+1*DIM_NUM] );
    }

    cout << "  " << setw(12) << d11
         << "  " << setw(12) << d21
         << "  " << setw(12) << d12
         << "  " << setw(12) << d22 << "\n";
  }

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

void test019 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST019 tests CIRCLE_PPR2IMP_2D.
//
//  Modified:
//
//    11 November 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double d11;
  double d12;
  double d21;
  double d22;
  double p_hi =  10.0;
  double p_lo = -10.0;
  double *pc;
  double *p1;
  double *p2;
  double r;
  double r_hi;
  double r_lo;
  int seed;
  int test;
  int test_num = 5;

  cout << "\n";
  cout << "TEST019\n";
  cout << "  CIRCLE_PPR2IMP_2D is given points P1 and P2,\n";
  cout << "  and a radius R,\n";
  cout << "  and determines the centers C of two circles\n";
  cout << "  of the given radius, passing through P1 and P2.\n";

  seed = 123456789;

  for ( test = 1; test <= test_num; test++ )
  {
    p1 = r8vec_uniform ( DIM_NUM, p_lo, p_hi, &seed );
    p2 = r8vec_uniform ( DIM_NUM, p_lo, p_hi, &seed );

    r_lo = r8vec_distance ( DIM_NUM, p1, p2 );
    r_hi = r_lo + 5.0;
    r = r8_uniform ( r_lo, r_hi, &seed );

    cout << "\n";
    cout << "  Radius R = " << r << "\n";

    cout << "  Point #1: ( "
         << setw(12) << p1[0] << ","
         << setw(12) << p1[1] << ")\n";
    cout << "  Point #2: ( "
         << setw(12) << p2[0] << ","
         << setw(12) << p2[1] << ")\n";

    pc = circle_ppr2imp_2d ( p1, p2, r );

    cout << "  Center #1: ( "
         << setw(12) << pc[0+0*DIM_NUM] << ","
         << setw(12) << pc[1+0*DIM_NUM] << ")\n";
    cout << "  Center #2: ( "
         << setw(12) << pc[0+1*DIM_NUM] << ","
         << setw(12) << pc[1+1*DIM_NUM] << ")\n";

    d11 = r8vec_distance ( DIM_NUM, p1, pc+0*DIM_NUM );
    d21 = r8vec_distance ( DIM_NUM, p2, pc+0*DIM_NUM );
    d12 = r8vec_distance ( DIM_NUM, p1, pc+1*DIM_NUM );
    d22 = r8vec_distance ( DIM_NUM, p2, pc+1*DIM_NUM );

    cout << "  " << setw(12) << d11
         << "  " << setw(12) << d21
         << "  " << setw(12) << d12
         << "  " << setw(12) << d22 << "\n";

    delete [] p1;
    delete [] p2;
    delete [] pc;
  }

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

void test020 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST020 tests CUBE_SIZE_3D and CUBE_SHAPE_3D.
//
//  Modified:
//
//    22 July 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  int edge_num;
  int face_num;
  int *face_order;
  int face_order_max;
  int *face_point;
  int point_num;
  double *point_coord;

  cout << "\n";
  cout << "TEST020\n";
  cout << "  For the cube,\n";
  cout << "  CUBE_SIZE_3D returns dimension information;\n";
  cout << "  CUBE_SHAPE_3D returns face and order information.\n";
  cout << "  SHAPE_PRINT_3D prints this information.\n";
//
//  Get the sizes.
//
  cube_size_3d ( &point_num, &edge_num, &face_num, &face_order_max );

  cout << "\n";
  cout << "    Number of vertices: " << point_num << "\n";
  cout << "    Number of edges   : " << edge_num << "\n";
  cout << "    Number of faces   : " << face_num << "\n";
  cout << "    Maximum face order: " << face_order_max << "\n";
//
//  Make room for the data.
//
  face_order = new int[face_num];
  face_point = new int[face_order_max*face_num];
  point_coord = new double[DIM_NUM*point_num];
//
//  Get the data.
//
  cube_shape_3d ( point_num, face_num, face_order_max, point_coord,
    face_order, face_point );
//
//  Print the data.
//
  shape_print_3d ( point_num, face_num, face_order_max,
    point_coord, face_order, face_point );

  delete [] face_order;
  delete [] face_point;
  delete [] point_coord;

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

void test0201 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0201 tests CYLINDER_POINT_DIST_3D.
//
//  Modified:
//
//    23 August 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3
# define TEST_NUM 6

  double dist;
  double dist_test[TEST_NUM] = {
    3.0, 0.5, 5.0, 8.0, 1.0, 0.25 };
  double *p;
  double p_test[DIM_NUM*TEST_NUM] = {
      4.0,    0.5,  0.0,
     -0.5,   -1.0,  0.0,
      4.0,    6.0,  0.0,
      0.75, -10.0,  0.0,
      0.0,    0.0,  0.0,
      0.25,   1.75, 0.0 };
  double p1[DIM_NUM] = { 0.0, -2.0, 0.0 };
  double p2[DIM_NUM] = { 0.0,  2.0, 0.0 };
  double r = 1.0;
  int test;

  cout << "\n";
  cout << "TEST0201\n";
  cout << "  CYLINDER_POINT_DIST_3D computes the distance\n";
  cout << "  to a cylinder.\n";

  cout << "\n";
  cout << "  Radius R = " << r << "\n";
  cout << "  Center of bottom disk ="
       << "  " << setw(10) << p1[0]
       << "  " << setw(10) << p1[1]
       << "  " << setw(10) << p1[2] << "\n";
  cout << "  Center of top disk =   "
       << "  " << setw(10) << p2[0]
       << "  " << setw(10) << p2[1]
       << "  " << setw(10) << p2[2] << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    p = p_test + test * DIM_NUM;

    cout << "\n";
    cout << "  P ="
         << "  " << setw(10) << p[0]
         << "  " << setw(10) << p[1]
         << "  " << setw(10) << p[2] << "\n";

    dist = cylinder_point_dist_3d ( p1, p2, r, p );

    cout << "  Distance (computed) ="
         << "  " << setw(10) << dist << "\n";
    cout << "  Distance (exact)    ="
         << "  " << setw(10) << dist_test[test] << "\n";
  }

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

void test02015 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST02015 tests CYLINDER_POINT_DIST_SIGNED_3D.
//
//  Modified:
//
//    26 August 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3
# define TEST_NUM 6

  double dist;
  double dist_test[TEST_NUM] = {
    3.0, -0.5, 5.0, 8.0, -1.0, -0.25 };
  double *p;
  double p_test[DIM_NUM*TEST_NUM] = {
      4.0,    0.5,  0.0,
     -0.5,   -1.0,  0.0,
      4.0,    6.0,  0.0,
      0.75, -10.0,  0.0,
      0.0,    0.0,  0.0,
      0.25,   1.75, 0.0 };
  double p1[DIM_NUM] = { 0.0, -2.0, 0.0 };
  double p2[DIM_NUM] = { 0.0,  2.0, 0.0 };
  double r = 1.0;
  int test;

  cout << "\n";
  cout << "TEST02015\n";
  cout << "  CYLINDER_POINT_DIST_SIGNED_3D computes the signed\n";
  cout << "  distance to a cylinder.\n";

  cout << "\n";
  cout << "  Radius R = " << r << "\n";
  cout << "  Center of bottom disk ="
       << "  " << setw(10) << p1[0]
       << "  " << setw(10) << p1[1]
       << "  " << setw(10) << p1[2] << "\n";
  cout << "  Center of top disk =   "
       << "  " << setw(10) << p2[0]
       << "  " << setw(10) << p2[1]
       << "  " << setw(10) << p2[2] << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    p = p_test + test * DIM_NUM;

    cout << "\n";
    cout << "  P ="
         << "  " << setw(10) << p[0]
         << "  " << setw(10) << p[1]
         << "  " << setw(10) << p[2] << "\n";

    dist = cylinder_point_dist_signed_3d ( p1, p2, r, p );

    cout << "  Signed distance (computed) ="
         << "  " << setw(10) << dist << "\n";
    cout << "  Signed distance (exact)    ="
         << "  " << setw(10) << dist_test[test] << "\n";
  }

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

void test0202 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0202 tests CYLINDER_POINT_INSIDE_3D.
//
//  Modified:
//
//    23 August 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3
# define TEST_NUM 6

  bool inside;
  bool inside_test[TEST_NUM] = {
    false, true, false, false, true, true };
  double *p;
  double p_test[DIM_NUM*TEST_NUM] = {
      4.0,    0.5,  0.0,
     -0.5,   -1.0,  0.0,
      4.0,    6.0,  0.0,
      0.75, -10.0,  0.0,
      0.0,    0.0,  0.0,
      0.25,   1.75, 0.0 };
  double p1[DIM_NUM] = { 0.0, -2.0, 0.0 };
  double p2[DIM_NUM] = { 0.0,  2.0, 0.0 };
  double r = 1.0;
  int test;

  cout << "\n";
  cout << "TEST0202\n";
  cout << "  CYLINDER_POINT_INSIDE_3D determines if a point\n";
  cout << "  is inside a cylinder.\n";

  cout << "\n";
  cout << "  Radius R = " << r << "\n";
  cout << "  Center of bottom disk ="
       << "  " << setw(10) << p1[0]
       << "  " << setw(10) << p1[1]
       << "  " << setw(10) << p1[2] << "\n";
  cout << "  Center of top disk =   "
       << "  " << setw(10) << p2[0]
       << "  " << setw(10) << p2[1]
       << "  " << setw(10) << p2[2] << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    p = p_test + test * DIM_NUM;

    cout << "\n";
    cout << "  P ="
         << "  " << setw(10) << p[0]
         << "  " << setw(10) << p[1]
         << "  " << setw(10) << p[2] << "\n";

    inside = cylinder_point_inside_3d ( p1, p2, r, p );

    cout << "  INSIDE (computed) ="
         << "  " << setw(1) << inside << "\n";
    cout << "  INSIDE (exact)    ="
         << "  " << setw(1) << inside_test[test] << "\n";
  }

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

void test0203 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0203 tests CYLINDER_POINT_NEAR_3D.
//
//  Modified:
//
//    23 August 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3
# define TEST_NUM 6

  double *p;
  double p_test[DIM_NUM*TEST_NUM] = {
      4.0,    0.5,  0.0,
     -0.5,   -1.0,  0.0,
      4.0,    6.0,  0.0,
      0.75, -10.0,  0.0,
      0.0,    0.0,  0.0,
      0.25,   1.75, 0.0 };
  double p1[DIM_NUM] = { 0.0, -2.0, 0.0 };
  double p2[DIM_NUM] = { 0.0,  2.0, 0.0 };
  double *pn;
  double pn_test[DIM_NUM*TEST_NUM] = {
      1.0,   0.5,  0.0,
     -1.0,  -1.0,  0.0,
      1.0,   2.0,  0.0,
      0.75, -2.0,  0.0,
      1.0,   0.0,  0.0,
      0.25,  2.0,  0.0 };
  double r = 1.0;
  int test;

  cout << "\n";
  cout << "TEST0203\n";
  cout << "  CYLINDER_POINT_NEAR_3D computes the nearest point\n";
  cout << "  on a cylinder.\n";

  cout << "\n";
  cout << "  Radius R = " << r << "\n";
  cout << "  Center of bottom disk ="
       << "  " << setw(10) << p1[0]
       << "  " << setw(10) << p1[1]
       << "  " << setw(10) << p1[2] << "\n";
  cout << "  Center of top disk =   "
       << "  " << setw(10) << p2[0]
       << "  " << setw(10) << p2[1]
       << "  " << setw(10) << p2[2] << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    p = p_test + test * DIM_NUM;

    cout << "\n";
    cout << "  P ="
         << "  " << setw(10) << p[0]
         << "  " << setw(10) << p[1]
         << "  " << setw(10) << p[2] << "\n";

    pn = cylinder_point_near_3d ( p1, p2, r, p );

    cout << "  PN (computed) ="
         << "  " << setw(10) << pn[0]
         << "  " << setw(10) << pn[1]
         << "  " << setw(10) << pn[2] << "\n";
    cout << "  PN (exact)    ="
         << "  " << setw(10) << pn_test[0+test*DIM_NUM]
         << "  " << setw(10) << pn_test[1+test*DIM_NUM]
         << "  " << setw(10) << pn_test[2+test*DIM_NUM] << "\n";
    delete [] pn;
  }

  cout << "\n";
  cout << "  (Note that case 5 is ambiguous.  The set of nearest\n";
  cout << "  points forms a circle, any of which will do.)\n";

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

void test02035 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST02035 tests CYLINDER_SAMPLE_3D.
//
//  Modified:
//
//    27 August 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3
# define N 20

  double *p;
  double p1[DIM_NUM] = { 0.0, -2.0, 0.0 };
  double p2[DIM_NUM] = { 0.0,  2.0, 0.0 };
  double r = 1.0;
  int seed = 123456789;

  cout << "\n";
  cout << "TEST02035\n";
  cout << "  CYLINDER_SAMPLE_3D samples points in a cylinder.\n";

  cout << "\n";
  cout << "  Radius R = " << r << "\n";
  cout << "  Center of bottom disk ="
       << "  " << setw(10) << p1[0]
       << "  " << setw(10) << p1[1]
       << "  " << setw(10) << p1[2] << "\n";
  cout << "  Center of top disk =   "
       << "  " << setw(10) << p2[0]
       << "  " << setw(10) << p2[1]
       << "  " << setw(10) << p2[2] << "\n";

  p = cylinder_sample_3d ( p1, p2, r, N, &seed );

  r8mat_transpose_print ( DIM_NUM, N, p, "  Sample points:" );

  delete [] p;

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

void test0204 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0204 tests CYLINDER_VOLUME_3D.
//
//  Modified:
//
//    21 August 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  double p1[DIM_NUM] = { 1.0, 2.0, 3.0 };
  double p2[DIM_NUM] = { 5.0, 6.0, 5.0 };
  double pi = 3.141592653589793;
  double r = 5.0;
  double volume;

  cout << "\n";
  cout << "TEST0204\n";
  cout << "  CYLINDER_VOLUME_3D computes the volume of a cylinder.\n";

  cout << "\n";
  cout << "  Radius R = " << r << "\n";
  cout << "  Center of bottom disk ="
       << "  " << setw(10) << p1[0]
       << "  " << setw(10) << p1[1]
       << "  " << setw(10) << p1[2] << "\n";
  cout << "  Center of top disk =   "
       << "  " << setw(10) << p2[0]
       << "  " << setw(10) << p2[1]
       << "  " << setw(10) << p2[2] << "\n";

  volume = cylinder_volume_3d ( p1, p2, r );

  cout << "\n";
  cout << "  Volume (computed) = " << volume << "\n";
  cout << "  Volume (exact)    = " << pi * 150.0 << "\n";

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

void test0205 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0205 tests DEGREES_TO_RADIANS and RADIANS_TO_DEGREES.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
  double angle_deg;
  double angle_deg2;
  double angle_rad;
  int i;

  cout << "\n";
  cout << "TEST0205\n";
  cout << "  DEGREES_TO_RADIANS converts an angle from degrees\n";
  cout << "    to radians;\n";
  cout << "  RADIANS_TO_DEGREES converts an angle from radians\n";
  cout << "    to degrees;\n";
  cout << "\n";
  cout << "  Degrees     Radians     Degrees\n";
  cout << "\n";

  for ( i = -2; i <= 14; i++ )
  {
    angle_deg = ( double ) ( 30 * i );
    angle_rad = degrees_to_radians ( angle_deg );
    angle_deg2 = radians_to_degrees ( angle_rad );
    cout << "  " << setw(10) << angle_deg
         << "  " << setw(10) << angle_rad
         << "  " << setw(10) << angle_deg2 << "\n";
  }

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

void test021 ( void )

//****************************************************************************80
//
//  Purpose:
//
//    TEST021 tests DIRECTION_PERT_3D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  int i;
  int itest;
  int seed;
  double sigma;
  double vbase[DIM_NUM] = { 1.0, 0.0, 0.0 };
  double *vran;

  cout << "\n";
  cout << "TEST021\n";
  cout << "  DIRECTION_PERT_3D perturbs a direction vector.\n";
  cout << "\n";

  seed = get_seed ( );

  cout << "\n";
  cout << "  We use SEED = " << seed << "\n";

  cout << "\n";
  cout << "  Base vector:\n";
  cout << "  " << vbase[0]
       << "  " << vbase[1]
       << "  " << vbase[2] << "\n";

  for ( itest = 0; itest < 3; itest++ )
  {
    if ( itest == 0 )
    {
      sigma = 0.99;
    }
    else if ( itest == 1 )
    {
      sigma = 0.5;
    }
    else
    {
      sigma = 0.1;
    }

    cout << "\n";
    cout << "  Using sigma = " << sigma << "\n";
    cout << "\n";

    for ( i = 0; i < 20; i++ )
    {
      vran = direction_pert_3d ( sigma, vbase, &seed );
      cout << "  " << setw(10) << vran[0]
           << "  " << setw(10) << vran[1]
           << "  " << setw(10) << vran[2] << "\n";
      delete [] vran;
    }

  }

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

void test022 ( void )

//****************************************************************************80
//
//  Purpose:
//
//    TEST022 tests DIRECTION_UNIFORM_3D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
  int i;
  int seed;
  double *vran;

  cout << "\n";
  cout << "TEST022\n";
  cout << "  DIRECTION_UNIFORM_3D picks a random direction vector.\n";

  seed = get_seed ( );

  cout << "\n";
  cout << "  We use SEED = " << seed << "\n";

  cout << "\n";

  for ( i = 0; i < 10; i++ )
  {
    vran = direction_uniform_3d ( &seed );
    cout << "  " << setw(10) << vran[0]
         << "  " << setw(10) << vran[1]
         << "  " << setw(10) << vran[2] << "\n";
    delete [] vran;
  }

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

void test023 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST023 tests DIRECTION_UNIFORM_ND and GET_SEED.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 4

  int i;
  int j;
  int seed = 123456789;
  double *vran;

  cout << "\n";
  cout << "TEST023\n";
  cout << "  DIRECTION_UNIFORM_ND picks a random direction vector.\n";
  cout << "\n";

  for ( i = 1; i <= 10; i++ )
  {
    vran = direction_uniform_nd ( DIM_NUM, &seed );
    for ( j = 0; j < DIM_NUM; j++ )
    {
      cout << "  " << setw(8) << vran[j];
    }
    cout << "\n";
    delete [] vran;
  }

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

void test0232 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0232 tests DISK_POINT_DIST_3D.
//
//  Modified:
//
//    17 August 2005
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3
# define TEST_NUM 5

  double axis[DIM_NUM] = { 0.0, 1.0, 1.0 };
  double dist;
  double dist_test[TEST_NUM] = { 2.0, 0.0, 0.0, 8.0, 10.0 };
  double *p;
  double pc[DIM_NUM] = { 0.0, 1.4142135, 1.4142135 };
  double p_test[DIM_NUM*TEST_NUM] = {
     0.0,  0.0,        0.0,
     0.0,  0.70710677, 2.1213202,
     2.0,  1.4142135,  1.4142135,
    10.0,  1.4142135,  1.4142135,
    10.0,  5.6568542,  5.6568542 };
  double r = 2.0;
  int test;

  cout << "\n";
  cout << "TEST0232\n";
  cout << "  DISK_POINT_DIST_3D finds the distance from\n";
  cout << "    a disk to a point in 3D.\n";

  cout << "\n";
  cout << "  Disk radius = " << r << "\n";
  r8vec_print ( DIM_NUM, pc, "  Disk center: " );
  r8vec_print ( DIM_NUM, axis, "  Disk axis: " );

  for ( test = 0; test < TEST_NUM; test++ )
  {
    p = p_test + test * DIM_NUM;

    r8vec_print ( DIM_NUM, p, "  Point: " );

    dist = disk_point_dist_3d ( pc, r, axis, p );

    cout << "\n";
    cout << "  Distance = " << dist
         << "  Expected = " << dist_test[test] << "\n";
  }

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

void test0234 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0234 tests R8MAT_SOLVE_2D.
//
//  Modified:
//
//    16 November 2005
//
//  Author:
//
//    John Burkardt
//
{
  double *a;
  double *b;
  double det;
  int i;
  int n = 2;
  int seed;
  int test;
  int test_num = 5;
  double *x;
  double *x2;

  cout << "\n";
  cout << "TEST0234\n";
  cout << "  R8MAT_SOLVE_2D solves 2D linear systems.\n";

  seed = 123456789;

  for ( test = 1; test <= test_num; test++ )
  {
    a = r8mat_uniform_01 ( n, n, &seed );
    x = r8vec_uniform_01 ( n, &seed );
    b = r8mat_mv ( n, n, a, x );

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

    cout << "\n";
    cout << "  Solution / Computed:\n";
    cout << "\n";

    for ( i = 0; i < n; i++ )
    {
      cout << "  " << setw(14) << x[i]
           << "  " << setw(14) << x2[i] << "\n";
    }

    delete [] a;
    delete [] b;
    delete [] x;
    delete [] x2;
  }

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

void test0235 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0235 tests DMS_TO_RADIANS and RADIANS_TO_DMS.
//
//  Modified:
//
//    17 January 2007
//
//  Author:
//
//    John Burkardt
//
{
  int angle_deg;
  int angle_min;
  double angle_rad;
  double angle_rad2;
  int angle_sec;
  int i;
  double pi = 3.141592653589793;

  cout << "\n";
  cout << "TEST0235\n";
  cout << "  DMS_TO_RADIANS converts an angle from\n";
  cout << "    degrees/minutes/seconds to radians;\n";
  cout << "  RADIANS_TO_DEGREES converts an angle from radians\n";
  cout << "    to degrees/minutes/seconds;\n";
  cout << "\n";
  cout << "  Radians     DMS     Radians\n";
  wcout << "\n";

  for ( i = -2; i <= 15; i++ )
  {
    angle_rad = pi * ( double ) ( i ) / 7.0;

    radians_to_dms ( angle_rad, &angle_deg, &angle_min, &angle_sec );

    angle_rad2 = dms_to_radians ( angle_deg, angle_min, angle_sec );

    cout << "  " << setw(10) << angle_rad
         << "  " << setw(4)  << angle_deg
         << "  " << setw(3)  << angle_min
         << "  " << setw(3)  << angle_sec
         << "  " << setw(10) << angle_rad2 << "\n";
  }

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

void test0236 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0236 tests DODEC_SIZE_3D, DODEC_SHAPE_3D, SHAPE_PRINT_3D.
//
//  Modified:
//
//    22 July 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  int edge_num;
  int face_num;
  int *face_order;
  int face_order_max;
  int *face_point;
  int point_num;
  double *point_coord;

  cout << "\n";
  cout << "TEST0236\n";
  cout << "  For the dodecahedron,\n";
  cout << "  DODEC_SIZE_3D returns dimension information;\n";
  cout << "  DODEC_SHAPE_3D returns face and order information.\n";
  cout << "  SHAPE_PRINT_3D prints this information.\n";
//
//  Get the sizes.
//
  dodec_size_3d ( &point_num, &edge_num, &face_num, &face_order_max );

  cout << "\n";
  cout << "    Number of vertices: " << point_num << "\n";
  cout << "    Number of edges   : " << edge_num << "\n";
  cout << "    Number of faces   : " << face_num << "\n";
  cout << "    Maximum face order: " << face_order_max << "\n";
//
//  Make room for the data.
//
  face_order = new int[face_num];
  face_point = new int[face_order_max*face_num];
  point_coord = new double[DIM_NUM*point_num];
//
//  Get the data.
//
  dodec_shape_3d ( point_num, face_num, face_order_max, point_coord,
    face_order, face_point );
//
//  Print the data.
//
  shape_print_3d ( point_num, face_num, face_order_max,
    point_coord, face_order, face_point );

  delete [] face_order;
  delete [] face_point;
  delete [] point_coord;

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

void test0238 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0238 tests DUAL_SIZE_3D;
//
//  Modified:
//
//    22 July 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  int edge_num1;
  int edge_num2;
  int *edge_point1;
  int face_num1;
  int face_num2;
  int *face_order1;
  int face_order_max1;
  int face_order_max2;
  int *face_point1;
  int point_num1;
  int point_num2;
  double *point_coord1;

  cout << "\n";
  cout << "TEST0238\n";
  cout << "  DUAL_SIZE_3D finds the sizes of the dual of a\n";
  cout << "    polyhedron;\n";
//
//  Get the CUBE shape.
//
  cube_size_3d ( &point_num1, &edge_num1, &face_num1, &face_order_max1 );

  cout << "\n";
  cout << "  The cube:\n";
  cout << "    Number of vertices: " << point_num1 << "\n";
  cout << "    Number of edges   : " << edge_num1 << "\n";
  cout << "    Number of faces   : " << face_num1 << "\n";
  cout << "    Maximum face order: " << face_order_max1 << "\n";

  face_order1 = new int[face_num1];
  face_point1 = new int[face_order_max1*face_num1];
  point_coord1 = new double[DIM_NUM*point_num1];

  cube_shape_3d ( point_num1, face_num1, face_order_max1, point_coord1,
    face_order1, face_point1 );

  dual_size_3d ( point_num1, edge_num1, face_num1, face_order_max1, 
    point_coord1, face_order1, face_point1, &point_num2, &edge_num2,
    &face_num2, &face_order_max2 );

  cout << "\n";
  cout << "  The dual of the cube:\n";
  cout << "    Number of vertices: " << point_num2 << "\n";
  cout << "    Number of edges   : " << edge_num2 << "\n";
  cout << "    Number of faces   : " << face_num2 << "\n";
  cout << "    Maximum face order: " << face_order_max2 << "\n";

  delete [] face_order1;
  delete [] face_point1;
  delete [] point_coord1;
//
//  Get the DODECAHEDRON shape.
//
  dodec_size_3d ( &point_num1, &edge_num1, &face_num1, &face_order_max1 );

  cout << "\n";
  cout << "  The dodecahedron:\n";
  cout << "    Number of vertices: " << point_num1 << "\n";
  cout << "    Number of edges   : " << edge_num1 << "\n";
  cout << "    Number of faces   : " << face_num1 << "\n";
  cout << "    Maximum face order: " << face_order_max1 << "\n";

  face_order1 = new int[face_num1];
  face_point1 = new int[face_order_max1*face_num1];
  point_coord1 = new double[DIM_NUM*point_num1];

  dodec_shape_3d ( point_num1, face_num1, face_order_max1, point_coord1,
    face_order1, face_point1 );

  dual_size_3d ( point_num1, edge_num1, face_num1, face_order_max1, 
    point_coord1, face_order1, face_point1, &point_num2, &edge_num2,
    &face_num2, &face_order_max2 );

  cout << "\n";
  cout << "  The dual of the dodecahedron:\n";
  cout << "    Number of vertices: " << point_num2 << "\n";
  cout << "    Number of edges   : " << edge_num2 << "\n";
  cout << "    Number of faces   : " << face_num2 << "\n";
  cout << "    Maximum face order: " << face_order_max2 << "\n";

  delete [] face_order1;
  delete [] face_point1;
  delete [] point_coord1;
//
//  Get the ICOSAHEDRON shape.
//
  icos_size_3d ( &point_num1, &edge_num1, &face_num1, &face_order_max1 );

  cout << "\n";
  cout << "  The icosahedron:\n";
  cout << "    Number of vertices: " << point_num1 << "\n";
  cout << "    Number of edges   : " << edge_num1 << "\n";
  cout << "    Number of faces   : " << face_num1 << "\n";
  cout << "    Maximum face order: " << face_order_max1 << "\n";

  edge_point1 = new int[2*edge_num1];
  face_order1 = new int[face_num1];
  face_point1 = new int[face_order_max1*face_num1];
  point_coord1 = new double[DIM_NUM*point_num1];

  icos_shape_3d ( point_num1, edge_num1, face_num1, face_order_max1, 
    point_coord1, edge_point1, face_order1, face_point1 );

  dual_size_3d ( point_num1, edge_num1, face_num1, face_order_max1, 
    point_coord1, face_order1, face_point1, &point_num2, &edge_num2,
    &face_num2, &face_order_max2 );

  cout << "\n";
  cout << "  The dual of the icosahedron:\n";
  cout << "    Number of vertices: " << point_num2 << "\n";
  cout << "    Number of edges   : " << edge_num2 << "\n";
  cout << "    Number of faces   : " << face_num2 << "\n";
  cout << "    Maximum face order: " << face_order_max2 << "\n";

  delete [] edge_point1;
  delete [] face_order1;
  delete [] face_point1;
  delete [] point_coord1;
//
//  Get the octahedron shape.
//
  octahedron_size_3d ( &point_num1, &edge_num1, &face_num1, &face_order_max1 );

  cout << "\n";
  cout << "  The octahedron:\n";
  cout << "    Number of vertices: " << point_num1 << "\n";
  cout << "    Number of edges   : " << edge_num1 << "\n";
  cout << "    Number of faces   : " << face_num1 << "\n";
  cout << "    Maximum face order: " << face_order_max1 << "\n";

  face_order1 = new int[face_num1];
  face_point1 = new int[face_order_max1*face_num1];
  point_coord1 = new double[DIM_NUM*point_num1];

  octahedron_shape_3d ( point_num1, face_num1, face_order_max1,
    point_coord1, face_order1, face_point1 );

  dual_size_3d ( point_num1, edge_num1, face_num1, face_order_max1, 
    point_coord1, face_order1, face_point1, &point_num2, &edge_num2,
    &face_num2, &face_order_max2 );

  cout << "\n";
  cout << "  The dual of the octahedron:\n";
  cout << "    Number of vertices: " << point_num2 << "\n";
  cout << "    Number of edges   : " << edge_num2 << "\n";
  cout << "    Number of faces   : " << face_num2 << "\n";
  cout << "    Maximum face order: " << face_order_max2 << "\n";

  delete [] face_order1;
  delete [] face_point1;
  delete [] point_coord1;
//
//  Get the soccer ball shape.
//
  soccer_size_3d ( &point_num1, &edge_num1, &face_num1, &face_order_max1 );

  cout << "\n";
  cout << "  The soccer ball:\n";
  cout << "    Number of vertices: " << point_num1 << "\n";
  cout << "    Number of edges   : " << edge_num1 << "\n";
  cout << "    Number of faces   : " << face_num1 << "\n";
  cout << "    Maximum face order: " << face_order_max1 << "\n";

  face_order1 = new int[face_num1];
  face_point1 = new int[face_order_max1*face_num1];
  point_coord1 = new double[DIM_NUM*point_num1];

  soccer_shape_3d ( point_num1, face_num1, face_order_max1,
    point_coord1, face_order1, face_point1 );

  dual_size_3d ( point_num1, edge_num1, face_num1, face_order_max1, 
    point_coord1, face_order1, face_point1, &point_num2, &edge_num2,
    &face_num2, &face_order_max2 );

  cout << "\n";
  cout << "  The dual of the soccer ball:\n";
  cout << "    Number of vertices: " << point_num2 << "\n";
  cout << "    Number of edges   : " << edge_num2 << "\n";
  cout << "    Number of faces   : " << face_num2 << "\n";
  cout << "    Maximum face order: " << face_order_max2 << "\n";

  delete [] face_order1;
  delete [] face_point1;
  delete [] point_coord1;
//
//  Get the tetrahedron shape.
//
  tetrahedron_size_3d ( &point_num1, &edge_num1, &face_num1, &face_order_max1 );

  cout << "\n";
  cout << "  The tetrahedron:\n";
  cout << "    Number of vertices: " << point_num1 << "\n";
  cout << "    Number of edges   : " << edge_num1 << "\n";
  cout << "    Number of faces   : " << face_num1 << "\n";
  cout << "    Maximum face order: " << face_order_max1 << "\n";

  face_order1 = new int[face_num1];
  face_point1 = new int[face_order_max1*face_num1];
  point_coord1 = new double[DIM_NUM*point_num1];

  tetrahedron_shape_3d ( point_num1, face_num1, face_order_max1,
    point_coord1, face_order1, face_point1 );

  dual_size_3d ( point_num1, edge_num1, face_num1, face_order_max1, 
    point_coord1, face_order1, face_point1, &point_num2, &edge_num2,
    &face_num2, &face_order_max2 );

  cout << "\n";
  cout << "  The dual of the tetrahedron:\n";
  cout << "    Number of vertices: " << point_num2 << "\n";
  cout << "    Number of edges   : " << edge_num2 << "\n";
  cout << "    Number of faces   : " << face_num2 << "\n";
  cout << "    Maximum face order: " << face_order_max2 << "\n";

  delete [] face_order1;
  delete [] face_point1;
  delete [] point_coord1;

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

void test024 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST024 tests DUAL_SHAPE_3D;
//
//  Modified:
//
//    22 July 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  int edge_num1;
  int edge_num2;
  int face_num1;
  int face_num2;
  int *face_order1;
  int *face_order2;
  int face_order_max1;
  int face_order_max2;
  int *face_point1;
  int *face_point2;
  int point_num1;
  int point_num2;
  double *point_coord1;
  double *point_coord2;

  cout << "\n";
  cout << "TEST024\n";
  cout << "  DUAL_SHAPE_3D finds the dual of a polyhedron.\n";
//
//  Get the dodecahedron shape.
//
  cout << "\n";
  cout << "  The dodecahedron:\n";

  dodec_size_3d ( &point_num1, &edge_num1, &face_num1, &face_order_max1 );

  cout << "\n";
  cout << "    Number of vertices: " << point_num1 << "\n";
  cout << "    Number of edges   : " << edge_num1 << "\n";
  cout << "    Number of faces   : " << face_num1 << "\n";
  cout << "    Maximum face order: " << face_order_max1 << "\n";

  face_order1 = new int[face_num1];
  face_point1 = new int [face_order_max1*face_num1];
  point_coord1 = new double[DIM_NUM*point_num1];

  dodec_shape_3d ( point_num1, face_num1, face_order_max1, point_coord1,
    face_order1, face_point1 );

  shape_print_3d ( point_num1, face_num1, face_order_max1,
    point_coord1, face_order1, face_point1 );
//
//  Get the dual.
//
  cout << "\n";
  cout << "  The dual of the dodecahedron:\n";

  dual_size_3d ( point_num1, edge_num1, face_num1, face_order_max1, 
    point_coord1, face_order1, face_point1, &point_num2, &edge_num2,
    &face_num2, &face_order_max2 );

  cout << "\n";
  cout << "    Number of vertices: " << point_num2 << "\n";
  cout << "    Number of edges   : " << edge_num2 << "\n";
  cout << "    Number of faces   : " << face_num2 << "\n";
  cout << "    Maximum face order: " << face_order_max2 << "\n";

  face_order2 = new int[face_num2];
  face_point2 = new int[face_order_max2*face_num2];
  point_coord2 = new double[DIM_NUM*point_num2];

  dual_shape_3d ( point_num1, face_num1, face_order_max1, point_coord1,
    face_order1, face_point1, point_num2, face_num2, face_order_max2,
    point_coord2, face_order2, face_point2 );

  shape_print_3d ( point_num2, face_num2, face_order_max2,
    point_coord2, face_order2, face_point2 );

  delete [] face_order1;
  delete [] face_order2;
  delete [] face_point1;
  delete [] face_point2;
  delete [] point_coord1;
  delete [] point_coord2;

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

void test0243 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0243 tests R8VEC_ANY_NORMAL.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 10
# define TEST_NUM 5

  int seed;
  int test;
  double *v1;
  double v1_length;
  double v1v2_dot;
  double *v2;
  double v2_length;

  cout << "\n";
  cout << "TEST0243\n";
  cout << "  R8VEC_ANY_NORMAL computes a vector V2 that is normal\n";
  cout << "  to a given vector V1.\n";
  cout << "\n";
  cout << "    Test    ||V1||      ||V2||        V1.V2\n";
  cout << "\n";

  seed = 123456789;

  for ( test = 0; test < TEST_NUM; test++ )
  {
    v1 = r8vec_uniform_01 ( DIM_NUM, &seed );
    v1_length = r8vec_length ( DIM_NUM, v1 );
    v2 = r8vec_any_normal ( DIM_NUM, v1 );
    v2_length = r8vec_length ( DIM_NUM, v2 );
    v1v2_dot = r8vec_dot ( DIM_NUM, v1, v2 );
    cout << "  " << setw(6)  << test
         << "  " << setw(10) << v1_length
         << "  " << setw(10) << v2_length
         << "  " << setw(10) << v1v2_dot << "\n";
  }

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

void test0245 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0245 tests R8VEC_NORMAL_01.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define N_MAX 1000

  int i;
  int j;
  int n;
  int seed = 123456789;
  double *x;
  double x_max;
  double x_mean;
  double x_min;
  double x_var;

  cout << "\n";
  cout << "TEST0245\n";
  cout << "  R8VEC_NORMAL_01 computes a vector of normally\n";
  cout << "  distributed random numbers.\n";
  cout << "  Using initial random number seed = " << seed << "\n";
//
//  Test 1:
//  Simply call 5 times for 1 value, and print.
//
  cout << "\n";
  cout << "  Test #1: Call 5 times, 1 value each time.\n";
  cout << "\n";

  for ( i = 1; i <= 5; i++ )
  {
    n = 1;
    x = r8vec_normal_01 ( n, &seed );
    for ( j = 0; j < n; j++ )
    {
      cout << "  " << setw(6)  << i
           << "  " << setw(10) << x[j] << "\n";
    }
    delete [] x;
  }
//
//  Test 2:
//  Restore the random number seed, and repeat.
//
  cout << "\n";
  cout << "  Test #2: Restore the random number seed.\n";
  cout << "  call 5 times, 1 value each time.\n";
  cout << "  The results should be identical.\n";
  cout << "\n";

  n = -1;
  x = r8vec_normal_01 ( n, &seed );
  delete [] x;

  seed = 123456789;

  for ( i = 1; i <= 5; i++ )
  {
    n = 1;
    x = r8vec_normal_01 ( n, &seed );
    for ( j = 0; j < n; j++ )
    {
      cout << "  " << setw(6)  << i
           << "  " << setw(10) << x[j] << "\n";
    }
    delete [] x;
  }
//
//  Test 3:
//  Restore the random number seed, compute all 5 values at once.
//
  cout << "\n";
  cout << "  Test #3: Restore the random number seed.\n";
  cout << "  call 1 time for 5 values.\n";
  cout << "  The results should be identical.\n";
  cout << "\n";

  n = -1;
  x = r8vec_normal_01 ( n, &seed );
  delete [] x;

  seed = 123456789;

  n = 5;
  x = r8vec_normal_01 ( n, &seed );

  i = 0;
  for ( j = 0; j < n; j++ )
  {
    i = i + 1;
    cout << "  " << setw(6)  << i
         << "  " << setw(10) << x[j] << "\n";
  }
  delete [] x;
//
//  Test 4:
//  Restore the random number seed, compute all 5 values at once.
//
  cout << "\n";
  cout << "  Test #4: Restore the random number seed.\n";
  cout << "  call for 2, 1, and 2 values.\n";
  cout << "  The results should be identical.\n";
  cout << "\n";

  n = -1;
  x = r8vec_normal_01 ( n, &seed );
  delete [] x;

  seed = 123456789;

  n = 2;
  x = r8vec_normal_01 ( n, &seed );
  i = 0;

  for ( j = 0; j < n; j++ )
  {
    i = i + 1;
    cout << "  " << setw(6)  << i
         << "  " << setw(10) << x[j] << "\n";
  }
  delete [] x;

  n = 1;
  x = r8vec_normal_01 ( n, &seed );

  for ( j = 0; j < n; j++ )
  {
    i = i + 1;
    cout << "  " << setw(6)  << i
         << "  " << setw(10) << x[j] << "\n";
  }
  delete [] x;

  n = 2;
  x = r8vec_normal_01 ( n, &seed );

  for ( j = 0; j < n; j++ )
  {
    i = i + 1;
    cout << "  " << setw(6)  << i
         << "  " << setw(10) << x[j] << "\n";
  }
  delete [] x;
//
//  Test 5:
//  Determine the minimum, maximum, mean and variance.
//
  n = N_MAX;
  x = r8vec_normal_01 ( n, &seed );

  x_min = r8vec_min ( n, x );
  x_max = r8vec_max ( n, x );
  x_mean = r8vec_mean ( n, x );
  x_var = r8vec_variance ( n, x );

  cout << "\n";
  cout << "  Test #5:\n";
  cout << "  Number of samples was " << n       << "\n";
  cout << "  Minimum value was     " << x_min   << "\n";
  cout << "  Maximum value was     " << x_max   << "\n";
  cout << "  Average value was     " << x_mean  << "\n";
  cout << "  Variance was          " << x_var   << "\n";
  cout << "  Expected average      " << 0.0 << "\n";
  cout << "  Expected variance     " << 1.0 << "\n";

  delete [] x;

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

void test025 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST025 tests ELLIPSE_POINT_DIST_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  double dist;
  int i;
  int n = 10;
  double p[DIM_NUM];
  double r1 = 3.0;
  double r2 = 2.0;

  cout << "\n";
  cout << "TEST025:\n";
  cout << "  ELLIPSE_POINT_DIST_2D is given a point P, and\n";
  cout << "  finds the distance to an ellipse in 2D.\n";
  cout << "\n";
  cout << "  The ellipse is (X/R1)^2 + (Y/R2)^2 = 1\n";
  cout << "\n";
  cout << "  R1 = " << r1 << "\n";
  cout << "  R2 = " << r2 << "\n";
  cout << "\n";
  cout << "           P            DIST\n";
  cout << "\n";

  for ( i = -3; i <= n + 3; i++ )
  {
    p[0] = ( double ( n - i ) * 0.0
           + double (     i ) * 4.0 )
           / double ( n     );

    p[1] = ( double ( n - i ) * 3.0
           + double (     i ) * 0.0 )
           / double ( n     );

    dist = ellipse_point_dist_2d ( r1, r2, p );

    cout << "  " << setw(8) << p[0]
         << "  " << setw(8) << p[1]
         << "  " << setw(8) << dist << "\n";
  }

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

void test0255 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST0255 tests ELLIPSE_POINT_NEAR_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2

  int i;
  int n = 10;
  double p[DIM_NUM];
  double *pn;
  double r1 = 3.0;
  double r2 = 2.0;

  cout << "\n";
  cout << "TEST0255:\n";
  cout << "  ELLIPSE_POINT_NEAR_2D is given a point P, and\n";
  cout << "  finds the nearest point on an ellipse in 2D.\n";
  cout << "\n";
  cout << "  The ellipse is (X/R1)^2 + (Y/R2)^2 = 1\n";
  cout << "\n";
  cout << "  R1 = " << r1 << "\n";
  cout << "  R2 = " << r2 << "\n";
  cout << "\n";
  cout << "           P            PN\n";
  cout << "\n";

  for ( i = -3; i <= n + 3; i++ )
  {
    p[0] = ( double ( n - i ) * 0.0
           + double (     i ) * 4.0 )
           / double ( n     );

    p[1] = ( double ( n - i ) * 3.0
           + double (     i ) * 0.0 )
           / double ( n     );

    pn = ellipse_point_near_2d ( r1, r2, p );

    cout << "  " << setw(8) << p[0]
         << "  " << setw(8) << p[1]
         << "  " << setw(8) << pn[0]
         << "  " << setw(8) << pn[1] << "\n";

    delete [] pn;
  }

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

void test026 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST026 tests ELLIPSE_POINTS_2D, ELLIPSE_AREA_2D, POLYGON_AREA_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define N_MAX 24

  double area;
  int n;
  double pc[DIM_NUM] = { 5.0, -2.0 };
  double pi = 3.141592653589793;
  double psi = 3.141592653589793 / 6.0;
  double r1 = 3.0;
  double r2 = 1.0;
  double result;
  double *v;

  cout << "\n";
  cout << "TEST026\n";
  cout << "  ELLIPSE_POINTS_2D returns points on an ellipse;\n";
  cout << "  ELLIPSE_AREA_2D returns the area of an ellipse;\n";
  cout << "  POLYGON_AREA_2D finds the area of a polygon.\n";

  r8vec_print ( DIM_NUM, pc, "  Ellipse center:" );

  cout << "\n";
  cout << "  radii R1 = " << r1 << " R2 = " << r2 << "\n";
  cout << "  and angle PSI = " << psi << "\n";

  area = ellipse_area_2d ( r1, r2 );

  cout << "  and area = " << area << "\n";

  n = 16;
  v = new double[DIM_NUM*n];

  ellipse_points_2d ( pc, r1, r2, psi, n, v );

  r8mat_transpose_print ( DIM_NUM, n, v, "  Sample points:" );

  delete [] v;

  cout << "\n";
  cout << "  For any N, the sampled points define a polygon\n";
  cout << "  whose area approximates the ellipse area.\n";
  cout << "\n";
  cout << "       N         Area\n";
  cout << "\n";

  for ( n = 3; n <= N_MAX; n++ )
  {
    v = new double[DIM_NUM*n];
    ellipse_points_2d ( pc, r1, r2, psi, n, v );
    result = polygon_area_2d ( n, v );
    cout << "  " << setw(6) << n
         << "  " << setw(12) << result << "\n";
    delete [] v;
  }

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

void test027 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST027 tests ELLIPSE_POINTS_ARC_2D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define N 13

  double p[DIM_NUM*N];
  double pc[DIM_NUM] = { 5.0, -2.0 };
  double pi = 3.141592653589793;
  double psi = 3.141592653589793 / 6.0;
  double r1 = 3.0;
  double r2 = 1.0;
  double theta1 = 3.141592653589793 / 2.0;
  double theta2 = 2.0 * 3.141592653589793;

  cout << "\n";
  cout << "TEST027\n";
  cout << "  ELLIPSE_POINTS_ARC_2D returns points on an\n";
  cout << "    elliptical arc.\n";
  cout << "\n";
  cout << "  The ellipse has center " << pc[0] << "  " << pc[1] << "\n";
  cout << "  radii R1 = " << r1 << " R2 = " << r2 << "\n";
  cout << "  and angle PSI = " << psi << "\n";
  cout << "\n";
  cout << "  The arc extends from THETA1 = " << theta1 << "\n";
  cout << "  to THETA2 = " << theta2 << "\n";

  ellipse_points_arc_2d ( pc, r1, r2, psi, theta1, theta2, N, p );

  r8mat_transpose_print ( DIM_NUM, N, p, "  Sample points:" );

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

void test028 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST028 tests HALFPLANE_CONTAINS_POINT_2D
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 2
# define TEST_NUM 4

  bool expected[TEST_NUM] = { true, false, true, false };
  double *p;
  double p_test[DIM_NUM*TEST_NUM] = {
    1.0,   1.0,
    1.0,  -1.0,
   -1.0,   1.0,
    2.0, 200.0 };
  double *p1;
  double p1_test[DIM_NUM*TEST_NUM] = {
    0.0,   0.0,
    0.0,   0.0,
   -5.0,  -5.0,
    3.0, 150.0 };
  double *p2;
  double p2_test[DIM_NUM*TEST_NUM] = {
    2.0,   0.0,
    2.0,   0.0,
   10.0,  10.0,
    1.0,  50.0 };
  bool temp;
  int test;

  cout << "\n";
  cout << "TEST028\n";
  cout << "  HALFPLANE_CONTAINS_POINT_2D determines whether a\n";
  cout << "  halfplane bounded by PA:PB contains the\n";
  cout << "  point P.\n";
  cout << "\n";

  for ( test = 0; test < TEST_NUM; test++ )
  {
    p1 = p1_test + test*DIM_NUM;
    p2 = p2_test + test*DIM_NUM;
    p = p_test + test*DIM_NUM;

    temp = halfplane_contains_point_2d ( p1, p2, p );

    cout << "\n";
    cout << "  P1 = " << setw(12) << p1[0]
         << "  "      << setw(12) << p1[1] << "\n";
    cout << "  P2 = " << setw(12) << p2[0]
         << "  "      << setw(12) << p2[1] << "\n";
    cout << "  P = "  << setw(12) << p[0]
         << "  "      << setw(12) << p[1] << "\n";
    cout << "  Contains? = " << temp
         << "  Correct = " << expected[test] << "\n";
  }

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

void test029 ( void )

//****************************************************************************80
//
//  Purpose:
//
//    TEST029 tests HALFSPACE_IMP_TRIANGLE_INT_3D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  double a;
  double b;
  double c;
  double d;
  int int_num;
  int j;
  int test;
  int test_num = 6;
  double p[DIM_NUM*4];
  double t[DIM_NUM*3];

  cout << "\n";
  cout << "TEST029\n";
  cout << "  HALFSPACE_IMP_TRIANGLE_INT_3D finds\n";
  cout << "    intersection points of an implicit\n";
  cout << "    halfspace and a triangle.\n";

  a =   1.0;
  b = - 2.0;
  c = - 3.0;
  d =   6.0;

  cout << "\n";
  cout << "  The implicitly defined bounding plane\n";
  cout << "  has the form: A*X + B*Y + C*Z + D = 0.\n";
  cout << "  A,B,C,D = " << a << "  " << b << "  " << c << "  " << d << "\n";

  for ( test = 0; test < test_num; test++ )
    {

    if ( test == 0 )
    {
      t[0+0*3] =  0.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] = -1.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] = -2.0;
    }
    else if ( test == 1 )
    {
      t[0+0*3] = -6.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] = -1.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] = -2.0;
    }
    else if ( test == 2 )
    {
      t[0+0*3] =  0.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] =  3.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] =  2.0;
    }
    else if ( test == 3 )
    {
      t[0+0*3] = -6.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] =  4.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] =  3.0;
    }
    else if ( test == 4 )
    {
      t[0+0*3] = -8.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] = -1.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] = -2.0;
    }
    else if ( test == 5 )
    {
      t[0+0*3] =  0.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] =  4.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] =  4.0;
    }

    cout << "\n";
    cout << "  Case " << test << "\n";

    r8mat_transpose_print ( DIM_NUM, 3, t, "  Triangle vertices:" );

    int_num = halfspace_imp_triangle_int_3d ( a, b, c, d, t, p );

    cout << "\n";
    cout << "  Number of intersection points is " << int_num << "\n";
    cout << "\n";

    if ( 0 < int_num )
    {
      for ( j = 0; j < int_num; j++ )
      {
        cout << "  " << setw(4)  << j
             << "  " << setw(10) << p[0+j*DIM_NUM]
             << "  " << setw(10) << p[1+j*DIM_NUM]
             << "  " << setw(10) << p[2+j*DIM_NUM] << "\n";
      }
    }
  }

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

void test030 ( void )

//****************************************************************************80
//
//  Purpose:
//
//    TEST030 tests HALFSPACE_NORM_TRIANGLE_INT_3D.
//
//  Modified:
//
//    14 February 2007
//
//  Author:
//
//    John Burkardt
//
{
# define DIM_NUM 3

  int int_num;
  int j;
  double p[DIM_NUM*4];
  double pn[DIM_NUM] = { 2.0, -4.0, -6.0 };
  double pp[DIM_NUM] = { -6.0, 0.0, 0.0 };
  double t[DIM_NUM*3];
  int test;
  int test_num = 6;

  cout << "\n";
  cout << "TEST030\n";
  cout << "  HALFSPACE_NORM_TRIANGLE_INT_3D finds\n";
  cout << "    intersection points of a normal form\n";
  cout << "    halfspace and a triangle.\n";

  r8vec_print ( DIM_NUM, pp, "  A point on the plane:" );
  r8vec_print ( DIM_NUM, pn, "  The normal vector:" );

  for ( test = 0; test < test_num; test++ )
  {
    if ( test == 0 )
    {
      t[0+0*3] =  0.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] = -1.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] = -2.0;
    }
    else if ( test == 1 )
    {
      t[0+0*3] = -6.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] = -1.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] = -2.0;
    }
    else if ( test == 2 )
    {
      t[0+0*3] =  0.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] =  3.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] =  2.0;
    }
    else if ( test == 3 )
    {
      t[0+0*3] = -6.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] =  4.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] =  3.0;
    }
    else if ( test == 4 )
    {
      t[0+0*3] = -8.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] = -1.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] = -2.0;
    }
    else if ( test == 5 )
    {
      t[0+0*3] =  0.0;
      t[1+0*3] =  0.0;
      t[2+0*3] =  0.0;

      t[0+1*3] =  0.0;
      t[1+1*3] =  4.0;
      t[2+1*3] =  0.0;

      t[0+2*3] =  0.0;
      t[1+2*3] =  0.0;
      t[2+2*3] =  4.0;
    }


    cout << "\n";
    cout << "  Case " << test << "\n";
    r8mat_transpose_print ( DIM_NUM, 3, t, "  Triangle vertices:" );

    int_num = halfspace_norm_triangle_int_3d ( pp, pn, t, p );

    cout << "\n";
    cout << "  Number of intersection points is " << int_num << "\n";
    cout << "\n";

    if ( 0 < int_num )
    {
      for ( j = 0; j < int_num; j++ )
      {
        cout << "  " << setw(4)  << j
             << "  " << setw(10) << p[0+j*3]
             << "  " << setw(10) << p[1+j*3]
             << "  " << setw(10) << p[2+j*3] << "\n";
      }
    }

  }

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

void test031 ( void )

//****************************************************************************80********80
//
//  Purpose:
//
//    TEST031 tests HAVERSINE.
//
{
  double d;
  d