subroutine r8_fake_use ( x )

!*****************************************************************************80
!
!! r8_fake_use() pretends to use an R8 variable.
!
!  Discussion:
!
!    Some compilers will issue a warning if a variable is unused.
!    Sometimes there's a good reason to include a variable in a program,
!    but not to use it.  Calling this function with that variable as
!    the argument will shut the compiler up.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 April 2020
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) X, the variable to be "used".
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) x

  if ( x /= x ) then
    write ( *, '(a)' ) '  r8_fake_use(): variable is NAN.'
  end if

  return
end
subroutine r8mat_transpose_print ( m, n, a, title )

!*****************************************************************************80
!
!! r8mat_transpose_print() prints an R8MAT, transposed.
!
!  Discussion:
!
!    An R8MAT is an array of R8 values.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    14 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, N, the number of rows and columns.
!
!    Input, real ( kind = rk ) A(M,N), an M by N matrix to be printed.
!
!    Input, character ( len = * ) TITLE, a title.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m
  integer n

  real ( kind = rk ) a(m,n)
  character ( len = * ) title

  call r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title )

  return
end
subroutine r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

!*****************************************************************************80
!
!! r8mat_transpose_print_some() prints some of an R8MAT, transposed.
!
!  Discussion:
!
!    An R8MAT is an array of R8 values.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    14 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, N, the number of rows and columns.
!
!    Input, real ( kind = rk ) A(M,N), an M by N matrix to be printed.
!
!    Input, integer ILO, JLO, the first row and column to print.
!
!    Input, integer IHI, JHI, the last row and column to print.
!
!    Input, character ( len = * ) TITLE, a title.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer, parameter :: incx = 5
  integer m
  integer n

  real ( kind = rk ) a(m,n)
  character ( len = 14 ) ctemp(incx)
  integer i
  integer i2
  integer i2hi
  integer i2lo
  integer ihi
  integer ilo
  integer inc
  integer j
  integer j2hi
  integer j2lo
  integer jhi
  integer jlo
  character ( len = * ) title

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( title )

  do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx

    i2hi = i2lo + incx - 1
    i2hi = min ( i2hi, m )
    i2hi = min ( i2hi, ihi )

    inc = i2hi + 1 - i2lo

    write ( *, '(a)' ) ' '

    do i = i2lo, i2hi
      i2 = i + 1 - i2lo
      write ( ctemp(i2), '(i8,6x)' ) i
    end do

    write ( *, '(''  Row   '',5a14)' ) ctemp(1:inc)
    write ( *, '(a)' ) '  Col'
    write ( *, '(a)' ) ' '

    j2lo = max ( jlo, 1 )
    j2hi = min ( jhi, n )

    do j = j2lo, j2hi

      do i2 = 1, inc
        i = i2lo - 1 + i2
        write ( ctemp(i2), '(g14.6)' ) a(i,j)
      end do

      write ( *, '(i5,1x,5a14)' ) j, ( ctemp(i), i = 1, inc )

    end do

  end do

  return
end
subroutine reference_to_physical_t3 ( t, n, ref, phy )

