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

using namespace std;

# include "sphere_stereograph.hpp"

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

int i4_max ( int i1, int i2 )

//****************************************************************************80
//
//  Purpose:
//
//    I4_MAX returns the maximum of two I4's.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 October 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I1, I2, are two integers to be compared.
//
//    Output, int I4_MAX, the larger of I1 and I2.
//
{
  int value;

  if ( i2 < i1 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}
//****************************************************************************80

int i4_min ( int i1, int i2 )

//****************************************************************************80
//
//  Purpose:
//
//    I4_MIN returns the minimum of two I4's.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 October 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I1, I2, two integers to be compared.
//
//    Output, int I4_MIN, the smaller of I1 and I2.
//
{
  int value;

  if ( i1 < i2 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}
//****************************************************************************80

void plane_normal_basis_3d ( double pp[3], double pn[3], double pq[3],
  double pr[3] )

//****************************************************************************80
//
//  Purpose:
//
//    PLANE_NORMAL_BASIS_3D finds two perpendicular vectors in a plane in 3D.
//
//  Discussion:
//
//    The normal form of a plane in 3D is:
//
//      PP is a point on the plane,
//      N is a normal vector to the plane.
//
//    The two vectors to be computed, PQ and PR, can be regarded as
//    the basis of a Cartesian coordinate system for points in the plane.
//    Any point in the plane can be described in terms of the "origin"
//    point PP plus a weighted sum of the two vectors PQ and PR:
//
//      P = PP + a * PQ + b * PR.
//
//    The vectors PQ and PR have unit length, and are perpendicular to N
//    and to each other.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    27 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double PP[3], a point on the plane.
//
//    Input, double PN[3], a normal vector to the plane.  The
//    vector must not have zero length, but it is not necessary for PN
//    to have unit length.
//
//    Output, double PQ[3], a vector of unit length, perpendicular
//    to the vector PN and the vector PR.
//
//    Output, double PR[3], a vector of unit length, perpendicular
//    to the vector PN and the vector PQ.
//
{
# define DIM_NUM 3

  int i;
  double normal_norm;
  double pr_norm;
  double *temp;
//
//  Compute the length of NORMAL.
//
  normal_norm = r8vec_norm ( DIM_NUM, pn );

  if ( normal_norm == 0.0 )
  {
    cerr << "\n";
    cerr << "PLANE_NORMAL_BASIS_3D - Fatal error!\n";
    cerr << "  The normal vector is 0.\n";
    exit ( 1 );
  }
//
//  Find a vector PQ that is normal to PN and has unit length.
//
  temp = r8vec_any_normal ( DIM_NUM, pn );
  r8vec_copy ( DIM_NUM, temp, pq );
  delete [] temp;
//
//  Now just take the cross product PN x PQ to get the PR vector.
//
  temp = r8vec_cross_product_3d ( pn, pq );

  pr_norm = r8vec_norm ( DIM_NUM, temp );

  for ( i = 0; i < DIM_NUM; i++ )
  {
    pr[i] = temp[i] / pr_norm;
  }
  delete [] temp;

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

void r8mat_transpose_print ( int m, int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 September 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Input, double A[M*N], an M by N matrix to be printed.
//
//    Input, string TITLE, a title.
//
{
  r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title );

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

void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo,
  int ihi, int jhi, string title )

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 August 2010
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Input, double A[M*N], an M by N matrix to be printed.
//
//    Input, int ILO, JLO, the first row and column to print.
//
//    Input, int IHI, JHI, the last row and column to print.
//
//    Input, string TITLE, a title.
//
{
# define INCX 5

  int i;
  int i2;
  int i2hi;
  int i2lo;
  int inc;
  int j;
  int j2hi;
  int j2lo;

  cout << "\n";
  cout << title << "\n";

  if ( m <= 0 || n <= 0 )
  {
    cout << "\n";
    cout << "  (None)\n";
    return;
  }

  for ( i2lo = i4_max ( ilo, 1 ); i2lo <= i4_min ( ihi, m ); i2lo = i2lo + INCX )
  {
    i2hi = i2lo + INCX - 1;
    i2hi = i4_min ( i2hi, m );
    i2hi = i4_min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    cout << "\n";
    cout << "  Row: ";
    for ( i = i2lo; i <= i2hi; i++ )
    {
      cout << setw(7) << i - 1 << "       ";
    }
    cout << "\n";
    cout << "  Col\n";
    cout << "\n";

    j2lo = i4_max ( jlo, 1 );
    j2hi = i4_min ( jhi, n );

    for ( j = j2lo; j <= j2hi; j++ )
    {
      cout << setw(5) << j - 1 << ":";
      for ( i2 = 1; i2 <= inc; i2++ )
      {
        i = i2lo - 1 + i2;
        cout << setw(14) << a[(i-1)+(j-1)*m];
      }
      cout << "\n";
    }
  }

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

double *r8mat_uniform_01_new ( int m, int n, int *seed )

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_UNIFORM_01_NEW returns a unit pseudorandom R8MAT.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8's,  stored as a vector
//    in column-major order.
//
//    This routine implements the recursion
//
//      seed = 16807 * seed mod ( 2**31 - 1 )
//      unif = seed / ( 2**31 - 1 )
//
//    The integer arithmetic never requires more than 32 bits,
//    including a sign bit.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    03 October 2005
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Paul Bratley, Bennett Fox, Linus Schrage,
//    A Guide to Simulation,
//    Springer Verlag, pages 201-202, 1983.
//
//    Bennett Fox,
//    Algorithm 647:
//    Implementation and Relative Efficiency of Quasirandom
//    Sequence Generators,
//    ACM Transactions on Mathematical Software,
//    Volume 12, Number 4, pages 362-376, 1986.
//
//    Philip Lewis, Allen Goodman, James Miller,
//    A Pseudo-Random Number Generator for the System/360,
//    IBM Systems Journal,
//    Volume 8, pages 136-143, 1969.
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Input/output, int *SEED, the "seed" value.  Normally, this
//    value should not be 0, otherwise the output value of SEED
//    will still be 0, and R8_UNIFORM will be 0.  On output, SEED has
//    been updated.
//
//    Output, double R8MAT_UNIFORM_01_NEW[M*N], a matrix of pseudorandom values.
//
{
  int i;
  int j;
  int k;
  double *r;

  r = new double[m*n];

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      k = *seed / 127773;

      *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;

      if ( *seed < 0 )
      {
        *seed = *seed + 2147483647;
      }
//
//  Although SEED can be represented exactly as a 32 bit integer,
//  it generally cannot be represented exactly as a 32 bit real number//
//
      r[i+j*m] = ( double ) ( *seed ) * 4.656612875E-10;
    }
  }

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

double *r8vec_any_normal ( int dim_num, double v1[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_ANY_NORMAL returns some normal vector to V1.
//
//  Discussion:
//
//    If DIM_NUM < 2, then no normal vector can be returned.
//
//    If V1 is the zero vector, then any unit vector will do.
//
//    No doubt, there are better, more robust algorithms.  But I will take
//    just about ANY reasonable unit vector that is normal to V1.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, double V1[DIM_NUM], the vector.
//
//    Output, double R8VEC_ANY_NORMAL[DIM_NUM], a vector that is
//    normal to V2, and has unit Euclidean length.
//
{
  int i;
  int j;
  int k;
  double *v2;
  double vj;
  double vk;

  if ( dim_num < 2 )
  {
    cerr << "\n";
    cerr << "R8VEC_ANY_NORMAL - Fatal error!\n";
    cerr << "  Called with DIM_NUM < 2.\n";
    exit ( 1 );
  }

  v2 = new double[dim_num];

  if ( r8vec_norm ( dim_num, v1 ) == 0.0 )
  {
    r8vec_zero ( dim_num, v2 );
    v2[0] = 1.0;
    return v2;
  }
//
//  Seek the largest entry in V1, VJ = V1(J), and the
//  second largest, VK = V1(K).
//
//  Since V1 does not have zero norm, we are guaranteed that
//  VJ, at least, is not zero.
//
  j = -1;
  vj = 0.0;

  k = -1;
  vk = 0.0;

  for ( i = 0; i < dim_num; i++ )
  {
    if ( fabs ( vk ) < fabs ( v1[i] ) || k == -1 )
    {
      if ( fabs ( vj ) < fabs ( v1[i] ) || j == -1 )
      {
        k = j;
        vk = vj;
        j = i;
        vj = v1[i];
      }
      else
      {
        k = i;
        vk = v1[i];
      }
    }
  }
//
//  Setting V2 to zero, except that V2(J) = -VK, and V2(K) = VJ,
//  will just about do the trick.
//
  r8vec_zero ( dim_num, v2 );

  v2[j] = -vk / sqrt ( vk * vk + vj * vj );
  v2[k] =  vj / sqrt ( vk * vk + vj * vj );

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

void r8vec_copy ( int n, double a1[], double a2[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_COPY copies an R8VEC.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    03 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vectors.
//
//    Input, double A1[N], the vector to be copied.
//
//    Input, double A2[N], the copy of A1.
//
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    a2[i] = a1[i];
  }
  return;
}
//****************************************************************************80

double *r8vec_cross_product_3d ( double v1[3], double v2[3] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_CROSS_PRODUCT_3D computes the cross product of two R8VEC's in 3D.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    07 August 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double V1[3], V2[3], the coordinates of the vectors.
//
//    Output, double R8VEC_CROSS_PRODUCT_3D[3], the cross product vector.
//
{
  double *v3;

  v3 = new double[3];

  v3[0] = v1[1] * v2[2] - v1[2] * v2[1];
  v3[1] = v1[2] * v2[0] - v1[0] * v2[2];
  v3[2] = v1[0] * v2[1] - v1[1] * v2[0];

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

double r8vec_norm ( int n, double a[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_NORM returns the L2 norm of an R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//    The vector L2 norm is defined as:
//
//      R8VEC_NORM = sqrt ( sum ( 1 <= I <= N ) A(I)^2 ).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 March 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in A.
//
//    Input, double A[N], the vector whose L2 norm is desired.
//
//    Output, double R8VEC_NORM, the L2 norm of A.
//
{
  int i;
  double v;

  v = 0.0;

  for ( i = 0; i < n; i++ )
  {
    v = v + a[i] * a[i];
  }
  v = sqrt ( v );

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

double r8vec_norm_affine ( int n, double v0[], double v1[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_NORM_AFFINE returns the affine L2 norm of an R8VEC.
//
//  Discussion:
//
//    The affine vector L2 norm is defined as:
//
//      R8VEC_NORM_AFFINE(V0,V1)
//        = sqrt ( sum ( 1 <= I <= N ) ( V1(I) - V0(I) )^2
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    27 October 2010
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the dimension of the vectors.
//
//    Input, double V0[N], the base vector.
//
//    Input, double V1[N], the vector whose affine L2 norm is desired.
//
//    Output, double R8VEC_NORM_AFFINE, the affine L2 norm of V1.
//
{
  int i;
  double value;

  value = 0.0;

  for ( i = 0; i < n; i++ )
  {
    value = value + ( v1[i] - v0[i] ) * ( v1[i] - v0[i] );
  }
  value = sqrt ( value );

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

void r8vec_normal_01 ( int n, int *seed, double x[] )

//*****************************************************************************
//
//  Purpose:
//
//    R8VEC_NORMAL_01 samples the standard normal probability distribution.
//
//  Discussion:
//
//    The standard normal probability distribution function (PDF) has
//    mean 0 and standard deviation 1.
//
//    This routine can generate a vector of values on one call.  It
//    has the feature that it should provide the same results
//    in the same order no matter how we break up the task.
//
//    Before calling this routine, the user may call RANDOM_SEED
//    in order to set the seed of the random number generator.
//
//    The Box-Muller method is used, which is efficient, but
//    generates an even number of values each time.  On any call
//    to this routine, an even number of new values are generated.
//    Depending on the situation, one value may be left over.
//    In that case, it is saved for the next call.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    18 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of values desired.  If N is negative,
//    then the code will flush its internal memory; in particular,
//    if there is a saved value to be used on the next call, it is
//    instead discarded.  This is useful if the user has reset the
//    random number seed, for instance.
//
//    Input/output, int *SEED, a seed for the random number generator.
//
//    Output, double X[N], a sample of the standard normal PDF.
//
{
  int i;
  int m;
  double pi = 3.141592653589793;
  double *r;
  static int made = 0;
  static int saved = 0;
  int xhi;
  int xlo;
  static double y = 0.0;
//
//  I'd like to allow the user to reset the random number seed.
//  But this won't work properly if we have a saved value Y.
//  I'm making a crock option that allows the user to signal
//  explicitly that any internal memory should be flushed,
//  by passing in a negative value for N.
//
  if ( n < 0 )
  {
    made = 0;
    saved = 0;
    y = 0.0;
    return;
  }
  else if ( n == 0 )
  {
    return;
  }
//
//  Record the range of X we need to fill in.
//
  xlo = 1;
  xhi = n;
//
//  Use up the old value, if we have it.
//
  if ( saved == 1 )
  {
    x[0] = y;
    saved = 0;
    xlo = 2;
  }
//
//  If we don't need any more values, return.
//
  if ( xhi - xlo + 1 == 0 )
  {
    return;
  }

  r = new double[n+1];
//
//  If we need just one new value, do that here to avoid null arrays.
//
  if ( xhi - xlo + 1 == 1 )
  {
    r8vec_uniform_01 ( 2, seed, r );

    x[xhi-1] = sqrt ( -2.0 * log ( r[0] ) ) * cos ( 2.0 * pi * r[1] );
    y =        sqrt ( -2.0 * log ( r[0] ) ) * sin ( 2.0 * pi * r[1] );

    saved = 1;

    made = made + 2;
  }
//
//  If we require an even number of values, that's easy.
//
  else if ( ( ( xhi-xlo+1) % 2 ) == 0 )
  {
    m = ( xhi-xlo+1 ) / 2;

    r8vec_uniform_01 ( 2*m, seed, r );

    for ( i = 0; i < 2*m; i = i + 2 )
    {
      x[xlo+i-1] = sqrt ( -2.0 * log ( r[i] ) ) * cos ( 2.0 * pi * r[i+1] );
      x[xlo+i]   = sqrt ( -2.0 * log ( r[i] ) ) * sin ( 2.0 * pi * r[i+1] );
    }

    made = made + xhi - xlo + 1;
  }
//
//  If we require an odd number of values, we generate an even number,
//  and handle the last pair specially, storing one in X(N), and
//  saving the other for later.
//
  else
  {
    xhi = xhi - 1;

    m = ( xhi-xlo+1 ) / 2 + 1;

    r8vec_uniform_01 ( 2*m, seed, r );

    for ( i = 0; i < 2*m-2; i = i + 2 )
    {
      x[xlo+i-1] = sqrt ( -2.0 * log ( r[i] ) ) * cos ( 2.0 * pi * r[i+1] );
      x[xlo+i  ] = sqrt ( -2.0 * log ( r[i] ) ) * sin ( 2.0 * pi * r[i+1] );
    }

    x[n-1] = sqrt ( -2.0 * log ( r[2*m-2] ) ) * cos ( 2.0 * pi * r[2*m-1] );
    y =      sqrt ( -2.0 * log ( r[2*m-2] ) ) * sin ( 2.0 * pi * r[2*m-1] );

    saved = 1;

    made = made + xhi - xlo + 2;

  }

  delete [] r;

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

void r8vec_transpose_print ( int n, double x[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_TRANSPOSE_PRINT prints a vector on one line.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 November 2010
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the dimension of the vector.
//
//    Input, double X[N], the vector.
//
//    Input, string TITLE, a title.
//
{
  int i;

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

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

void r8vec_uniform_01 ( int n, int *seed, double r[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_UNIFORM_01 fills a double precision vector with pseudorandom values.
//
//  Discussion:
//
//    This routine implements the recursion
//
//      seed = 16807 * seed mod ( 2**31 - 1 )
//      unif = seed / ( 2**31 - 1 )
//
//    The integer arithmetic never requires more than 32 bits,
//    including a sign bit.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Paul Bratley, Bennett Fox, L E Schrage,
//    A Guide to Simulation,
//    Springer Verlag, pages 201-202, 1983.
//
//    Bennett Fox,
//    Algorithm 647:
//    Implementation and Relative Efficiency of Quasirandom
//    Sequence Generators,
//    ACM Transactions on Mathematical Software,
//    Volume 12, Number 4, pages 362-376, 1986.
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Input/output, int *SEED, a seed for the random number generator.
//
//    Output, double R[N], the vector of pseudorandom values.
//
{
  int i;
  int k;

  for ( i = 0; i < n; i++ )
  {
    k = *seed / 127773;

    *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;

    if ( *seed < 0 )
    {
      *seed = *seed + 2147483647;
    }

    r[i] = ( double ) ( *seed ) * 4.656612875E-10;
  }

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

double *r8vec_uniform_01_new ( int n, int *seed )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_UNIFORM_01_NEW returns a new unit pseudorandom R8VEC.
//
//  Discussion:
//
//    This routine implements the recursion
//
//      seed = ( 16807 * seed ) mod ( 2^31 - 1 )
//      u = seed / ( 2^31 - 1 )
//
//    The integer arithmetic never requires more than 32 bits,
//    including a sign bit.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Paul Bratley, Bennett Fox, Linus Schrage,
//    A Guide to Simulation,
//    Second Edition,
//    Springer, 1987,
//    ISBN: 0387964673,
//    LC: QA76.9.C65.B73.
//
//    Bennett Fox,
//    Algorithm 647:
//    Implementation and Relative Efficiency of Quasirandom
//    Sequence Generators,
//    ACM Transactions on Mathematical Software,
//    Volume 12, Number 4, December 1986, pages 362-376.
//
//    Pierre L'Ecuyer,
//    Random Number Generation,
//    in Handbook of Simulation,
//    edited by Jerry Banks,
//    Wiley, 1998,
//    ISBN: 0471134031,
//    LC: T57.62.H37.
//
//    Peter Lewis, Allen Goodman, James Miller,
//    A Pseudo-Random Number Generator for the System/360,
//    IBM Systems Journal,
//    Volume 8, Number 2, 1969, pages 136-143.
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Input/output, int *SEED, a seed for the random number generator.
//
//    Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values.
//
{
  int i;
  int i4_huge = 2147483647;
  int k;
  double *r;

  if ( *seed == 0 )
  {
    cerr << "\n";
    cerr << "R8VEC_UNIFORM_01_NEW - Fatal error!\n";
    cerr << "  Input value of SEED = 0.\n";
    exit ( 1 );
  }

  r = new double[n];

  for ( i = 0; i < n; i++ )
  {
    k = *seed / 127773;

    *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;

    if ( *seed < 0 )
    {
      *seed = *seed + i4_huge;
    }

    r[i] = ( double ) ( *seed ) * 4.656612875E-10;
  }

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

void r8vec_zero ( int n, double a[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_ZERO zeroes an R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    03 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Output, double A[N], a vector of zeroes.
//
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    a[i] = 0.0;
  }
  return;
}
//****************************************************************************80

double *sphere_stereograph ( int m, int n, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    SPHERE_STEREOGRAPH computes the stereographic image of points on a sphere.
//
//  Discussion:
//
//    We start with a sphere of radius 1 and center (0,0,0).
//
//    The north pole N = (0,0,1) is the point of tangency to the sphere
//    of a plane, and the south pole S = (0,0,-1) is the focus for the
//    stereographic projection.
//
//    For any point P on the sphere, the stereographic projection Q of the
//    point is defined by drawing the line from S through P, and computing
//    Q as the intersection of this line with the plane.
//
//    Actually, we allow the spatial dimension M to be arbitrary.  Values
//    of M make sense starting with 2.  The north and south poles are
//    selected as the points (0,0,...,+1) and (0,0,...,-1).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 November 2010
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    C F Marcus,
//    The stereographic projection in vector notation,
//    Mathematics Magazine,
//    Volume 39, Number 2, March 1966, pages 100-102.
//
//  Parameters:
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, double P[M*N], a set of points on the unit sphere.
//
//    Output, double SPHERE_STEREOGRAPH[M*N], the coordinates of the
//    image points.
//
{
  int i;
  int j;
  double *q;

  q = new double[m*n];

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m - 1; i++ )
    {
      q[i+j*m] = 2.0 * p[i+j*m] / ( 1.0 + p[m-1+j*m] );
    }
    q[m-1+j*m] = 1.0;
  }

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

double *sphere_stereograph_inverse ( int m, int n, double q[] )

//****************************************************************************80
//
//  Purpose:
//
//    SPHERE_STEREOGRAPH_INVERSE computes stereographic preimages of points.
//
//  Discussion:
//
//    We start with a sphere of radius 1 and center (0,0,0).
//
//    The north pole N = (0,0,1) is the point of tangency to the sphere
//    of a plane, and the south pole S = (0,0,-1) is the focus for the
//    stereographic projection.
//
//    For any point Q on the plane, the stereographic inverse projection
//    P of the point is defined by drawing the line from S through Q, and
//    computing P as the intersection of this line with the sphere.
//
//    Actually, we allow the spatial dimension M to be arbitrary.  Values
//    of M make sense starting with 2.  The north and south poles are
//    selected as the points (0,0,...,+1) and (0,0,...,-1).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 November 2010
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    C F Marcus,
//    The stereographic projection in vector notation,
//    Mathematics Magazine,
//    Volume 39, Number 2, March 1966, pages 100-102.
//
//  Parameters:
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, double Q[M*N], the points, which are presumed to lie
//    on the plane Z = 1.
//
//    Output, double SPHERE_STEREOGRAPH_INVERSE[M*N], the stereographic
//    inverse projections of the points.
//
{
  int i;
  int j;
  double *p;
  double qn;

  p = new double[m*n];

  for ( j = 0; j < n; j++ )
  {
    qn = 0.0;
    for ( i = 0; i < m - 1; i++ )
    {
      qn = qn + pow ( q[i+j*m], 2 );
    }
    for ( i = 0; i < m - 1; i++ )
    {
      p[i+j*m] = 4.0 * q[i+j*m] / ( 4.0 + qn );
    }
    p[m-1+j*m] = ( 4.0 - qn ) / ( 4.0 + qn );
  }

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

double *sphere_stereograph2 ( int m, int n, double p[], double focus[],
  double center[] )

//****************************************************************************80
//
//  Purpose:
//
//    SPHERE_STEREOGRAPH2 computes the stereographic image of points on a sphere.
//
//  Discussion:
//
//    We start with a sphere of center C.
//
//    F is a point on the sphere which is the focus of the mapping,
//    and the antipodal point 2*C-F is the point of tangency
//    to the sphere of a plane.
//
//    For any point P on the sphere, the stereographic projection Q of the
//    point is defined by drawing the line from F through P, and computing
//    Q as the intersection of this line with the plane.
//
//    The spatial dimension M is arbitrary, but should be at least 2.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 November 2010
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    C F Marcus,
//    The stereographic projection in vector notation,
//    Mathematics Magazine,
//    Volume 39, Number 2, March 1966, pages 100-102.
//
//  Parameters:
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, double P[M*N], a set of points on the unit sphere.
//
//    Input, double FOCUS[M], the coordinates of the focus point.
//
//    Input, double CENTER[M], the coordinates of the center of the sphere.
//
//    Output, double SPHERE_STEREOGRAPH2[M*N], the coordinates of the
//    image points,
//
{
  double cf_dot_pf;
  double cf_normsq;
  int i;
  int j;
  double *q;
  double s;

  q = new double[m*n];

  for ( j = 0; j < n; j++ )
  {
    cf_normsq = 0.0;
    cf_dot_pf = 0.0;
    for ( i = 0; i < m; i++ )
    {
      cf_normsq = cf_normsq + pow ( center[i] - focus[i], 2 );
      cf_dot_pf = cf_dot_pf + ( center[i] - focus[i] ) * ( p[i+j*m] - focus[i] );
    }
    s = 2.0 * cf_normsq / cf_dot_pf;
    for ( i = 0; i < m; i++ )
    {
      q[i+j*m] = s * p[i+j*m] + ( 1.0 - s ) * focus[i];
    }
  }
  return q;
}
//****************************************************************************80

double *sphere_stereograph2_inverse ( int m, int n, double q[], double focus[],
  double center[] )

//****************************************************************************80
//
//  Purpose:
//
//    SPHERE_STEREOGRAPH2_INVERSE computes stereographic preimages of points.
//
//  Discussion:
//
//    We start with a sphere of center C.
//
//    F is a point on the sphere which is the focus of the mapping,
//    and the antipodal point 2*C-F is the point of tangency
//    to the sphere of a plane.
//
//    For any point Q on the plane, the stereographic inverse projection
//    P of the point is defined by drawing the line from F through Q, and
//    computing P as the intersection of this line with the sphere.
//
//    The spatial dimension M is arbitrary, but should be at least 2.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 November 2010
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    C F Marcus,
//    The stereographic projection in vector notation,
//    Mathematics Magazine,
//    Volume 39, Number 2, March 1966, pages 100-102.
//
//  Parameters:
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, double Q[M*N], the points, which are presumed to lie
//    on the plane.
//
//    Input, double FOCUS[M], the coordinates of the focus point.
//
//    Input, double CENTER[M], the coordinates of the center of the sphere.
//
//    Output, double SPHERE_STEREOGRAPH2_INVERSE[M*N], the stereographic
//    inverse projections of the points.
//
{
  double cf_dot_qf;
  int i;
  int j;
  double *p;
  double qf_normsq;
  double s;

  p = new double[m*n];

  for ( j = 0; j < n; j++ )
  {
    cf_dot_qf = 0.0;
    qf_normsq = 0.0;
    for ( i = 0; i < m; i++ )
    {
      cf_dot_qf = cf_dot_qf + ( center[i] - focus[i] ) * ( q[i+j*m] - focus[i] );
      qf_normsq = qf_normsq + pow ( q[i+j*m] - focus[i], 2 );
    }
    s = 2.0 * cf_dot_qf / qf_normsq;
    for ( i = 0; i < m; i++ )
    {
      p[i+j*m] = s * q[i+j*m] + ( 1.0 - s ) * focus[i];
    }
  }
  return p;
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 July 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

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

double *uniform_on_sphere01_map ( int dim_num, int n, int *seed )

//****************************************************************************80
//
//  Purpose:
//
//    UNIFORM_ON_SPHERE01_MAP maps uniform points onto the unit sphere.
//
//  Discussion:
//
//    The sphere has center 0 and radius 1.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Russell Cheng,
//    Random Variate Generation,
//    in Handbook of Simulation,
//    edited by Jerry Banks,
//    Wiley, 1998, pages 168.
//
//    Reuven Rubinstein,
//    Monte Carlo Optimization, Simulation, and Sensitivity
//    of Queueing Networks,
//    Krieger, 1992,
//    ISBN: 0894647644,
//    LC: QA298.R79.
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the space.
//
//    Input, int N, the number of points.
//
//    Input/output, int *SEED, a seed for the random number generator.
//
//    Output, double UNIFORM_ON_SPHERE01_MAP[DIM_NUM*N], the points.
//
{
  int i;
  int j;
  double norm;
  double *u;
  double *x;

  u = new double[dim_num];
  x = new double[dim_num*n];

  for ( j = 0; j < n; j++ )
  {
//
//  Fill a vector with normally distributed values.
//
    r8vec_normal_01 ( dim_num, seed, u );
//
//  Compute the length of the vector.
//
    norm = 0.0;
    for ( i = 0; i < dim_num; i++ )
    {
      norm = norm + u[i] * u[i];
    }
    norm = sqrt ( norm );
//
//  Normalize the vector.
//
    for ( i = 0; i < dim_num; i++ )
    {
      x[i+j*dim_num] = u[i] / norm;
    }

  }

  delete [] u;
  return x;
}