subroutine i4_to_pascal ( k, i, j )

!*****************************************************************************80
!
!! I4_TO_PASCAL converts a linear index to Pascal triangle coordinates.
!
!  Discussion:
!
!    We describe the grid points in Pascal's triangle in two ways:
!
!    As a linear index K:
!
!                     1
!                   2   3
!                 4   5   6
!               7   8   9   10
!
!    As elements (I,J) of Pascal's triangle:
!
!                     0,0
!                  1,0   0,1
!               2,0   1,1    0,2
!            3,0   2,1   1,2    0,3
!
!  Example:
!
!     K  I  J
!
!     1  0  0
!     2  1  0
!     3  0  1
!     4  2  0
!     5  1  1
!     6  0  2
!     7  3  0
!     8  2  1
!     9  1  2
!    10  0  3
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    14 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer K, the linear index of the (I,J) element.
!    1 <= K.
!
!    Output, integer I, J, the Pascal indices.
!
  implicit none

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

  integer d
  integer i
  integer j
  integer k

  if ( k <= 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'I4_TO_PASCAL - Fatal error!'
    write ( *, '(a)' ) '  K must be positive.'
    stop 1
  end if

  call i4_to_pascal_degree ( k, d )

  j = k - ( d * ( d + 1 ) ) / 2 - 1
  i = d - j

  return
end
subroutine i4_to_pascal_degree ( k, d )

!*****************************************************************************80
!
!! I4_TO_PASCAL_DEGREE converts a linear index to a Pascal triangle degree.
!
!  Discussion:
!
!    We describe the grid points in Pascal's triangle in two ways:
!
!    As a linear index K:
!
!                     1
!                   2   3
!                 4   5   6
!               7   8   9   10
!
!    As elements (I,J) of Pascal's triangle:
!
!                     0,0
!                  1,0   0,1
!               2,0   1,1    0,2
!            3,0   2,1   1,2    0,3
!
!    The quantity D represents the "degree" of the corresponding monomial,
!    that is, D = I + J.
!
!    We can compute D directly from K using the quadratic formula.
!
!  Example:
!
!     K  I  J  D
!
!     1  0  0  0
!
!     2  1  0  1
!     3  0  1  1
!
!     4  2  0  2
!     5  1  1  2
!     6  0  2  2
!
!     7  3  0  3
!     8  2  1  3
!     9  1  2  3
!    10  0  3  3
!
!    11  4  0  4
!    12  3  1  4
!    13  2  2  4
!    14  1  3  4
!    15  0  4  4
!
!    16  5  0  5
!    17  4  1  5
!    18  3  2  5
!    19  2  3  5
!    20  1  4  5
!    21  0  5  5
!
!    22  6  0  6
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    14 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer K, the linear index of the (I,J) element.
!    1 <= K.
!
!    Output, integer D, the degree (sum) of the corresponding 
!    Pascal indices.
!
  implicit none

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

  real ( kind = rk ) arg
  integer d
  integer k

  if ( k <= 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'I4_TO_PASCAL_DEGREE - Fatal error!'
    write ( *, '(a)' ) '  K must be positive.'
    stop 1
  end if

  arg = real ( 1 + 8 * ( k - 1 ), kind = rk )

  d = int ( 0.5D+00 * ( -1.0D+00 + sqrt ( arg ) ) )

  return
end
subroutine pascal_to_i4 ( i, j, k )

!*****************************************************************************80
!
!! PASCAL_TO_I4 converts Pacal triangle coordinates to a linear index.
!
!  Discussion:
!
!    We describe the grid points in a Pascal triangle in two ways:
!
!    As a linear index K:
!
!                     1
!                   2   3
!                 4   5   6
!               7   8   9   10
!
!    As elements (I,J) of Pascal's triangle:
!
!                     0,0
!                  1,0   0,1
!               2,0   1,1    0,2
!            3,0   2,1   1,2    0,3
!
!  Example:
!
!     K  I  J
!
!     1  0  0
!     2  1  0
!     3  0  1
!     4  2  0
!     5  1  1
!     6  0  2
!     7  3  0
!     8  2  1
!     9  1  2
!    10  0  3
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    14 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, J, the row and column indices.  I and J 
!    must be nonnegative.
!
!    Output, integer K, the linear index of the (I,J) element.
!
  implicit none

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

  integer d
  integer i
  integer j
  integer k

  if ( i < 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PASCAL_TO_I4 - Fatal error!'
    write ( *, '(a)' ) '  I < 0.'
    write ( *, '(a,i8)' ) '  I = ', i
    stop 1
  else if ( j < 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PASCAL_TO_I4 - Fatal error!'
    write ( *, '(a)' ) '  J < 0.'
    write ( *, '(a,i8)' ) '  J = ', j
    stop 1
  end if

  d = i + j

  k = ( d * ( d + 1 ) ) / 2 + j + 1

  return
end
subroutine poly_power ( d1, p1, n, d2, p2 )

!*****************************************************************************80
!
!! POLY_POWER computes a power of a polynomial.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer D1, the degree of the polynomial.
!
!    Input, real ( kind = rk ) P1(M1), the polynomial coefficients.
!    M1 = ((D1+1)*(D1+2))/2.
!
!    Input, integer N, the nonnegative integer power.
!
!    Input, integer D2, the degree of the power polynomial.
!    D2 = N * D1.
!
!    Output, real ( kind = rk ) P2(M2), the polynomial power.
!    M2 = ((D2+1)*(D2+2))/2.
!
  implicit none

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

  integer d1
  integer d2
  integer n

  integer d3
  integer d4
  integer i
  integer m2
  integer m3
  integer m4
  real ( kind = rk ) p1(( ( d1 + 1 ) * ( d1 + 2 ) ) / 2)
  real ( kind = rk ) p2(( ( d2 + 1 ) * ( d2 + 2 ) ) / 2)
  real ( kind = rk ), allocatable :: p3(:)
  real ( kind = rk ), allocatable :: p4(:)

  m2 = ( ( d2 + 1 ) * ( d2 + 2 ) ) / 2
!
!  Create P3, a polynomial representation of 1, that is
!  big enough to store the final result.
!
  d3 = n * d1
  m3 = ( ( d3 + 1 ) * ( d3 + 2 ) ) / 2
  allocate ( p3(1:m3) )
  p3(1) = 1.0D+00
  p3(2:m3) = 0.0D+00
!
!  Create P4, big enough to hold the result.
!
  d4 = n * d1
  m4 = ( ( d4 + 1 ) * ( d4 + 2 ) ) / 2
  allocate ( p4(1:m4) )
  p4(1:m4) = 0.0D+00
!
!  Now set D3 to 0, to indicate that P3 currently contains only
!  a constant term.
!
  d3 = 0
!
!  Iterate N times:
!    P <= P1 * P
!
  do i = 1, n
    d4 = d1 + d3
    call poly_product ( d1, p1, d3, p3, d4, p4 )
    d3 = d4
    p3(1:m3) = p4(1:m4)
  end do
!
!  Copy the result to the user.
!
  p2(1:m2) = p3(1:m3)
!
!  Free memory.
!
  deallocate ( p3 )
  deallocate ( p4 )

  return
end
subroutine poly_power_linear ( d1, p1, n, d2, p2 )

!*****************************************************************************80
!
!! POLY_POWER_LINEAR computes the polynomial ( a + b*x + c*y ) ^ n.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer D1, the degree of the linear polynomial,
!    which should be 1 (or possibly 0).
!
!    Input, real ( kind = rk ) P1(M1), the coefficients of the linear polynomial.
!    M1 = ( (D1+1)*(D1+2) ) / 2, which should be 3.
!
!    Input, integer N, the power to which the polynomial is to be 
!    raised.  0 <= N.
!
!    Input, integer D2, the degree of the power polyynomial.
!    D2 = N * D1 = N.
!
!    Output, real ( kind = rk ) P2(M2), the coefficients of the power polynomial.
!    M2 = ( (D2+1)*(D2+2) ) / 2, which should be ((N+1)*(N+2))/2.
!
  implicit none

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

  integer d1
  integer d2
  integer n

  integer i
  integer j
  integer k
  integer l
  real ( kind = rk ) p1(( ( d1 + 1 ) * ( d1 + 2 ) ) / 2 )
  real ( kind = rk ) p2(( ( d2 + 1 ) * ( d2 + 2 ) ) / 2 )
  integer trinomial

  if ( d1 < 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'POLY_POWER_LINEAR - Fatal error!'
    write ( *, '(a)' ) '  D1 < 0.'
    stop 1
  end if

  if ( n < 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'POLY_POWER_LINEAR - Fatal error!'
    write ( *, '(a)' ) '  N < 0.'
    stop 1
  end if

  if ( d1 == 0 ) then
    p2(1) = p1(1) ** n
    return
  end if

  if ( n == 0 ) then
    p2(1) = 1.0D+00
    return
  end if
!
!  Use the Trinomial formula.
!
  do i = 0, n
    do j = 0, n - i
      do k = 0, n - i - j
!
!  We store X^J Y^K in location L.
!
        call pascal_to_i4 ( j, k, l )
        p2(l) = real ( trinomial ( i, j, k ), kind = rk ) &
          * p1(1) ** i * p1(2) ** j * p1(3) ** k
      end do
    end do
  end do

  return
end
subroutine poly_print ( d, p, title )

!*****************************************************************************80
!
!! POLY_PRINT prints an XY polynomial.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer D, the degree of the polynomial.
!
!    Output, real ( kind = rk ) P(M), the coefficients of all monomials of 
!    degree 0 through D.  P must contain ((D+1)*(D+2))/2 entries.
!
!    Input, character ( len = * ) TITLE, a title string.
!
  implicit none

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

  integer d

  integer i
  integer j
  integer k
  integer m
  real ( kind = rk ) p(( ( d + 1 ) * ( d + 2 ) ) / 2)
  character ( len = * ) title

  m = ( ( d + 1 ) * ( d + 2 ) ) / 2

  if ( all ( p(1:m) == 0.0D+00 ) ) then

    write ( *, '(a)' ) trim ( title ) // ' = 0'

  else

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

    do k = 1, m

      call i4_to_pascal ( k, i, j )

      if ( p(k) /= 0.0D+00 ) then

        if ( p(k) < 0.0D+00 ) then
          write ( *, '(2x,a,g14.6)', advance = 'no' ) '-', abs ( p(k) )
        else
          write ( *, '(2x,a,g14.6)', advance = 'no' ) '+', p(k)
        end if

        if ( i + j /= 0 ) then
          write ( *, '(a)', advance = 'no' ) ' '
        end if

        if ( i == 0 ) then
        else if ( i == 1 ) then
          write ( *, '(a)', advance = 'no' ) 'x'
        else
          write ( *, '(a,i2)', advance = 'no' ) 'x^', i
        end if

        if ( j == 0 ) then
        else if ( j == 1 ) then
          write ( *, '(a)', advance = 'no' ) 'y'
        else
          write ( *, '(a,i2)', advance = 'no' ) 'y^', j
        end if

        write ( *, '(a)' )

      end if

    end do

  end if

  return
end
subroutine poly_product ( d1, p1, d2, p2, d3, p3 )

!*****************************************************************************80
!
!! POLY_PRODUCT computes P3(x,y) = P1(x,y) * P2(x,y) for polynomials.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer D1, the degree of factor 1.
!
!    Input, real ( kind = rk ) P1(M1), the factor 1 coefficients.
!    M1 = ((D1+1)*(D1+2))/2.
!
!    Input, integer D2, the degree of factor 2.
!
!    Input, real ( kind = rk ) P2(M2), the factor2 coefficients.
!    M2 = ((D2+1)*(D2+2))/2.
!
!    Input, integer D3, the degree of the result.
!    D3 = D1 + D2.
!
!    Output, real ( kind = rk ) P3(M3), the result coefficients.
!    M3 = ((D3+1)*(D3+2))/2.
!
  implicit none

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

  integer d1
  integer d2
  integer d3

  integer i1
  integer i2
  integer i3
  integer j1
  integer j2
  integer j3
  integer k1
  integer k2
  integer k3
  integer m1
  integer m2
  integer m3
  real ( kind = rk ) p1(( ( d1 + 1 ) * ( d1 + 2 ) ) / 2)
  real ( kind = rk ) p2(( ( d2 + 1 ) * ( d2 + 2 ) ) / 2)
  real ( kind = rk ) p3(( ( d3 + 1 ) * ( d3 + 2 ) ) / 2)

  m1 = ( ( d1 + 1 ) * ( d1 + 2 ) ) / 2
  m2 = ( ( d2 + 1 ) * ( d2 + 2 ) ) / 2
  m3 = ( ( d3 + 1 ) * ( d3 + 2 ) ) / 2
!
!  Consider each entry in P1:
!    P1(K1) * X^I1 * Y^J1
!  and multiply it by each entry in P2:
!    P2(K2) * X^I2 * Y^J2
!  getting 
!    P3(K3) = P3(K3) + P1(K1) * P2(X2) * X^(I1+I2) * Y(J1+J2)
!
  p3(1:m3) = 0.0D+00

  do k1 = 1, m1
    call i4_to_pascal ( k1, i1, j1 )
    do k2 = 1, m2
      call i4_to_pascal ( k2, i2, j2 )
      i3 = i1 + i2
      j3 = j1 + j2
      call pascal_to_i4 ( i3, j3, k3 )
      p3(k3) = p3(k3) + p1(k1) * p2(k2)
    end do
  end do

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

!*****************************************************************************80
!
!! R8MAT_PRINT prints an R8MAT.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    12 September 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the number of rows in A.
!
!    Input, integer N, the number of columns in A.
!
!    Input, real ( kind = rk ) A(M,N), the matrix.
!
!    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_print_some ( m, n, a, 1, 1, m, n, title )

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

!*****************************************************************************80
!
!! R8MAT_PRINT_SOME prints some of an R8MAT.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    10 September 2009
!
!  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 i2hi
  integer i2lo
  integer ihi
  integer ilo
  integer inc
  integer j
  integer j2
  integer j2hi
  integer j2lo
  integer jhi
  integer jlo
  character ( len = * ) title

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

  if ( m <= 0 .or. n <= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  (None)'
    return
  end if

  do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx

    j2hi = j2lo + incx - 1
    j2hi = min ( j2hi, n )
    j2hi = min ( j2hi, jhi )

    inc = j2hi + 1 - j2lo

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

    do j = j2lo, j2hi
      j2 = j + 1 - j2lo
      write ( ctemp(j2), '(i8,6x)' ) j
    end do

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

    i2lo = max ( ilo, 1 )
    i2hi = min ( ihi, m )

    do i = i2lo, i2hi

      do j2 = 1, inc

        j = j2lo - 1 + j2

        if ( a(i,j) == real ( int ( a(i,j) ), kind = rk ) ) then
          write ( ctemp(j2), '(f8.0,6x)' ) a(i,j)
        else
          write ( ctemp(j2), '(g14.6)' ) a(i,j)
        end if

      end do

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

    end do

  end do

  return
end
subroutine rs_to_xy_map ( t, a, b, c, d, e, f )

!*****************************************************************************80
!
!! RS_TO_XY_MAP returns the linear map from reference to physical triangle.
!
!  Discussion:
!
!    This function returns the coefficients of the linear map that sends
!    the vertices of the reference triangle, (0,0), (1,0) and (0,1), to
!    the vertices of a physical triangle T, of the form:
!
!      X = A + B * R + C * S;
!      Y = D + E * R + F * S.
!
!  Reference Element:
!
!    |
!    1  3
!    |  |\
!    |  | \
!    S  |  \
!    |  |   \
!    |  |    \
!    0  1-----2
!    |
!    +--0--R--1-->
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 April 2015
!
!  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.
!
!    Output, real ( kind = rk ) A, B, C, D, E, F, the mapping coefficients.
!
  implicit none

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

  real ( kind = rk ) a
  real ( kind = rk ) b
  real ( kind = rk ) c
  real ( kind = rk ) d
  real ( kind = rk ) e
  real ( kind = rk ) f
  real ( kind = rk ) t(2,3)

  a = t(1,1)
  b = t(1,2) - t(1,1)
  c = t(1,3) - t(1,1)

  d = t(2,1)
  e = t(2,2) - t(2,1)
  f = t(2,3) - t(2,1)

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

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

  character ( len = 8 ) ampm
  integer d
  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
  integer values(8)
  integer y

  call date_and_time ( values = 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 ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end
function triangle01_monomial_integral ( i, j )

!*****************************************************************************80
!
!! TRIANGLE01_MONOMIAL_INTEGRAL: monomial integrals in the unit triangle.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    19 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, J, the exponents.  
!    Each exponent must be nonnegative.
!
!    Output, real ( kind = rk ) TRIANGLE01_MONOMIAL_INTEGRAL, the integral.
!
  implicit none

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

  integer i
  integer j
  integer k
  integer l
  real ( kind = rk ) q
  real ( kind = rk ) triangle01_monomial_integral

  k = 0
  q = 1.0D+00

  do l = 1, i
    k = k + 1
    q = q * real ( l, kind = rk ) / real ( k, kind = rk )
  end do

  do l = 1, j
    k = k + 1
    q = q * real ( l, kind = rk ) / real ( k, kind = rk )
  end do

  do l = 1, 2
    k = k + 1
    q = q / real ( k, kind = rk )
  end do

  triangle01_monomial_integral = q

  return
end
function triangle01_poly_integral ( d, p )

!*****************************************************************************80
!
!! TRIANGLE01_POLY_INTEGRAL: polynomial integral over the unit triangle.
!
!  Discussion:
!
!    The unit triangle is T = ( (0,0), (1,0), (0,1) ).
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    19 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer D, the degree of the polynomial.
!
!    Input, real ( kind = rk ) P(M), the polynomial coefficients.
!    M = ((D+1)*(D+2))/2.
!
!    Output, real ( kind = rk ) TRIANGLE01_POLY_INTEGRAL, the integral.
!
  implicit none

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

  integer d

  integer i
  integer j
  integer k
  integer m
  real ( kind = rk ) p(( ( d + 1 ) * ( d + 2 ) ) / 2)
  real ( kind = rk ) q
  real ( kind = rk ) triangle01_poly_integral
  real ( kind = rk ) triangle01_monomial_integral

  m = ( ( d + 1 ) * ( d + 2 ) ) / 2

  q = 0.0D+00
  do k = 1, m
    call i4_to_pascal ( k, i, j )
    q = q + p(k) * triangle01_monomial_integral ( i, j )
  end do

  triangle01_poly_integral = q

  return
end
function triangle_area ( t )

!*****************************************************************************80
!
!! TRIANGLE_AREA returns the area of a triangle.
!
!  Discussion:
!
!    If the vertices are given in counter clockwise order, the area
!    will be positive.
!
!  Licensing:
!
!    This code is distributed under the GNU GPL license.
!
!  Modified:
!
!    18 April 2015
!
!  Author:
!
!    John Burkardt.
!
!  Parameters:
!
!    Input, real ( kind = rk ) T(2,3), the vertices of the triangle.
!
!    Output, real ( kind = rk ) TRIANGLE_AREA, the area of the triangle.
!
  implicit none

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

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

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

  return
end
function triangle_monomial_integral ( i, j, t )

!*****************************************************************************80
!
!! TRIANGLE_MONOMIAL_INTEGRAL integrates a monomial over an arbitrary triangle.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    19 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, J, the exponents of X and Y in the monomial.
!    0 <= I, J.
!
!    Input, real ( kind = rk ) T(2,3), the vertices of the triangle.
!
!    Output, real ( kind = rk ) TRIANGLE_MONOMIAL_INTEGRAL, the integral 
!    of X^I * Y^J over triangle T.
!
  implicit none

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

  integer, parameter :: d1 = 1
  integer, parameter :: d2 = 1

  integer, parameter :: m1 = ( ( d1 + 1 ) * ( d1 + 2 ) ) / 2
  integer, parameter :: m2 = ( ( d2 + 1 ) * ( d2 + 2 ) ) / 2

  real ( kind = rk ) a
  real ( kind = rk ) b
  real ( kind = rk ) c
  real ( kind = rk ) d
  integer d3
  integer d4
  integer d5
  real ( kind = rk ) e
  real ( kind = rk ) f
  integer i
  integer j
  integer m3
  integer m4
  integer m5
  real ( kind = rk ) p1(m1)
  real ( kind = rk ) p2(m2)
  real ( kind = rk ), allocatable :: p3(:)
  real ( kind = rk ), allocatable :: p4(:)
  real ( kind = rk ), allocatable :: p5(:)
  real ( kind = rk ) q
  real ( kind = rk ) t(2,3)
  real ( kind = rk ) triangle_area
  real ( kind = rk ) triangle_monomial_integral
  real ( kind = rk ) triangle01_poly_integral
!
!  Get map coefficients from reference RS triangle to general XY triangle.
!    R = a+b*X+c*Y
!    S = d+e*X+f*Y
!
  call rs_to_xy_map ( t, a, b, c, d, e, f )
!
!  Set
!    P1(R,S) = a+b*R+c*S
!    P2(R,S) = d+e*R+f*S
!
  p1(1) = a
  p1(2) = b
  p1(3) = c

  p2(1) = d
  p2(2) = e
  p2(3) = f
!
!  Exponentiate:
!    P3(R,S) = P1(R,S)^i
!    P4(R,S) = P2(R,S)^j
!
  d3 = i * d1
  m3 = ( ( d3 + 1 ) * ( d3 + 2 ) ) / 2
  allocate ( p3(1:m3) )
  call poly_power_linear ( d1, p1, i, d3, p3 )

  d4 = j * d2
  m4 = ( ( d4 + 1 ) * ( d4 + 2 ) ) / 2
  allocate ( p4(1:m4) )
  call poly_power_linear ( d2, p2, j, d4, p4 )
!
!  Compute the product 
!    P5(R,S) = P3(R,S) * P4(R,S)
!
  d5 = d3 + d4
  m5 = ( ( d5 + 1 ) * ( d5 + 2 ) ) / 2
  allocate ( p5(1:m5) )
  call poly_product ( d3, p3, d4, p4, d5, p5 )
!
!  Compute the integral of P5(R,S) over the reference triangle.
!
  q = triangle01_poly_integral ( d5, p5 )
!
!  Multiply by the area of the physical triangle T(X,Y) divided by
!  the area of the reference triangle.
!
  q = q * triangle_area ( t ) / 0.5D+00

  triangle_monomial_integral = q

  deallocate ( p3 )
  deallocate ( p4 )
  deallocate ( p5 )

  return
end
function triangle_poly_integral ( d, p, t )

!*****************************************************************************80
!
!! TRIANGLE_POLY_INTEGRAL: polynomial integral over a triangle.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    19 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer D, the degree of the polynomial.
!
!    Input, real ( kind = rk ) P(M), the polynomial coefficients.
!    M = ((D+1)*(D+2))/2.
!
!    Input, real ( kind = rk ) T(2,3), the vertices of the triangle.
!
!    Output, real ( kind = rk ) TRIANGLE_POLY_INTEGRAL, the integral.
!
  implicit none

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

  integer d

  integer i
  integer j
  integer k
  integer m
  real ( kind = rk ) p(( ( d + 1 ) * ( d + 2 ) ) / 2)
  real ( kind = rk ) q
  real ( kind = rk ) t(2,3)
  real ( kind = rk ) triangle_monomial_integral
  real ( kind = rk ) triangle_poly_integral

  m = ( ( d + 1 ) * ( d + 2 ) ) / 2

  q = 0.0D+00
  do k = 1, m
    call i4_to_pascal ( k, i, j )
    q = q + p(k) * triangle_monomial_integral ( i, j, t )
  end do

  triangle_poly_integral = q

  return
end
function triangle_xy_integral ( x1, y1, x2, y2, x3, y3 )

!*****************************************************************************80
!
!! TRIANGLE_XY_INTEGRAL computes the integral of XY over a triangle.
!
!  Discussion:
!
!    This function was written as a special test case for the general
!    problem of integrating a monomial x^alpha * y^beta over a general 
!    triangle.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) X1, Y1, X2, Y2, X3, Y3, the coordinates of the
!    triangle vertices.
!
!    Output, real ( kind = rk ) TRIANGLE_XY_INTEGRAL, the integral of X*Y 
!    over the triangle.
!
  implicit none

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

  real ( kind = rk ) det
  real ( kind = rk ) p00
  real ( kind = rk ) p01
  real ( kind = rk ) p02
  real ( kind = rk ) p10
  real ( kind = rk ) p11
  real ( kind = rk ) p20
  real ( kind = rk ) q
  real ( kind = rk ) triangle01_monomial_integral
  real ( kind = rk ) triangle_xy_integral
  real ( kind = rk ) x1
  real ( kind = rk ) x2
  real ( kind = rk ) x3
  real ( kind = rk ) y1
  real ( kind = rk ) y2
  real ( kind = rk ) y3
!
!  x = x1 * ( 1 - xi - eta )
!    + x2 *       xi
!    + x3 *            eta
!
!  y = y1 * ( 1 - xi - eta )
!    + y2 *       xi
!    + y3 *            eta
!
!  Rewrite as linear polynomials in (xi,eta):
!
!  x = x1 + ( x2 - x1 ) * xi + ( x3 - x1 ) * eta
!  y = y1 + ( y2 - y1 ) * xi + ( y3 - y1 ) * eta
!
!  Jacobian:
!
!    J = [ ( x2 - x1 )  ( x3 - x1 ) ]
!        [ ( y2 - y1 )  ( y3 - y1 ) ]
!
!    det J = ( x2 - x1 ) * ( y3 - y1 ) - ( y2 - y1 ) * ( x3 - x1 )
!
!  Integrand
!
!    x * y = ( x1 + ( x2 - x1 ) * xi + ( x3 - x1 ) * eta )
!          * ( y1 + ( y2 - y1 ) * xi + ( y3 - y1 ) * eta )
!
!  Rewrite as linear combination of monomials:
!
!    x * y = 1      * x1 * y1
!          + eta    * ( x1 * ( y3 - y1 ) + ( x3 - x1 ) * y1 )
!          + xi     * ( x1 * ( y2 - y1 ) + ( x2 - x1 ) * y1 )
!          + eta^2  * ( x3 - x1 ) * ( y3 - y1 )
!          + xi*eta * ( ( x2 - x1 ) * ( y3 - y1 ) + ( x3 - x1 ) * ( y2 - y1 ) )
!          + xi^2   * ( x2 - x1 ) * ( y2 - y1 )
!
  det = ( x2 - x1 ) * ( y3 - y1 ) - ( y2 - y1 ) * ( x3 - x1 )

  p00 = x1 * y1

  p01 = x1 * ( y3 - y1 ) + ( x3 - x1 ) * y1
  p10 = x1 * ( y2 - y1 ) + ( x2 - x1 ) * y1

  p02 = ( x3 - x1 ) * ( y3 - y1 )
  p11 = ( x2 - x1 ) * ( y3 - y1 ) + ( x3 - x1 ) * ( y2 - y1 )
  p20 = ( x2 - x1 ) * ( y2 - y1 )

  q = 0.0D+00
  q = q + p00 * triangle01_monomial_integral ( 0, 0 )
  q = q + p10 * triangle01_monomial_integral ( 1, 0 )
  q = q + p01 * triangle01_monomial_integral ( 0, 1 )
  q = q + p20 * triangle01_monomial_integral ( 2, 0 )
  q = q + p11 * triangle01_monomial_integral ( 1, 1 )
  q = q + p02 * triangle01_monomial_integral ( 0, 2 )

  q = q * det

  triangle_xy_integral = q

  return
end
function trinomial ( i, j, k )

!*****************************************************************************80
!
!! TRINOMIAL computes a trinomial coefficient.
!
!  Discussion:
!
!    The trinomial coefficient is a generalization of the binomial
!    coefficient.  It may be interpreted as the number of combinations of
!    N objects, where I objects are of type 1, J of type 2, and K of type 3.
!    and N = I + J + K.
!
!    T(I,J,K) = N! / ( I! J! K! )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    11 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, J, K, the factors.
!    All should be nonnegative.
!
!    Output, integer TRINOMIAL, the trinomial coefficient.
!
  implicit none

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

  integer i
  integer j
  integer k
  integer l
  integer t
  integer trinomial
  integer value
!
!  Each factor must be nonnegative.
!
  if ( i < 0 .or. j < 0 .or. k < 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TRINOMIAL - Fatal error!'
    write ( *, '(a)' ) '  Negative factor encountered.'
    stop 1
  end if

  value = 1

  t = 1

  do l = 1, i
!   value = value * t / l
    t = t + 1
  end do

  do l = 1, j
    value = value * t / l
    t = t + 1
  end do

  do l = 1, k
    value = value * t / l
    t = t + 1
  end do
  
  trinomial = value

  return
end
subroutine xy_to_rs_map ( t, a, b, c, d, e, f )

!*****************************************************************************80
!
!! XY_TO_RS_MAP returns the linear map from physical to reference triangle.
!
!  Discussion:
!
!    Given the vertices T of an arbitrary triangle in the (X,Y) coordinate
!    system, this function returns the coefficients of the linear map
!    that sends the vertices of T to (0,0), (1,0) and (0,1) respectively
!    in the reference triangle with coordinates (R,S):
!
!      R = A + B * X + C * Y;
!      S = D + E * X + F * Y.
!
!  Reference Element T3:
!
!    |
!    1  3
!    |  |\
!    |  | \
!    S  |  \
!    |  |   \
!    |  |    \
!    0  1-----2
!    |
!    +--0--R--1-->
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = rk ) T(2,3), the X and Y coordinates
!    of the vertices.  The vertices are assumed to be the images of
!    (0,0), (1,0) and (0,1) respectively.
!
!    Output, real ( kind = rk ) A, B, C, D, E, F, the mapping coefficients.
!
  implicit none

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

  real ( kind = rk ) a
  real ( kind = rk ) b
  real ( kind = rk ) c
  real ( kind = rk ) d
  real ( kind = rk ) e
  real ( kind = rk ) f
  real ( kind = rk ) g
  real ( kind = rk ) t(2,3)

  g =   ( ( t(2,3) - t(2,1) ) * ( t(1,2) - t(1,1) )   &
        - ( t(1,3) - t(1,1) ) * ( t(2,2) - t(2,1) ) )

  a = ( - ( t(2,3) - t(2,1) ) * t(1,1)  &
        + ( t(1,3) - t(1,1) ) * t(2,1) ) / g

  b =     ( t(2,3) - t(2,1) ) / g

  c =   - ( t(1,3) - t(1,1) ) / g

  d = (   ( t(2,2) - t(2,1) ) * t(1,1) &
        - ( t(1,2) - t(1,1) ) * t(2,1) ) / g

  e =   - ( t(2,2) - t(2,1) ) / g

  f =     ( t(1,2) - t(1,1) ) / g

  return
end