!*****************************************************************************80
!
!! reference_to_physical_t3() maps T3 reference points to physical points.
!
!  Discussion:
!
!    Given the vertices of an order 3 physical triangle and a point 
!    (XSI,ETA) in the reference triangle, the routine computes the value 
!    of the corresponding image point (X,Y) in physical space.
!
!    This routine is also appropriate for an order 4 triangle,
!    as long as the fourth node is the centroid of the triangle.
!
!    This routine may also be appropriate for an order 6
!    triangle, if the mapping between reference and physical space
!    is linear.  This implies, in particular, that the sides of the
!    image triangle are straight and that the "midside" nodes in the
!    physical triangle are halfway along the sides of
!    the physical triangle.
!
!  Reference Element T3:
!
!    |
!    1  3
!    |  |\
!    |  | \
!    S  |  \
!    |  |   \
!    |  |    \
!    0  1-----2
!    |
!    +--0--R--1-->
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 June 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) T(2,3), the coordinates of the vertices.  
!    The vertices are assumed to be the images of (0,0), (1,0) and 
!    (0,1) respectively.
!
!    Input, integer N, the number of points to transform.
!
!    Input, real ( kind = rk ) REF(2,N), points in the reference element.
!
!    Output, real ( kind = rk ) PHY(2,N), corresponding points in the
!    physical element.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  integer i
  real ( kind = rk ) phy(2,n)
  real ( kind = rk ) ref(2,n)
  real ( kind = rk ) t(2,3)

  do i = 1, 2
    phy(i,1:n) = t(i,1) * ( 1.0D+00 - ref(1,1:n) - ref(2,1:n) ) &
               + t(i,2) *             ref(1,1:n)                &
               + t(i,3) *                          ref(2,1:n)
  end do

  return
end
subroutine triangle_area ( t, area )

!*****************************************************************************80
!
!! triangle_area() computes the area of a triangle.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    17 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) T(2,3), the triangle vertices.
!
!    Output, real ( kind = rk ) AREA, the absolute area of the triangle.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) area
  real ( kind = rk ) t(2,3)

  area = 0.5D+00 * abs ( &
      t(1,1) * ( t(2,2) - t(2,3) ) &
    + t(1,2) * ( t(2,3) - t(2,1) ) &
    + t(1,3) * ( t(2,1) - t(2,2) ) )

  return
end
subroutine triangle_integrand_01 ( p_num, p, f_num, fp )

!*****************************************************************************80
!
!! triangle_integrand_01() evaluates 1 integrand function.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    15 August 2009
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer P_NUM, the number of points.
!
!    Input, real ( kind = rk ) P(2,P_NUM), the evaluation points.
!
!    Input, integer F_NUM, the number of integrands.
!
!    Output, real ( kind = rk ) FP(F_NUM,P_NUM), the integrand values.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer f_num
  integer p_num

  real ( kind = rk ) fp(f_num,p_num)
  real ( kind = rk ) p(2,p_num)

  call r8_fake_use ( p(1,1) )

  fp(1,1:p_num) = 1.0D+00

  return
end
subroutine triangle_integrand_02 ( p_num, p, f_num, fp )

!*****************************************************************************80
!
!! triangle_integrand_02() evaluates 2 integrand functions.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    15 August 2009
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer P_NUM, the number of points.
!
!    Input, real ( kind = rk ) P(2,P_NUM), the evaluation points.
!
!    Input, integer F_NUM, the number of integrands.
!
!    Output, real ( kind = rk ) FP(F_NUM,P_NUM), the integrand values.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer f_num
  integer p_num

  real ( kind = rk ) fp(f_num,p_num)
  real ( kind = rk ) p(2,p_num)

  fp(1,1:p_num) = p(1,1:p_num)
  fp(2,1:p_num) = p(2,1:p_num)

  return
end
subroutine triangle_integrand_03 ( p_num, p, f_num, fp )

!*****************************************************************************80
!
!! triangle_integrand_03() evaluates 3 integrand functions.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    15 August 2009
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer P_NUM, the number of points.
!
!    Input, real ( kind = rk ) P(2,P_NUM), the evaluation points.
!
!    Input, integer F_NUM, the number of integrands.
!
!    Output, real ( kind = rk ) FP(F_NUM,P_NUM), the integrand values.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer f_num
  integer p_num

  real ( kind = rk ) fp(f_num,p_num)
  real ( kind = rk ) p(2,p_num)

  fp(1,1:p_num) = p(1,1:p_num) * p(1,1:p_num)
  fp(2,1:p_num) = p(1,1:p_num) * p(2,1:p_num)
  fp(3,1:p_num) = p(2,1:p_num) * p(2,1:p_num)

  return
end
subroutine triangle_integrand_04 ( p_num, p, f_num, fp )

