subroutine random_initialize ( seed )

!*****************************************************************************80
!
!! RANDOM_INITIALIZE initializes the FORTRAN90 random number seed.
!
!  Discussion:
!
!    If you don't initialize the random number generator, its behavior
!    is not specified.  If you initialize it simply by:
!
!      call random_seed ( )
!
!    its behavior is not specified.  On the DEC ALPHA, if that's all you
!    do, the same random number sequence is returned.  In order to actually
!    try to scramble up the random number generator a bit, this routine
!    goes through the tedious process of getting the size of the random
!    number seed, making up values based on the current time, and setting
!    the random number seed.
!
!  Modified:
!
!    19 December 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer SEED.
!    If SEED is zero on input, then you're asking this routine to come up
!    with a seed value, which is returned as output.
!    If SEED is nonzero on input, then you're asking this routine to
!    use the input value of SEED to initialize the random number generator,
!    and SEED is not changed on output.
!
  implicit none

  integer count
  integer count_max
  integer count_rate
  logical, parameter :: debug = .false.
  integer i
  integer seed
  integer, allocatable :: seed_vector(:)
  integer seed_size
  real ( kind = 8 ) t
!
!  Initialize the random number seed.
!
  call random_seed ( )
!
!  Determine the size of the random number seed.
!
  call random_seed ( size = seed_size )
!
!  Allocate a seed of the right size.
!
  allocate ( seed_vector(seed_size) )

  if ( seed /= 0 ) then

    if ( debug ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'RANDOM_INITIALIZE'
      write ( *, '(a,i20)' ) '  Initialize RANDOM_NUMBER, user SEED = ', seed
    end if

  else

    call system_clock ( count, count_rate, count_max )

    seed = count

    if ( debug ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'RANDOM_INITIALIZE'
      write ( *, '(a,i20)' ) '  Initialize RANDOM_NUMBER, arbitrary SEED = ', &
        seed
    end if

  end if
!
!  Now set the seed.
!
  seed_vector(1:seed_size) = seed

  call random_seed ( put = seed_vector(1:seed_size) )
!
!  Free up the seed space.
!
  deallocate ( seed_vector )
!
!  Call the random number routine a bunch of times.
!
  do i = 1, 100
    call random_number ( harvest = t )
  end do

  return
end
subroutine rtp_to_xyz ( r, theta, phi, v )

!*****************************************************************************80
!
!! RTP_TO_XYZ converts spherical coordinates to XYZ coordinates.
!
!  Modified:
!
!    18 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, THETA, PHI, the spherical coordinates 
!    of a point.
!
!    Output, ( kind = 8 ) real V(3), the XYZ coordinates.
!
  implicit none

  real ( kind = 8 ) phi
  real ( kind = 8 ) r
  real ( kind = 8 ) theta
  real ( kind = 8 ) v(3)

  v(1) = r * cos ( theta ) * sin ( phi )
  v(2) = r * sin ( theta ) * sin ( phi )
  v(3) = r *                 cos ( phi )

  return
end
subroutine s_cat ( s1, s2, s3 )

!*****************************************************************************80
!
!! S_CAT concatenates two strings to make a third string.
!
!  Modified:
!
!    18 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S1, the "prefix" string.
!
!    Input, character ( len = * ) S2, the "postfix" string.
!
!    Output, character ( len = * ) S3, the string made by
!    concatenating S1 and S2, ignoring any trailing blanks.
!
  implicit none

  character ( len = * ) s1
  character ( len = * ) s2
  character ( len = * ) s3

  if ( s1 == ' ' .and. s2 == ' ' ) then
    s3 = ' '
  else if ( s1 == ' ' ) then
    s3 = s2
  else if ( s2 == ' ' ) then
    s3 = s1
  else
    s3 = trim ( s1 ) // trim ( s2 )
  end if

  return
end
subroutine sphere_area_3d ( r, area )

!*****************************************************************************80
!
!! SPHERE_AREA_3D computes the surface area of a sphere in 3D.
!
!  Formula:
!
!    A sphere in 3D satisfies the equation:
!
!      ( X - XC )**2 + ( Y - YC )**2 + ( Z - ZC )**2 = R**2
!
!  Modified:
!
!    26 January 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Output, real ( kind = 8 ) AREA, the area of the sphere.
!
  implicit none

  real ( kind = 8 ) area
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) r

  area = 4.0D+00 * pi * r * r

  return
end
subroutine sphere_quad_00 ( r, f, h, result )

!*****************************************************************************80
!
!! SPHERE_QUAD_00 approximates an integral over a sphere in 3D.
!
!  Discussion:
!
!    A number of points N are chosen at random on the sphere, with N
!    being determined so that, if the points were laid out on a regular
!    grid, the average spacing would be no more than H.
!
!  Modified:
!
!    20 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real ( kind = 8 ), external :: F, evaluates the function 
!    to be integrated, of the form:
!      function f ( v )
!      real ( kind = 8 ) f
!      real ( kind = 8 ) v(3)
!
!    Input, real ( kind = 8 ) H, the maximum length of a side of the spherical
!    quadrilaterals.
!
!    Output, real ( kind = 8 ) RESULT, the approximate integral.
!
  implicit none

  logical, parameter :: debug = .false.
  real ( kind = 8 ), external :: f
  real ( kind = 8 ) h
  integer i
  integer n
  integer phi_num
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) quad
  real ( kind = 8 ) r
  real ( kind = 8 ) result
  real ( kind = 8 ) sphere_area
  integer theta_num
  real ( kind = 8 ) v(3)
