#! /usr/bin/env python
#
def tetrahedron_grid_points ( n, xv, ng ):

#*****************************************************************************80
#
## TETRAHEDRON_GRID_POINTs computes points on a tetrahedral grid.
#
#  Discussion:
#
#    The grid is defined by specifying the coordinates of an enclosing
#    tetrahedron T, and the number of subintervals each edge of the 
#    tetrahedron should be divided into.
#
#    Choosing N = 10, for instance, breaks each side into 10 subintervals,
#    and produces a grid of ((10+1)*(10+2)*(10+3))/6 = 286 points.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    15 April 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the number of subintervals.
#
#    Input, real XV[4,3], the vertices of the tetrahedron.
#
#    Input, integer NG, the number of grid points.
#
#    Output, real XG[NG,3], the tetrahedron grid points.
#
  import numpy as np

  xg = np.zeros ( [ ng, 3 ] )

  p = 0

  for i in range ( 0, n + 1 ):
    for j in range ( 0, n + 1 - i ):
      for k in range ( 0, n + 1 - i - j ):
        l = n - i - j - k
        xg[p,0] = ( float ( i ) * xv[0,0] + float ( j ) * xv[1,0] \
                  + float ( k ) * xv[2,0] + float ( l ) * xv[3,0] ) / float ( n )
        xg[p,1] = ( float ( i ) * xv[0,1] + float ( j ) * xv[1,1] \
                  + float ( k ) * xv[2,1] + float ( l ) * xv[3,1] ) / float ( n )
        xg[p,2] = ( float ( i ) * xv[0,2] + float ( j ) * xv[1,2] \
                  + float ( k ) * xv[2,2] + float ( l ) * xv[3,2] ) / float ( n )
        p = p + 1

  return xg

def tetrahedron_grid_points_test ( ):

#*****************************************************************************80
#
## TETRAHEDRON_GRID_POINTS_TEST tests TETRAHEDRON_GRID_POINTS.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    10 November 2011
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from r83col_print_part import r83col_print_part
  from r8mat_write import r8mat_write
  from tetrahedron_grid_count import tetrahedron_grid_count
  from tetrahedron_grid_display import tetrahedron_grid_display

  print ( '' )
  print ( 'TETRAHEDRON_GRID_POINTS_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  TETRAHEDRON_GRID_POINTS can define a tetrahedral grid' )
  print ( '  with N+1 points on a side, based on any tetrahedron.' )

  n = 5
  print ( '  N = %d' % ( n ) )

  ng = tetrahedron_grid_count ( n )
#
#  Define the tetrahedron vertices.
#
  xv = np.array ( [ \
    [ 0.0, 0.0, 0.0 ], \
    [ 2.0, 1.0, 0.0 ], \
    [ 1.0, 4.0, 0.0 ], \
    [ 3.0, 3.0, 3.0 ] ] )

  r83col_print_part ( 4, xv, 20, '  Tetrahedron vertices:' )

  xg = tetrahedron_grid_points ( n, xv, ng )

  r83col_print_part ( ng, xg, 20, '  Tetrahedron grid points:' )
#
#  Write the data to a file.
#
  filename = 'tetrahedron_grid_points.xyz'

  r8mat_write ( filename, ng, 3, xg )

  print ( '' )
  print ( '  Data written to the file "%s".' % ( filename ) )
#
#  Plot the data.
#
  filename = 'tetrahedron_grid_points.png'
  tetrahedron_grid_display ( xv, ng, xg, filename )

  print ( '' )
  print ( '  Plot written to the file "%s".' % ( filename ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'TETRAHEDRON_GRID_POINTS_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  tetrahedron_grid_points_test ( )
  timestamp ( )