!*****************************************************************************80
!
!! triangle_integrand_04() evaluates 4 integrand functions.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    15 August 2009
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer P_NUM, the number of points.
!
!    Input, real ( kind = rk ) P(2,P_NUM), the evaluation points.
!
!    Input, integer F_NUM, the number of integrands.
!
!    Output, real ( kind = rk ) FP(F_NUM,P_NUM), the integrand values.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer f_num
  integer p_num

  real ( kind = rk ) fp(f_num,p_num)
  real ( kind = rk ) p(2,p_num)

  fp(1,1:p_num) = p(1,1:p_num) * p(1,1:p_num) * p(1,1:p_num)
  fp(2,1:p_num) = p(1,1:p_num) * p(1,1:p_num) * p(2,1:p_num)
  fp(3,1:p_num) = p(1,1:p_num) * p(2,1:p_num) * p(2,1:p_num)
  fp(4,1:p_num) = p(2,1:p_num) * p(2,1:p_num) * p(2,1:p_num)

  return
end
subroutine triangle_integrand_05 ( p_num, p, f_num, fp )

!*****************************************************************************80
!
!! triangle_integrand_05() evaluates 5 integrand functions.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    15 August 2009
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer P_NUM, the number of points.
!
!    Input, real ( kind = rk ) P(2,P_NUM), the evaluation points.
!
!    Input, integer F_NUM, the number of integrands.
!
!    Output, real ( kind = rk ) FP(F_NUM,P_NUM), the integrand values.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer f_num
  integer p_num

  real ( kind = rk ) fp(f_num,p_num)
  real ( kind = rk ) p(2,p_num)

  fp(1,1:p_num) = p(1,1:p_num)**4
  fp(2,1:p_num) = p(1,1:p_num)**3 * p(2,1:p_num)
  fp(3,1:p_num) = p(1,1:p_num)**2 * p(2,1:p_num)**2
  fp(4,1:p_num) = p(1,1:p_num)    * p(2,1:p_num)**3
  fp(5,1:p_num) =                   p(2,1:p_num)**4
  
  return
end
subroutine triangle_monte_carlo ( t, p_num, f_num, triangle_unit_sample, &
  triangle_integrand, result )

!*****************************************************************************80
!
!! triangle_monte_carlo() applies the Monte Carlo rule to integrate a function.
!
!  Discussion:
!
!    The function f(x,y) is to be integrated over a triangle T.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    15 August 2009
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) T(2,3), the triangle vertices.
!
!    Input, integer P_NUM, the number of sample points.
!
!    Input, integer F_NUM, the number of functions to integrate.
!
!    Input, external TRIANGLE_UNIT_SAMPLE, the sampling routine.
!
!    Input, external TRIANGLE_INTEGRAND, the integrand routine.
!
!    Output, real ( kind = rk ) RESULT(F_NUM), the approximate integrals.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer f_num
  integer p_num

  real ( kind = rk ) area
  real ( kind = rk ) fp(f_num,p_num)
  integer i
  real ( kind = rk ) p(2,p_num)
  real ( kind = rk ) p2(2,p_num)
  real ( kind = rk ) result(f_num)
  real ( kind = rk ) t(2,3)
  external triangle_sample
  external triangle_integrand

  call triangle_area ( t, area )

  call triangle_unit_sample ( p_num, p )

  call reference_to_physical_t3 ( t, p_num, p, p2 )

  call triangle_integrand ( p_num, p2, f_num, fp )

  do i = 1, f_num
    result(i) = area * sum ( fp(i,1:p_num) ) / real ( p_num, kind = rk )
  end do

  return
end
subroutine triangle_unit_sample_01 ( p_num, p )