!
!  Choose PHI and THETA counts that make short sides.
!
  phi_num = int ( pi / h )

  if ( h * real ( phi_num, kind = 8 ) < pi ) then
    phi_num = phi_num + 1
  end if

  theta_num = int ( 2.0D+00 * pi / h )

  if ( h * real ( theta_num, kind = 8 ) < pi ) then
    theta_num = theta_num + 1
  end if

  n = phi_num * theta_num

  if ( debug ) then
    write ( *, * ) ' '
    write ( *, * ) '  PHI_NUM =   ', phi_num
    write ( *, * ) '  THETA_NUM = ', theta_num
    write ( *, * ) '  N =         ', n
  end if

  call sphere_area_3d ( r, sphere_area )

  quad = 0.0

  do i = 1, n

    call sphere_sample_3d ( r, v )

    quad = quad + f ( v )

  end do

  result = quad * sphere_area / real ( n, kind = 8 )

  return
end
subroutine sphere_quad_01 ( r, f, h, result )

!*****************************************************************************80
!
!! SPHERE_QUAD_01 approximates an integral over a sphere in 3D.
!
!  Discussion:
!
!    The sphere is broken up into spherical triangles, whose sides
!    do not exceed the length H.  Then a centroid rule is used on
!    each spherical triangle.
!
!  Modified:
!
!    18 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real ( kind = 8 ) , external :: F, evaluates the 
!    function to be integrated, of the form:
!      function f ( v )
!      real ( kind = 8 ) f
!      real ( kind = 8 ) v(3)
!
!    Input, real ( kind = 8 ) H, the maximum length of a side of the spherical
!    quadrilaterals.
!
!    Output, real ( kind = 8 ) RESULT, the approximate integral.
!
  implicit none

  real ( kind = 8 ) area
  logical, parameter :: debug = .false.
  real ( kind = 8 ), external :: f
  real ( kind = 8 ) h
  integer i
  integer j
  real ( kind = 8 ) phi
  integer phi_num
  real ( kind = 8 ) phi1
  real ( kind = 8 ) phi2
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) quad
  real ( kind = 8 ) r
  real ( kind = 8 ) result
  real ( kind = 8 ) sector_area
  real ( kind = 8 ) sphere_area
  real ( kind = 8 ) theta
  integer theta_num
  real ( kind = 8 ) theta1
  real ( kind = 8 ) theta2
  real ( kind = 8 ) v(3)
  real ( kind = 8 ) v1(3)
  real ( kind = 8 ) v11(3)
  real ( kind = 8 ) v12(3)
  real ( kind = 8 ) v2(3)
  real ( kind = 8 ) v21(3)
  real ( kind = 8 ) v22(3)
!
!  Choose PHI and THETA counts that make short sides.
!
  phi_num = int ( pi / h )

  if ( h * real ( phi_num, kind = 8 ) < pi ) then
    phi_num = phi_num + 1
  end if

  theta_num = int ( 2.0D+00 * pi / h )

  if ( h * real ( theta_num, kind = 8 ) < pi ) then
    theta_num = theta_num + 1
  end if

  if ( debug ) then
    write ( *, * ) ' '
    write ( *, * ) '  PHI_NUM =   ', phi_num
    write ( *, * ) '  THETA_NUM = ', theta_num
  end if

  result = 0.0D+00
!
!  Only one THETA (and hence, only one PHI.)
!
  if ( theta_num == 1 ) then

    call sphere_area_3d ( r, sphere_area )

    theta = 0.0D+00
    phi = pi / 2.0D+00
    call rtp_to_xyz ( r, theta, phi, v )

    quad = f ( v )
    result = quad * sphere_area
!
!  Several THETA's, one PHI.
!
  else if ( phi_num == 1 ) then

    call sphere_area_3d ( r, sphere_area )
    sector_area = sphere_area / real ( theta_num, kind = 8 )

    result = 0.0D+00

    do j = 1, theta_num

      theta = real ( ( j - 1 ) * 2, kind = 8 ) * pi &
            / real ( theta_num, kind = 8 )
      phi = pi / 2.0D+00
      call rtp_to_xyz ( r, theta, phi, v )

      quad = f ( v )

      result = result + quad * sector_area

    end do
