subroutine triangle_grid ( n, t, tg )  

!*****************************************************************************80
!
!! triangle_grid() computes points on a triangular grid.
!
!  Discussion:
!
!    The grid is defined by specifying the coordinates of an enclosing
!    triangle T, and the number of subintervals each side of the triangle
!    should be divided into.
!
!    Choosing N = 10, for instance, breaks each side into 10 subintervals,
!    and produces a grid of ((10+1)*(10+2))/2 = 66 points.
!
!              X
!             9 X
!            8 9 X
!           7 8 9 X
!          6 7 8 9 X
!         5 6 7 8 9 X
!        4 5 6 7 8 9 X
!       3 4 5 6 7 8 9 X
!      2 3 4 5 6 7 8 9 X
!     1 2 3 4 5 6 7 8 9 X
!    0 1 2 3 4 5 6 7 8 9 X
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 September 2010
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of subintervals.
!
!    Input, real ( kind = rk ) T(2,3), the coordinates of the points
!    defining the triangle.
!
!    Output, real ( kind = rk ) TG(2,((N+1)*(N+2))/2), the coordinates
!    of the points in the triangle.
!
  implicit none

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

  integer n

  integer i
  real ( kind = rk ) ir
  integer j
  real ( kind = rk ) jr
  integer k
  real ( kind = rk ) kr
  real ( kind = rk ) nr
  integer p
  real ( kind = rk ) t(2,3)
  real ( kind = rk ) tg(2,((n+1)*(n+2))/2)

  p = 0
  nr = real ( n, kind = rk )

  do i = 0, n
    ir = real ( i, kind = rk )
    do j = 0, n - i
      jr = real ( j, kind = rk )
      k = n - i - j
      kr = real ( k, kind = rk )
      p = p + 1
      tg(1:2,p) = ( ir * t(1:2,1) + jr * t(1:2,2) + kr * t(1:2,3) ) / nr
    end do
  end do

  return
end