!*****************************************************************************80
!
!! triangle_unit_sample_01() selects points from the unit triangle.
!
!  Discussion:
!
!    The unit triangle has vertices (1,0), (0,1), (0,0).
!
!    Any point in the unit triangle CAN be chosen by this algorithm.
!
!    However, the points that are chosen tend to be clustered near
!    the centroid.
!
!    This routine is supplied as an example of "bad" sampling.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    14 August 2009
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer P_NUM, the number of points.
!
!    Output, real ( kind = rk ) P(2,P_NUM), the points.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer p_num

  real ( kind = rk ) e(3)
  real ( kind = rk ) e_sum
  integer j
  real ( kind = rk ) p(2,p_num)

  do j = 1, p_num

    call random_number ( harvest = e(1:3) )

    e_sum = sum ( e(1:3) )

    e(1:3) = e(1:3) / e_sum
!
!  We take E(1:3) as the barycentric coordinates of the point.
!
    p(1:2,j) = e(1:2)

  end do

  return
end
subroutine triangle_unit_sample_02 ( p_num, p )

!*****************************************************************************80
!
!! triangle_unit_sample_02() selects points from the unit triangle.
!
!  Discussion:
!
!    The unit triangle has vertices (1,0), (0,1), (0,0).
!
!    The sampling is uniform.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    19 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer P_NUM, the number of points.
!
!    Output, real ( kind = rk ) P(2,P_NUM), the points.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer p_num

  integer j
  real ( kind = rk ) r(2)
  real ( kind = rk ) p(2,p_num)
!
!  Generate the points using barycentric coordinates.
!
  do j = 1, p_num

    call random_number ( harvest = r(1:2) )

    if ( 1.0D+00 < sum ( r(1:2) ) ) then
      r(1:2) = 1.0D+00 - r(1:2)
    end if

    p(1:2,j) = r(1:2)

  end do

  return
end
subroutine triangle_unit_sample_03 ( p_num, p )

!*****************************************************************************80
!
!! triangle_unit_sample_03() selects points from the unit triangle.
!
!  Discussion:
!
!    The unit triangle has vertices (1,0), (0,1), (0,0).
!
!    This routine uses Turk's rule #1.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    19 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Greg Turk,
!    Generating Random Points in a Triangle,
!    in Graphics Gems,
!    edited by Andrew Glassner,
!    AP Professional, 1990, pages 24-28.
!
!  Parameters:
!
!    Input, integer P_NUM, the number of points.
!
!    Output, real ( kind = rk ) P(2,P_NUM), the points.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer p_num

  real ( kind = rk ) a
  real ( kind = rk ) b
  real ( kind = rk ) c
  integer j
  real ( kind = rk ) p(2,p_num)
  real ( kind = rk ) r(2)
!
!  Generate the points using Turk's rule 1.
!
  do j = 1, p_num

    call random_number ( harvest = r(1:2) )

    a = 1.0D+00            - sqrt ( r(2) )
    b = ( 1.0D+00 - r(1) ) * sqrt ( r(2) )
    c =             r(1)   * sqrt ( r(2) )

    p(1,j) = a
    p(2,j) = b

  end do

  return
end
subroutine triangle_unit_sample_04 ( p_num, p )

!*****************************************************************************80
!
!! triangle_unit_sample_04() selects points from the unit triangle.
!
!  Discussion:
!
!    The unit triangle has vertices (1,0), (0,1), (0,0).
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    08 July 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Reuven Rubinstein,
!    Monte Carlo Optimization, Simulation, and Sensitivity 
!    of Queueing Networks,
!    Krieger, 1992,
!    ISBN: 0894647644,
!    LC: QA298.R79.
!
!  Parameters:
!
!    Input, integer P_NUM, the number of points.
!
!    Output, real ( kind = rk ) P(2,P_NUM), the points.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer p_num

  real ( kind = rk ) e(3)
  integer j
  real ( kind = rk ) p(2,p_num)
!
!  The construction begins by sampling DIM_NUM+1 points from the
!  exponential distribution with parameter 1.
!
  do j = 1, p_num

    call random_number ( harvest = e(1:3) )

    e(1:3) = - log ( e(1:3) )

    p(1:2,j) = e(1:2) / sum ( e(1:3) )

  end do

  return
end