!
!  At least two PHI's.
!
  else

    result = 0.0
!
!  Picture:
!
!        V1
!       /  \
!      /    \
!    V12----V22
!
    phi1 = 0.0D+00
    phi2 = pi / real ( phi_num, kind = 8 )

    do j = 1, theta_num

      theta1 = real ( j - 1, kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )
      theta2 = real ( j    , kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )

      call rtp_to_xyz ( r, theta1, phi1, v1 )
      call rtp_to_xyz ( r, theta1, phi2, v12 )
      call rtp_to_xyz ( r, theta2, phi2, v22 )

      call stri_vertices_to_area_3d ( r, v1, v12, v22, area )

      call stri_vertices_to_centroid_3d ( r, v1, v12, v22, v )

      result = result + f ( v ) * area

    end do
!
!  Picture:
!
!    V11--V21
!     | \  |
!     |  \ |
!    V12--V22
!
    do i = 2, phi_num-1

      phi1 = real ( i - 1, kind = 8 ) * pi / real ( phi_num, kind = 8 )
      phi2 = real ( i,     kind = 8 ) * pi / real ( phi_num, kind = 8 )

      do j = 1, theta_num

        theta1 = real ( j - 1, kind = 8 ) * 2.0D+00 * pi &
               / real ( theta_num, kind = 8 )
        theta2 = real ( j,     kind = 8 ) * 2.0D+00 * pi &
               / real ( theta_num, kind = 8 )

        call rtp_to_xyz ( r, theta1, phi1, v11 )
        call rtp_to_xyz ( r, theta2, phi1, v21 )
        call rtp_to_xyz ( r, theta1, phi2, v12 )
        call rtp_to_xyz ( r, theta2, phi2, v22 )

        call stri_vertices_to_area_3d ( r, v11, v12, v22, area )

        call stri_vertices_to_centroid_3d ( r, v11, v12, v22, v )

        result = result + f ( v ) * area

        call stri_vertices_to_area_3d ( r, v22, v21, v11, area )

        call stri_vertices_to_centroid_3d ( r, v22, v21, v11, v )

        result = result + f ( v ) * area

      end do

    end do
!
!  Picture:
!
!    V11----V21
!      \    /
!       \  /
!        V2
!
    phi1 = real ( phi_num - 1, kind = 8 ) * pi &
         / real ( phi_num, kind = 8 )
    phi2 =                                  pi

    do j = 1, theta_num

      theta1 = real ( j - 1, kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )
      theta2 = real ( j,     kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )

      call rtp_to_xyz ( r, theta1, phi1, v11 )
      call rtp_to_xyz ( r, theta2, phi1, v21 )
      call rtp_to_xyz ( r, theta2, phi2, v2 )

      call stri_vertices_to_area_3d ( r, v11, v2, v21, area )

      call stri_vertices_to_centroid_3d ( r, v11, v2, v21, v )

      result = result + f ( v ) * area

    end do

  end if

  return
end
subroutine sphere_quad_02 ( r, f, h, result )

!*****************************************************************************80
!
!! SPHERE_QUAD_02 approximates an integral over a sphere in 3D.
!
!  Discussion:
!
!    The sphere is broken up into spherical triangles, whose sides
!    do not exceed the length H.  Then the function is evaluated
!    at the vertices, and the average is multiplied by the area.
!
!  Modified:
!
!    20 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real ( kind = 8 ) , external :: F, evaluates the 
!    function to be integrated, of the form:
!      function f ( v )
!      real ( kind = 8 ) f
!      real ( kind = 8 ) v(3)
!
!    Input, real ( kind = 8 ) H, the maximum length of a side of the spherical
!    quadrilaterals.
!
!    Output, real ( kind = 8 ) RESULT, the approximate integral.
!
  implicit none

  real ( kind = 8 ) area
  logical, parameter :: debug = .false.
  real ( kind = 8 ), external :: f
  real ( kind = 8 ) h
  integer i
  integer j
  real ( kind = 8 ) phi
  integer phi_num
  real ( kind = 8 ) phi1
  real ( kind = 8 ) phi2
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) quad
  real ( kind = 8 ) r
  real ( kind = 8 ) result
  real ( kind = 8 ) sector_area
  real ( kind = 8 ) sphere_area
  real ( kind = 8 ) theta
  integer theta_num
  real ( kind = 8 ) theta1
  real ( kind = 8 ) theta2
  real ( kind = 8 ) v(3)
  real ( kind = 8 ) v1(3)
  real ( kind = 8 ) v11(3)
  real ( kind = 8 ) v12(3)
  real ( kind = 8 ) v2(3)
  real ( kind = 8 ) v21(3)
  real ( kind = 8 ) v22(3)

  if ( debug ) then
    write ( *, * ) ' '
    write ( *, * ) 'SPHERE_QUAD_02 - Enter.'
  end if
!
!  Choose PHI and THETA counts that make short sides.
!
  phi_num = int ( pi / h )

  if ( h * real ( phi_num, kind = 8 ) < pi ) then
    phi_num = phi_num + 1
  end if

  theta_num = int ( 2.0D+00 * pi / h )

  if ( h * real ( theta_num, kind = 8 ) < pi ) then
    theta_num = theta_num + 1
  end if

  if ( debug ) then
    write ( *, * ) ' '
    write ( *, * ) '  PHI_NUM =   ', phi_num
    write ( *, * ) '  THETA_NUM = ', theta_num
  end if

  result = 0.0D+00
!
!  Only one THETA (and hence, only one PHI.)
!
  if ( theta_num == 1 ) then

    call sphere_area_3d ( r, sphere_area )

    theta = 0.0D+00
    phi = pi / 2.0D+00
    call rtp_to_xyz ( r, theta, phi, v )

    quad = f ( v )
    result = quad * sphere_area
!
!  Several THETA's, one PHI.
!
  else if ( phi_num == 1 ) then

    call sphere_area_3d ( r, sphere_area )
    sector_area = sphere_area / real ( theta_num, kind = 8 )

    result = 0.0D+00

    do j = 1, theta_num

      theta = real ( ( j - 1 ) * 2, kind = 8 ) * pi &
        / real ( theta_num, kind = 8 )
      phi = pi / 2.0D+00
      call rtp_to_xyz ( r, theta, phi, v )

      quad = f ( v )

      result = result + quad * sector_area

    end do
!
!  At least two PHI's.
!
  else

    result = 0.0D+00
!
!  Picture:
!
!        V1
!       /  \
!      /    \
!    V12----V22
!
    phi1 = 0.0D+00
    phi2 = pi / real ( phi_num, kind = 8 )

    if ( debug ) then
      write ( *, * ) 'PHI1 = ', phi1
      write ( *, * ) 'PHI2 = ', phi2
    end if

    do j = 1, theta_num

      theta1 = real ( j - 1, kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )
      theta2 = real ( j,     kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )

      if ( debug ) then
        write ( *, * ) 'THETA1 = ', theta1
        write ( *, * ) 'THETA2 = ', theta2
      end if

      call rtp_to_xyz ( r, theta1, phi1, v1 )
      call rtp_to_xyz ( r, theta1, phi2, v12 )
      call rtp_to_xyz ( r, theta2, phi2, v22 )

      if ( debug ) then
        write ( *, * ) 'V1 =  ', v1(1:3)
        write ( *, * ) 'V12 = ', v12(1:3)
        write ( *, * ) 'V22 = ', v22(1:3)
      end if

      call stri_vertices_to_area_3d ( r, v1, v12, v22, area )

      if ( debug ) then
        write ( *, * ) 'AREA = ', area
      end if

      result = result + ( f ( v1 ) + f ( v12 ) + f ( v22 ) ) * area &
        / 3.0D+00

    end do
!
!  Picture:
!
!    V11--V21
!     | \  |
!     |  \ |
!    V12--V22
!
    do i = 2, phi_num-1

      phi1 = real ( i - 1, kind = 8 ) * pi &
           / real ( phi_num, kind = 8 )
      phi2 = real ( i,     kind = 8 ) * pi &
           / real ( phi_num, kind = 8 )

      do j = 1, theta_num

        theta1 = real ( j - 1, kind = 8 ) * 2.0D+00 * pi &
               / real ( theta_num, kind = 8 )
        theta2 = real ( j,     kind = 8 ) * 2.0D+00 * pi &
               / real ( theta_num, kind = 8 )

        call rtp_to_xyz ( r, theta1, phi1, v11 )
        call rtp_to_xyz ( r, theta2, phi1, v21 )
        call rtp_to_xyz ( r, theta1, phi2, v12 )
        call rtp_to_xyz ( r, theta2, phi2, v22 )

        call stri_vertices_to_area_3d ( r, v11, v12, v22, area )

        result = result + ( f ( v11 ) + f ( v12 ) + f ( v22 ) ) * area &
          / 3.0D+00

        call stri_vertices_to_area_3d ( r, v22, v21, v11, area )

        result = result + ( f ( v22 ) + f ( v21 ) + f ( v11 ) ) * area &
          / 3.0D+00

      end do

    end do
!
!  Picture:
!
!    V11----V21
!      \    /
!       \  /
!        V2
!
    phi1 = real ( phi_num - 1, kind = 8 ) * pi &
         / real ( phi_num, kind = 8 )
    phi2 =                                  pi

    do j = 1, theta_num

      theta1 = real ( j - 1, kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )
      theta2 = real ( j,     kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )

      call rtp_to_xyz ( r, theta1, phi1, v11 )
      call rtp_to_xyz ( r, theta2, phi1, v21 )
      call rtp_to_xyz ( r, theta2, phi2, v2 )

      call stri_vertices_to_area_3d ( r, v11, v2, v21, area )

      result = result + ( f ( v11 ) + f ( v2 ) + f ( v21 ) ) * area &
        / 3.0D+00

    end do

  end if

  if ( debug ) then
    write ( *, * ) ' '
    write ( *, * ) 'SPHERE_QUAD_02 - Exit.'
  end if

  return
end
subroutine sphere_quad_03 ( r, f, h, result )

!*****************************************************************************80
!
!! SPHERE_QUAD_03 approximates an integral over a sphere in 3D.
!
!  Discussion:
!
!    The sphere is broken up into spherical triangles, whose sides
!    do not exceed the length H.  Then the function is evaluated
!    at the midpoints, and the average is multiplied by the area.
!
!  Modified:
!
!    20 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real( kind = 8 ) , external :: F, evaluates the function 
!    to be integrated, of the form:
!      function f ( v )
!      real ( kind = 8 ) f
!      real ( kind = 8 ) v(3)
!
!    Input, real ( kind = 8 ) H, the maximum length of a side of the spherical
!    quadrilaterals.
!
!    Output, real ( kind = 8 ) RESULT, the approximate integral.
!
  implicit none

  real ( kind = 8 ) area
  logical, parameter :: debug = .false.
  real ( kind = 8 ), external :: f
  real ( kind = 8 ) h
  integer i
  integer j
  real ( kind = 8 ) m1(3)
  real ( kind = 8 ) m2(3)
  real ( kind = 8 ) m3(3)
  real ( kind = 8 ) phi
  integer phi_num
  real ( kind = 8 ) phi1
  real ( kind = 8 ) phi2
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) quad
  real ( kind = 8 ) r
  real ( kind = 8 ) result
  real ( kind = 8 ) sector_area
  real ( kind = 8 ) sphere_area
  real ( kind = 8 ) theta
  integer theta_num
  real ( kind = 8 ) theta1
  real ( kind = 8 ) theta2
  real ( kind = 8 ) v(3)
  real ( kind = 8 ) v1(3)
  real ( kind = 8 ) v11(3)
  real ( kind = 8 ) v12(3)
  real ( kind = 8 ) v2(3)
  real ( kind = 8 ) v21(3)
  real ( kind = 8 ) v22(3)
!
!  Choose PHI and THETA counts that make short sides.
!
  phi_num = int ( pi / h )

  if ( h * real ( phi_num, kind = 8 ) < pi ) then
    phi_num = phi_num + 1
  end if

  theta_num = int ( 2.0D+00 * pi / h )

  if ( h * real ( theta_num, kind = 8 ) < pi ) then
    theta_num = theta_num + 1
  end if

  if ( debug ) then
    write ( *, * ) ' '
    write ( *, * ) '  PHI_NUM =   ', phi_num
    write ( *, * ) '  THETA_NUM = ', theta_num
  end if

  result = 0.0D+00
!
!  Only one THETA (and hence, only one PHI.)
!
  if ( theta_num == 1 ) then

    call sphere_area_3d ( r, sphere_area )

    theta = 0.0D+00
    phi = pi / 2.0D+00
    call rtp_to_xyz ( r, theta, phi, v )

    quad = f ( v )
    result = quad * sphere_area
!
!  Several THETA's, one PHI.
!
  else if ( phi_num == 1 ) then

    call sphere_area_3d ( r, sphere_area )
    sector_area = sphere_area / real ( theta_num, kind = 8 )

    result = 0.0D+00

    do j = 1, theta_num

      theta = real ( ( j - 1 ) * 2, kind = 8 ) * pi &
            / real ( theta_num, kind = 8 )
      phi = pi / 2.0D+00
      call rtp_to_xyz ( r, theta, phi, v )

      quad = f ( v )

      result = result + quad * sector_area

    end do
!
!  At least two PHI's.
!
  else

    result = 0.0D+00
!
!  Picture:
!
!        V1
!       /  \
!      /    \
!    V12----V22
!
    phi1 = 0.0D+00
    phi2 = pi / real ( phi_num, kind = 8 )

    do j = 1, theta_num

      theta1 = real ( j - 1, kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )
      theta2 = real ( j,     kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )

      call rtp_to_xyz ( r, theta1, phi1, v1 )
      call rtp_to_xyz ( r, theta1, phi2, v12 )
      call rtp_to_xyz ( r, theta2, phi2, v22 )

      call stri_vertices_to_area_3d ( r, v1, v12, v22, area )

      call stri_vertices_to_midpoints_3d ( r, v1, v12, v22, m1, m2, m3 )

      result = result + ( f ( m1 ) + f ( m2 ) + f ( m3 ) ) * area &
        / 3.0D+00

    end do
!
!  Picture:
!
!    V11--V21
!     | \  |
!     |  \ |
!    V12--V22
!
    do i = 2, phi_num-1

      phi1 = real ( i - 1, kind = 8 ) * pi &
           / real ( phi_num, kind = 8 )
      phi2 = real ( i,     kind = 8 ) * pi &
           / real ( phi_num, kind = 8 )

      do j = 1, theta_num

        theta1 = real ( j - 1, kind = 8 ) * 2.0D+00 * pi &
               / real ( theta_num, kind = 8 )
        theta2 = real ( j,     kind = 8 ) * 2.0D+00 * pi &
               / real ( theta_num, kind = 8 )

        call rtp_to_xyz ( r, theta1, phi1, v11 )
        call rtp_to_xyz ( r, theta2, phi1, v21 )
        call rtp_to_xyz ( r, theta1, phi2, v12 )
        call rtp_to_xyz ( r, theta2, phi2, v22 )

        call stri_vertices_to_area_3d ( r, v11, v12, v22, area )

        call stri_vertices_to_midpoints_3d ( r, v11, v12, v22, m1, m2, m3 )

        result = result + ( f ( m1 ) + f ( m2 ) + f ( m3 ) ) * area &
          / 3.0D+00

        call stri_vertices_to_area_3d ( r, v22, v21, v11, area )

        call stri_vertices_to_midpoints_3d ( r, v22, v21, v11, m1, m2, m3 )

        result = result + ( f ( m1 ) + f ( m2 ) + f ( m3 ) ) * area &
          / 3.0D+00

      end do

    end do
!
!  Picture:
!
!    V11----V21
!      \    /
!       \  /
!        V2
!
    phi1 = real ( phi_num - 1, kind = 8 ) * pi &
         / real ( phi_num, kind = 8 )
    phi2 =                                  pi

    do j = 1, theta_num

      theta1 = real ( j - 1, kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )
      theta2 = real ( j,     kind = 8 ) * 2.0D+00 * pi &
             / real ( theta_num, kind = 8 )

      call rtp_to_xyz ( r, theta1, phi1, v11 )
      call rtp_to_xyz ( r, theta2, phi1, v21 )
      call rtp_to_xyz ( r, theta2, phi2, v2 )

      call stri_vertices_to_area_3d ( r, v11, v2, v21, area )

      call stri_vertices_to_midpoints_3d ( r, v11, v2, v21, m1, m2, m3 )

      result = result + ( f ( m1 ) + f ( m2 ) + f ( m3 ) ) * area &
        / 3.0D+00

    end do

  end if

  return
end
subroutine sphere_sample_3d ( r, x )

!*****************************************************************************80
!
!! SPHERE_SAMPLE_3D picks a random point on a sphere in 3D.
!
!  Modified:
!
!    20 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Output, real ( kind = 8 ) X(3), the sample point.
!
  implicit none

  real ( kind = 8 ) phi
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) r
  real ( kind = 8 ) theta
  real ( kind = 8 ) vdot
  real ( kind = 8 ) x(3)
!
!  Pick a uniformly random VDOT, which must be between -1 and 1.
!  This represents the dot product of the random vector with the Z unit vector.
!
!  Note: this works because the surface area of the sphere between
!  Z and Z + dZ is independent of Z.  So choosing Z uniformly chooses
!  a patch of area uniformly.
!
  call random_number ( harvest = vdot )
  vdot = 2.0D+00 * vdot - 1.0D+00

  phi = acos ( vdot )
!
!  Pick a uniformly random rotation between 0 and 2 Pi around the
!  axis of the Z vector.
!
  call random_number ( harvest = theta )
  theta = 2.0D+00 * pi * theta

  x(1) = r * cos ( theta ) * sin ( phi )
  x(2) = r * sin ( theta ) * sin ( phi )
  x(3) = r * cos ( phi )

  return
end
subroutine stri_angles_to_area_3d ( r, a, b, c, area )

!*****************************************************************************80
!
!! STRI_ANGLES_TO_AREA_3D computes the area of a spherical triangle.
!
!  Formula:
!
!    A sphere in 3D satisfies the equation:
!
!      X**2 + Y**2 + Z**2 = R**2
!
!    A spherical triangle is specified by three points on the surface
!    of the sphere.
!
!    The area formula is known as Girard's formula.
!
!    The area of a spherical triangle is:
!
!      AREA = ( A + B + C - PI ) * R**2
!
!    where A, B and C are the (surface) angles of the triangle.
!
!  Modified:
!
!    09 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real ( kind = 8 ) A, B, C, the angles of the triangle.
!
!    Output, real ( kind = 8 ) AREA, the area of the sphere.
!
  implicit none

  real ( kind = 8 ) area
  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) c
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  real ( kind = 8 ) r
!
!  Apply Girard's formula.
!
  area = r * r * ( a + b + c - pi )

  return
end
subroutine stri_quad_01 ( r, v1, v2, v3, f, result )

!*****************************************************************************80
!
!! STRI_QUAD_01 approximates an integral over a spherical triangle in 3D.
!
!  Discussion:
!
!    The integral is approximated by the value at the centroid,
!    multiplied by the area.
!
!  Modified:
!
!    18 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real ( kind = 8 ) V1(3), V2(3), V3(3), the vertices of the triangle.
!
!    Input, real ( kind = 8 ), external :: F, evaluates the function
!    to be integrated, of the form:
!      function f ( v )
!      real ( kind = 8 ) f
!      real ( kind = 8 ) v(3)
!
!    Output, real ( kind = 8 ) RESULT, the approximate integral.
!
  implicit none

  real ( kind = 8 ) area
  real ( kind = 8 ), external :: f
  real ( kind = 8 ) quad
  real ( kind = 8 ) r
  real ( kind = 8 ) result
  real ( kind = 8 ) v1(3)
  real ( kind = 8 ) v2(3)
  real ( kind = 8 ) v3(3)
  real ( kind = 8 ) vc(3)

  call stri_vertices_to_area_3d ( r, v1, v2, v3, area )

  call stri_vertices_to_centroid_3d ( r, v1, v2, v3, vc )

  quad = f ( vc )
  result = quad * area

  return
end
subroutine stri_sides_to_angles_3d ( r, as, bs, cs, a, b, c )

!*****************************************************************************80
!
!! STRI_SIDES_TO_ANGLES_3D computes spherical triangle angles in 3D.
!
!  Modified:
!
!    09 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real ( kind = 8 ) AS, BS, CS, the (geodesic) length of the 
!    sides of the triangle.
!
!    Output, real ( kind = 8 ) A, B, C, the spherical angles of the triangle.
!    Angle A is opposite the side of length AS, and so on.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) as
  real ( kind = 8 ) asu
  real ( kind = 8 ) b
  real ( kind = 8 ) bs
  real ( kind = 8 ) bsu
  real ( kind = 8 ) c
  real ( kind = 8 ) cs
  real ( kind = 8 ) csu
  real ( kind = 8 ) r
  real ( kind = 8 ) ssu
  real ( kind = 8 ) tan_a2
  real ( kind = 8 ) tan_b2
  real ( kind = 8 ) tan_c2

  asu = as / r
  bsu = bs / r
  csu = cs / r
  ssu = ( asu + bsu + csu ) / 2.0D+00

  tan_a2 = sqrt ( ( sin ( ssu - bsu ) * sin ( ssu - csu ) ) / &
                  ( sin ( ssu ) * sin ( ssu - asu )     ) )

  a = 2.0D+00 * atan ( tan_a2 )

  tan_b2 = sqrt ( ( sin ( ssu - asu ) * sin ( ssu - csu ) ) / &
                  ( sin ( ssu ) * sin ( ssu - bsu )     ) )

  b = 2.0D+00 * atan ( tan_b2 )

  tan_c2 = sqrt ( ( sin ( ssu - asu ) * sin ( ssu - bsu ) ) / &
                  ( sin ( ssu ) * sin ( ssu - csu )     ) )

  c = 2.0D+00 * atan ( tan_c2 )

  return
end
subroutine stri_vertices_to_area_3d ( r, v1, v2, v3, area )

!*****************************************************************************80
!
!! STRI_VERTICES_TO_AREA_3D computes the area of a spherical triangle in 3D.
!
!  Formula:
!
!    A sphere in 3D satisfies the equation:
!
!      X**2 + Y**2 + Z**2 = R**2
!
!    A spherical triangle is specified by three points on the surface
!    of the sphere.
!
!    The area formula is known as Girard's formula.
!
!    The area of a spherical triangle is:
!
!      AREA = ( A + B + C - PI ) * R**2
!
!    where A, B and C are the (surface) angles of the triangle.
!
!  Modified:
!
!    09 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real ( kind = 8 ) V1(3), V2(3), V3(3), the vertices of the triangle.
!
!    Output, real ( kind = 8 ) AREA, the area of the sphere.
!
  implicit none

  real ( kind = 8 ) area
  real ( kind = 8 ) a
  real ( kind = 8 ) as
  real ( kind = 8 ) b
  real ( kind = 8 ) bs
  real ( kind = 8 ) c
  real ( kind = 8 ) cs
  real ( kind = 8 ) r
  real ( kind = 8 ) v1(3)
  real ( kind = 8 ) v2(3)
  real ( kind = 8 ) v3(3)
!
!  Compute the lengths of the sides of the spherical triangle.
!
  call stri_vertices_to_sides ( r, v1, v2, v3, as, bs, cs )
!
!  Get the spherical angles.
!
  call stri_sides_to_angles_3d ( r, as, bs, cs, a, b, c )
!
!  Get the area
!
  call stri_angles_to_area_3d ( r, a, b, c, area )

  return
end
subroutine stri_vertices_to_centroid_3d ( r, v1, v2, v3, vs )

!*****************************************************************************80
!
!! STRI_VERTICES_TO_CENTROID_3D gets a spherical triangle "centroid" in 3D.
!
!  Formula:
!
!    A sphere in 3D satisfies the equation:
!
!      X**2 + Y**2 + Z**2 = R**2
!
!    A spherical triangle is specified by three points on the sphere.
!
!    The (true) centroid of a spherical triangle is the point
!
!      VT = (XT,YT,ZT) = Integral ( X, Y, Z ) dArea / Integral 1 dArea
!
!    Note that the true centroid does NOT, in general, lie on the sphere.
!
!    The "flat" centroid VF is the centroid of the planar triangle defined by
!    the vertices of the spherical triangle.
!
!    The "spherical" centroid VS of a spherical triangle is computed by
!    the intersection of the geodesic bisectors of the triangle angles.
!    The spherical centroid lies on the sphere.
!
!    VF, VT and VS lie on a line through the center of the sphere.  We can
!    easily calculate VF by averaging the vertices, and from this determine
!    VS by normalizing.
!
!    (Of course, we still will not have actually computed VT, which lies
!    somewhere between VF and VS!)
!
!  Modified:
!
!    18 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real ( kind = 8 ) V1(3), V2(3), V3(3), the vertices of the triangle.
!
!    Output, real ( kind = 8 ) VS(3), the coordinates of the "spherical
!    centroid" of the spherical triangle.
!
  implicit none

  real ( kind = 8 ) norm
  real ( kind = 8 ) r
  real ( kind = 8 ) v1(3)
  real ( kind = 8 ) v2(3)
  real ( kind = 8 ) v3(3)
  real ( kind = 8 ) vs(3)

  vs(1:3) = ( v1(1:3) + v2(1:3) + v3(1:3) ) / 3.0D+00

  norm = sqrt ( sum ( vs(1:3)**2 ) )

  vs(1:3) = r * vs(1:3) / norm

  return
end
subroutine stri_vertices_to_midpoints_3d ( r, v1, v2, v3, m1, m2, m3 )

!*****************************************************************************80
!
!! STRI_VERTICES_TO_MIDPOINTS_3D gets a spherical triangles midsides in 3D.
!
!  Modified:
!
!    20 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real ( kind = 8 ) V1(3), V2(3), V3(3), the vertices of the triangle.
!
!    Output, real ( kind = 8 ) M1(3), M2(3), M3(3), the coordinates of 
!    the midpoints of the sides of the spherical triangle.
!
  implicit none

  real ( kind = 8 ) m1(3)
  real ( kind = 8 ) m2(3)
  real ( kind = 8 ) m3(3)
  real ( kind = 8 ) norm
  real ( kind = 8 ) r
  real ( kind = 8 ) v1(3)
  real ( kind = 8 ) v2(3)
  real ( kind = 8 ) v3(3)

  m1(1:3) = ( v1(1:3) + v2(1:3) ) / 2.0D+00
  norm = sqrt ( sum ( m1(1:3)**2 ) )
  m1(1:3) = r * m1(1:3) / norm

  m2(1:3) = ( v2(1:3) + v3(1:3) ) / 2.0D+00
  norm = sqrt ( sum ( m2(1:3)**2 ) )
  m2(1:3) = r * m2(1:3) / norm

  m3(1:3) = ( v3(1:3) + v1(1:3) ) / 2.0D+00
  norm = sqrt ( sum ( m3(1:3)**2 ) )
  m3(1:3) = r * m3(1:3) / norm

  return
end
subroutine stri_vertices_to_sides ( r, v1, v2, v3, as, bs, cs )

!*****************************************************************************80
!
!! STRI_VERTICES_TO_SIDES_3D computes spherical triangle sides in 3D.
!
!  Modified:
!
!    09 June 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) R, the radius of the sphere.
!
!    Input, real ( kind = 8 ) V1(3), V2(3), V3(3), the vertices of the spherical
!    triangle.
!
!    Output, real ( kind = 8 ) AS, BS, CS, the (geodesic) length of the 
!    sides of the triangle.
!
  implicit none

  real ( kind = 8 ) as
  real ( kind = 8 ) bs
  real ( kind = 8 ) cs
  real ( kind = 8 ) r
  real ( kind = 8 ) v1(3)
  real ( kind = 8 ) v2(3)
  real ( kind = 8 ) v3(3)

  as = r * acos ( dot_product ( v2(1:3), v3(1:3) ) / r**2 )
  bs = r * acos ( dot_product ( v3(1:3), v1(1:3) ) / r**2 )
  cs = r * acos ( dot_product ( v1(1:3), v2(1:3) ) / r**2 )

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end