#! /usr/bin/env python3
#
def padua_order ( l ):

#*****************************************************************************80
#
## padua_order() returns the size of the Padua set of given level.
#
#  Discussion:
#
#    The Padua sets are indexed by a level that starts at 0.
#    This function returns the number of points in each level.
#
#  Example:
#
#    Level   Size
#    -----   ----
#        0      1
#        1      3
#        2      6
#        3     10
#        4     15
#        5     21
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Marco Caliari, Stefano de Marchi, Marco Vianello,
#    Bivariate interpolation on the square at new nodal sets,
#    Applied Mathematics and Computation,
#    Volume 165, Number 2, 2005, pages 261-274.
#
#  Input:
#
#    integer L, the level of the set.
#    0 <= L
#
#  Output:
#
#    integer N, the order (number of points) in the set.
#
  n = 0

  for i in range ( 0, l + 1 ):
    n = n + ( l // 2 ) + 1
    if ( ( l % 2 ) == 1 and ( i % 2 ) == 1 ):
      n = n + 1

  return n

def padua_order_test ( ):

#*****************************************************************************80
#
## padua_order_test() tests padua_order().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'padua_order_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  padua_order converts the level L into the order N' )
  print ( '  of any Padua rule.' )
  print ( '' )
  print ( '     L         N' )
  print ( '' )

  for l in range ( 0, 11 ):
    n = padua_order ( l )
    print ( '  %4d  %8d' % ( l, n ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'padua_order_test:' )
  print ( '  Normal end of execution.' )
  return

def padua_plot ( l, filename ):

#*****************************************************************************80
#
## padua_plot() plots the Padua points of given level.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    01 September 2016
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer L, the level of the set.
#    0 <= L
#
#    string FILENAME, the name for the PNG file to be created.
#
  import matplotlib.pyplot as plt

  n = padua_order ( l )

  xy = padua_points ( l )

  plt.plot ( xy[0,0:n], xy[1,0:n], 'b.', markersize = 30 )
  plt.axis ( [ -1.05, +1.05, -1.05, +1.05 ] )
  plt.grid ( True )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  label = 'Padua Points, Level %d' % ( l )
  plt.title ( label )
  plt.gca().set_aspect ( 'equal', adjustable = 'box' )

  plt.savefig ( filename )
  plt.show ( block = False )
  plt.close ( )

  return

def padua_plot_test ( ):

#*****************************************************************************80
#
## padua_plot_test() tests padua_plot().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    01 September 2016
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'padua_plot_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  padua_plot plots the Padua points.' )
  print ( '' )

  for l in range ( 0, 11 ):
    filename = 'padua_' + str ( l ) + '.png'
    padua_plot ( l, filename )
    print ( '  Plot file stored as "%s"' % ( filename ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'padua_plot_test' )
  print ( '  Normal end of execution' )
  return

def padua_points ( l ):

#*****************************************************************************80
#
## padua_points() returns the Padua points of level L.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Marco Caliari, Stefano de Marchi, Marco Vianello,
#    Bivariate interpolation on the square at new nodal sets,
#    Applied Mathematics and Computation,
#    Volume 165, Number 2, 2005, pages 261-274.
#
#  Input:
#
#    integer L, the level of the set.
#    0 <= L
#
#  Output:
#
#    real XY(2,N), the Padua points of level L.
#
  import numpy as np

  n = ( ( l + 1 ) * ( l + 2 ) ) // 2

  xy = np.zeros ( [ 2, n ] )

  if ( 0 < l ):

    k = 0

    for i in range ( 0, l + 1 ):

      j_hi = ( l // 2 ) + 1

      if ( ( l % 2 ) == 1 and ( i % 2 ) == 1 ):
        j_hi = j_hi + 1

      for j in range ( 1, j_hi + 1 ):

        if ( i * 2 == l ):
          xy[0,k] = 0.0
        else:
          angle1 = float ( i ) * np.pi / float ( l )
          xy[0,k] = np.cos ( angle1 )

        if ( ( i % 2 ) == 0 ):
          if ( 2 * ( 2 * j - 1 ) == l + 1 ):
            xy[1,k] = 0.0
          else:
            angle2 = float ( 2 * j - 1 ) * np.pi / float ( l + 1 )
            xy[1,k] = np.cos ( angle2 )
        else:
          if ( 2 * ( 2 * j - 2 ) == l + 1 ):
            xy[1,k] = 0.0
          else:
            angle2 = float ( 2 * j - 2 ) * np.pi / float ( l + 1 )
            xy[1,k] = np.cos ( angle2 )

        k = k + 1

  return xy

def padua_points_test ( ):

#*****************************************************************************80
#
## padua_points_test() tests padua_points().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'padua_points_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  padua_points returns the points of a Padua rule.' )

  for l in range ( 0, 11 ):
    n = padua_order ( l )
    xy = padua_points ( l )
    label = '  Level %d Padua points:' % ( l )
    r8mat_transpose_print ( 2, n, xy, label )
#
#  Terminate.
#
  print ( '' )
  print ( 'padua_points_test:' )
  print ( '  Normal end of execution.' )
  return

def padua_points_set ( l ):

#*****************************************************************************80
#
## padua_points_set() sets the Padua points.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 August 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Marco Caliari, Stefano de Marchi, Marco Vianello,
#    Bivariate interpolation on the square at new nodal sets,
#    Applied Mathematics and Computation,
#    Volume 165, Number 2, 2005, pages 261-274.
#
#  Input:
#
#    integer L, the level.
#    0 <= L <= 10.
#
#  Output:
#
#    real XY[2,N], the Padua points.
#
  import numpy as np

  if ( l == 0 ):

    xy = np.array ( [ \
      [ 0.000000000000000, 0.000000000000000 ] \
    ])
 
  elif ( l == 1 ):

    xy = np.array ( [ \
      [ -1.000000000000000, -1.000000000000000 ], \
      [ -1.000000000000000,  1.000000000000000 ], \
      [  1.000000000000000,  0.000000000000000 ] \
    ] )

  elif ( l == 2 ):

    xy = np.array ( [ \
      [ -1.000000000000000, -1.000000000000000 ], \
      [ -1.000000000000000,  0.5000000000000001 ], \
      [  0.000000000000000, -0.4999999999999998 ], \
      [  0.000000000000000,  1.000000000000000 ], \
      [  1.000000000000000, -1.000000000000000 ], \
      [  1.000000000000000,  0.5000000000000001 ] \
    ] )
 
  elif ( l == 3 ):

    xy = np.array ( [ \
      [ -1.000000000000000,  -1.000000000000000 ], \
      [ -1.000000000000000,   0.000000000000000 ], \
      [ -1.000000000000000,   1.000000000000000 ], \
      [ -0.4999999999999998, -0.7071067811865475 ], \
      [ -0.4999999999999998,  0.7071067811865476 ], \
      [  0.5000000000000001, -1.000000000000000 ], \
      [  0.5000000000000001,  0.000000000000000 ], \
      [  0.5000000000000001,  1.000000000000000 ], \
      [  1.000000000000000,  -0.7071067811865475 ], \
      [  1.000000000000000,   0.7071067811865476 ] \
    ] )
  
  elif ( l == 4 ):

    xy = np.array ( [ \
      [ -1.000000000000000,  -1.000000000000000 ], \
      [ -1.000000000000000,  -0.3090169943749473 ], \
      [ -1.000000000000000,   0.8090169943749475 ], \
      [ -0.7071067811865475, -0.8090169943749473 ], \
      [ -0.7071067811865475,  0.3090169943749475 ], \
      [ -0.7071067811865475,  1.000000000000000 ], \
      [  0.000000000000000,  -1.000000000000000 ], \
      [  0.000000000000000,  -0.3090169943749473 ], \
      [  0.000000000000000,   0.8090169943749475 ], \
      [  0.7071067811865476, -0.8090169943749473 ], \
      [  0.7071067811865476,  0.3090169943749475 ], \
      [  0.7071067811865476,  1.000000000000000 ], \
      [  1.000000000000000,  -1.000000000000000 ], \
      [  1.000000000000000,  -0.3090169943749473 ], \
      [  1.000000000000000,   0.8090169943749475 ] \
    ] )

  elif ( l == 5 ):

    xy = np.array ( [ \
      [ -1.000000000000000,  -1.000000000000000 ], \
      [ -1.000000000000000,  -0.4999999999999998 ], \
      [ -1.000000000000000,   0.5000000000000001 ], \
      [ -1.000000000000000,   1.000000000000000 ], \
      [ -0.8090169943749473, -0.8660254037844387 ], \
      [ -0.8090169943749473,  0.000000000000000 ], \
      [ -0.8090169943749473,  0.8660254037844387 ], \
      [ -0.3090169943749473, -1.000000000000000 ], \
      [ -0.3090169943749473, -0.4999999999999998 ], \
      [ -0.3090169943749473,  0.5000000000000001 ], \
      [ -0.3090169943749473,  1.000000000000000 ], \
      [  0.3090169943749475, -0.8660254037844387 ], \
      [  0.3090169943749475,  0.000000000000000 ], \
      [  0.3090169943749475,  0.8660254037844387 ], \
      [  0.8090169943749475, -1.000000000000000 ], \
      [  0.8090169943749475, -0.4999999999999998 ], \
      [  0.8090169943749475,  0.5000000000000001 ], \
      [  0.8090169943749475,  1.000000000000000 ], \
      [  1.000000000000000,  -0.8660254037844387 ], \
      [  1.000000000000000,   0.000000000000000 ], \
      [  1.000000000000000,   0.8660254037844387 ] \
    ] )

  elif ( l == 6 ):   

    xy = np.array ( [ \
      [ -1.000000000000000,  -1.000000000000000 ], \
      [ -1.000000000000000,  -0.6234898018587335 ], \
      [ -1.000000000000000,   0.2225209339563144 ], \
      [ -1.000000000000000,   0.9009688679024191 ], \
      [ -0.8660254037844387, -0.9009688679024190 ], \
      [ -0.8660254037844387, -0.2225209339563143 ], \
      [ -0.8660254037844387,  0.6234898018587336 ], \
      [ -0.8660254037844387,  1.000000000000000 ], \
      [ -0.4999999999999998, -1.000000000000000 ], \
      [ -0.4999999999999998, -0.6234898018587335 ], \
      [ -0.4999999999999998,  0.2225209339563144 ], \
      [ -0.4999999999999998,  0.9009688679024191 ], \
      [  0.000000000000000,  -0.9009688679024190 ], \
      [  0.000000000000000,  -0.2225209339563143 ], \
      [  0.000000000000000,   0.6234898018587336 ], \
      [  0.000000000000000,   1.000000000000000 ], \
      [  0.5000000000000001, -1.000000000000000 ], \
      [  0.5000000000000001, -0.6234898018587335 ], \
      [  0.5000000000000001,  0.2225209339563144 ], \
      [  0.5000000000000001,  0.9009688679024191 ], \
      [  0.8660254037844387, -0.9009688679024190 ], \
      [  0.8660254037844387, -0.2225209339563143 ], \
      [  0.8660254037844387,  0.6234898018587336 ], \
      [  0.8660254037844387,  1.000000000000000 ], \
      [  1.000000000000000,  -1.000000000000000 ], \
      [  1.000000000000000,  -0.6234898018587335 ], \
      [  1.000000000000000,   0.2225209339563144 ], \
      [  1.000000000000000,   0.9009688679024191 ] \
    ] )

  elif ( l == 7 ):    

    xy = np.array ( [ \
      [ -1.000000000000000, -1.000000000000000 ], \
      [ -1.000000000000000, -0.7071067811865475 ], \
      [ -1.000000000000000,  0.000000000000000 ], \
      [ -1.000000000000000,  0.7071067811865476 ], \
      [ -1.000000000000000,  1.000000000000000 ], \
      [ -0.9009688679024190, -0.9238795325112867 ], \
      [ -0.9009688679024190, -0.3826834323650897 ], \
      [ -0.9009688679024190,  0.3826834323650898 ], \
      [ -0.9009688679024190,  0.9238795325112867 ], \
      [ -0.6234898018587335, -1.000000000000000 ], \
      [ -0.6234898018587335, -0.7071067811865475 ], \
      [ -0.6234898018587335,  0.000000000000000 ], \
      [ -0.6234898018587335,  0.7071067811865476 ], \
      [ -0.6234898018587335, 1.000000000000000 ], \
      [ -0.2225209339563143, -0.9238795325112867 ], \
      [ -0.2225209339563143, -0.3826834323650897 ], \
      [ -0.2225209339563143,  0.3826834323650898 ], \
      [ -0.2225209339563143,  0.9238795325112867 ], \
      [  0.2225209339563144, -1.000000000000000 ], \
      [  0.2225209339563144, -0.7071067811865475 ], \
      [  0.2225209339563144,  0.000000000000000 ], \
      [  0.2225209339563144,  0.7071067811865476 ], \
      [  0.2225209339563144,  1.000000000000000 ], \
      [  0.6234898018587336, -0.9238795325112867 ], \
      [  0.6234898018587336, -0.3826834323650897 ], \
      [  0.6234898018587336,  0.3826834323650898 ], \
      [  0.6234898018587336,  0.9238795325112867 ], \
      [  0.9009688679024191, -1.000000000000000 ], \
      [  0.9009688679024191, -0.7071067811865475 ], \
      [  0.9009688679024191,  0.000000000000000 ], \
      [  0.9009688679024191,  0.7071067811865476 ], \
      [  0.9009688679024191,  1.000000000000000 ], \
      [  1.000000000000000, -0.9238795325112867 ], \
      [  1.000000000000000, -0.3826834323650897 ], \
      [  1.000000000000000,  0.3826834323650898 ], \
      [  1.000000000000000,  0.9238795325112867 ] \
    ] )

  elif ( l == 8 ): 

    xy = np.array ( [ \
      [ -1.000000000000000, -1.000000000000000 ], \
      [ -1.000000000000000, -0.7660444431189779 ], \
      [ -1.000000000000000, -0.1736481776669303 ], \
      [ -1.000000000000000,  0.5000000000000001 ], \
      [ -1.000000000000000,  0.9396926207859084 ], \
      [ -0.9238795325112867, -0.9396926207859083 ], \
      [ -0.9238795325112867, -0.4999999999999998 ], \
      [ -0.9238795325112867,  0.1736481776669304 ], \
      [ -0.9238795325112867,  0.7660444431189780 ], \
      [ -0.9238795325112867,  1.000000000000000 ], \
      [ -0.7071067811865475, -1.000000000000000 ], \
      [ -0.7071067811865475, -0.7660444431189779 ], \
      [ -0.7071067811865475, -0.1736481776669303 ], \
      [ -0.7071067811865475,  0.5000000000000001 ], \
      [ -0.7071067811865475,  0.9396926207859084 ], \
      [ -0.3826834323650897, -0.9396926207859083 ], \
      [ -0.3826834323650897, -0.4999999999999998 ], \
      [ -0.3826834323650897,  0.1736481776669304 ], \
      [ -0.3826834323650897,  0.7660444431189780 ], \
      [ -0.3826834323650897,  1.000000000000000 ], \
      [  0.000000000000000, -1.000000000000000 ], \
      [  0.000000000000000, -0.7660444431189779 ], \
      [  0.000000000000000, -0.1736481776669303 ], \
      [  0.000000000000000,  0.5000000000000001 ], \
      [  0.000000000000000,  0.9396926207859084 ], \
      [  0.3826834323650898, -0.9396926207859083 ], \
      [  0.3826834323650898, -0.4999999999999998 ], \
      [  0.3826834323650898,  0.1736481776669304 ], \
      [  0.3826834323650898,  0.7660444431189780 ], \
      [  0.3826834323650898,  1.000000000000000 ], \
      [  0.7071067811865476, -1.000000000000000 ], \
      [  0.7071067811865476, -0.7660444431189779 ], \
      [  0.7071067811865476, -0.1736481776669303 ], \
      [  0.7071067811865476,  0.5000000000000001 ], \
      [  0.7071067811865476,  0.9396926207859084 ], \
      [  0.9238795325112867, -0.9396926207859083 ], \
      [  0.9238795325112867, -0.4999999999999998 ], \
      [  0.9238795325112867,  0.1736481776669304 ], \
      [  0.9238795325112867,  0.7660444431189780 ], \
      [  0.9238795325112867,  1.000000000000000 ], \
      [  1.000000000000000, -1.000000000000000 ], \
      [  1.000000000000000, -0.7660444431189779 ], \
      [  1.000000000000000, -0.1736481776669303 ], \
      [  1.000000000000000,  0.5000000000000001 ], \
      [  1.000000000000000,  0.9396926207859084 ] \
    ] )

  elif ( l == 9 ): 

    xy = np.array ( [ \
      [ -1.000000000000000, -1.000000000000000 ], \
      [ -1.000000000000000, -0.8090169943749473 ], \
      [ -1.000000000000000, -0.3090169943749473 ], \
      [ -1.000000000000000,  0.3090169943749475 ], \
      [ -1.000000000000000,  0.8090169943749475 ], \
      [ -1.000000000000000,  1.000000000000000 ], \
      [ -0.9396926207859083, -0.9510565162951535 ], \
      [ -0.9396926207859083, -0.5877852522924730 ], \
      [ -0.9396926207859083,  0.000000000000000 ], \
      [ -0.9396926207859083,  0.5877852522924731 ], \
      [ -0.9396926207859083,  0.9510565162951535 ], \
      [ -0.7660444431189779, -1.000000000000000 ], \
      [ -0.7660444431189779, -0.8090169943749473 ], \
      [ -0.7660444431189779, -0.3090169943749473 ], \
      [ -0.7660444431189779,  0.3090169943749475 ], \
      [ -0.7660444431189779,  0.8090169943749475 ], \
      [ -0.7660444431189779,  1.000000000000000 ], \
      [ -0.4999999999999998, -0.9510565162951535 ], \
      [ -0.4999999999999998, -0.5877852522924730 ], \
      [ -0.4999999999999998,  0.000000000000000 ], \
      [ -0.4999999999999998,  0.5877852522924731 ], \
      [ -0.4999999999999998,  0.9510565162951535 ], \
      [ -0.1736481776669303, -1.000000000000000 ], \
      [ -0.1736481776669303, -0.8090169943749473 ], \
      [ -0.1736481776669303, -0.3090169943749473 ], \
      [ -0.1736481776669303,  0.3090169943749475 ], \
      [ -0.1736481776669303,  0.8090169943749475 ], \
      [ -0.1736481776669303,  1.000000000000000 ], \
      [  0.1736481776669304, -0.9510565162951535 ], \
      [  0.1736481776669304, -0.5877852522924730 ], \
      [  0.1736481776669304,  0.000000000000000 ], \
      [  0.1736481776669304,  0.5877852522924731 ], \
      [  0.1736481776669304,  0.9510565162951535 ], \
      [  0.5000000000000001, -1.000000000000000 ], \
      [  0.5000000000000001, -0.8090169943749473 ], \
      [  0.5000000000000001, -0.3090169943749473 ], \
      [  0.5000000000000001,  0.3090169943749475 ], \
      [  0.5000000000000001,  0.8090169943749475 ], \
      [  0.5000000000000001,  1.000000000000000 ], \
      [  0.7660444431189780, -0.9510565162951535 ], \
      [  0.7660444431189780, -0.5877852522924730 ], \
      [  0.7660444431189780,  0.000000000000000 ], \
      [  0.7660444431189780,  0.5877852522924731 ], \
      [  0.7660444431189780,  0.9510565162951535 ], \
      [  0.9396926207859084, -1.000000000000000 ], \
      [  0.9396926207859084, -0.8090169943749473 ], \
      [  0.9396926207859084, -0.3090169943749473 ], \
      [  0.9396926207859084,  0.3090169943749475 ], \
      [  0.9396926207859084,  0.8090169943749475 ], \
      [  0.9396926207859084,  1.000000000000000 ], \
      [  1.000000000000000, -0.9510565162951535 ], \
      [  1.000000000000000, -0.5877852522924730 ], \
      [  1.000000000000000,  0.000000000000000 ], \
      [  1.000000000000000,  0.5877852522924731 ], \
      [  1.000000000000000,  0.9510565162951535 ] \
    ] )

  elif ( l == 10 ):

    xy = np.array ( [ \
      [ -1.000000000000000, -1.000000000000000 ], \
      [ -1.000000000000000, -0.8412535328311811 ], \
      [ -1.000000000000000, -0.4154150130018863 ], \
      [ -1.000000000000000,  0.1423148382732851 ], \
      [ -1.000000000000000,  0.6548607339452851 ], \
      [ -1.000000000000000,  0.9594929736144974 ], \
      [ -0.9510565162951535, -0.9594929736144974 ], \
      [ -0.9510565162951535, -0.6548607339452850 ], \
      [ -0.9510565162951535, -0.1423148382732850 ], \
      [ -0.9510565162951535,  0.4154150130018864 ], \
      [ -0.9510565162951535,  0.8412535328311812 ], \
      [ -0.9510565162951535,  1.000000000000000 ], \
      [ -0.8090169943749473, -1.000000000000000 ], \
      [ -0.8090169943749473, -0.8412535328311811 ], \
      [ -0.8090169943749473, -0.4154150130018863 ], \
      [ -0.8090169943749473,  0.1423148382732851 ], \
      [ -0.8090169943749473,  0.6548607339452851 ], \
      [ -0.8090169943749473,  0.9594929736144974 ], \
      [ -0.5877852522924730, -0.9594929736144974 ], \
      [ -0.5877852522924730, -0.6548607339452850 ], \
      [ -0.5877852522924730, -0.1423148382732850 ], \
      [ -0.5877852522924730,  0.4154150130018864 ], \
      [ -0.5877852522924730,  0.8412535328311812 ], \
      [ -0.5877852522924730,  1.000000000000000 ], \
      [ -0.3090169943749473, -1.000000000000000 ], \
      [ -0.3090169943749473, -0.8412535328311811 ], \
      [ -0.3090169943749473, -0.4154150130018863 ], \
      [ -0.3090169943749473,  0.1423148382732851 ], \
      [ -0.3090169943749473,  0.6548607339452851 ], \
      [ -0.3090169943749473, 0.9594929736144974 ], \
      [  0.000000000000000, -0.9594929736144974 ], \
      [  0.000000000000000, -0.6548607339452850 ], \
      [  0.000000000000000, -0.1423148382732850 ], \
      [  0.000000000000000,  0.4154150130018864 ], \
      [  0.000000000000000,  0.8412535328311812 ], \
      [  0.000000000000000,  1.000000000000000 ], \
      [  0.3090169943749475, -1.000000000000000 ], \
      [  0.3090169943749475, -0.8412535328311811 ], \
      [  0.3090169943749475, -0.4154150130018863 ], \
      [  0.3090169943749475,  0.1423148382732851 ], \
      [  0.3090169943749475,  0.6548607339452851 ], \
      [  0.3090169943749475,  0.9594929736144974 ], \
      [  0.5877852522924731, -0.9594929736144974 ], \
      [  0.5877852522924731, -0.6548607339452850 ], \
      [  0.5877852522924731, -0.1423148382732850 ], \
      [  0.5877852522924731,  0.4154150130018864 ], \
      [  0.5877852522924731,  0.8412535328311812 ], \
      [  0.5877852522924731, 1.000000000000000 ], \
      [  0.8090169943749475, -1.000000000000000 ], \
      [  0.8090169943749475, -0.8412535328311811 ], \
      [  0.8090169943749475, -0.4154150130018863 ], \
      [  0.8090169943749475,  0.1423148382732851 ], \
      [  0.8090169943749475,  0.6548607339452851 ], \
      [  0.8090169943749475,  0.9594929736144974 ], \
      [  0.9510565162951535, -0.9594929736144974 ], \
      [  0.9510565162951535, -0.6548607339452850 ], \
      [  0.9510565162951535, -0.1423148382732850 ], \
      [  0.9510565162951535,  0.4154150130018864 ], \
      [  0.9510565162951535,  0.8412535328311812 ], \
      [  0.9510565162951535,  1.000000000000000 ], \
      [  1.000000000000000, -1.000000000000000 ], \
      [  1.000000000000000, -0.8412535328311811 ], \
      [  1.000000000000000, -0.4154150130018863 ], \
      [  1.000000000000000,  0.1423148382732851 ], \
      [  1.000000000000000,  0.6548607339452851 ], \
      [  1.000000000000000,  0.9594929736144974 ] \
    ] )

  else:

    print ( '' )
    print ( 'padua_points_set - Fatal error!' )
    print ( '  Illegal value of L = %d' % ( l ) )
    print ( '  Legal values are 1 through 10.' )
    raise Exception ( 'padua_points_set - Fatal error!' )
#
#  Transpose the array.
#
  xy = xy.transpose ( )
#
#  Reverse the indexing to match the published data.
#
  n = ( ( l + 1 ) * ( l + 2 ) ) // 2
  j1_hi = ( ( n - 1 ) // 2 )

  for j1 in range ( 0, j1_hi + 1 ):

    j2 = n - 1 - j1
    t        = xy[0,j1]
    xy[0,j1] = xy[0,j2]
    xy[0,j2] = t
    t        = xy[1,j1]
    xy[1,j1] = xy[1,j2]
    xy[1,j2] = t

  return xy

def padua_points_set_test ( ):

#*****************************************************************************80
#
## padua_points_set_test() tests padua_points_set().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 August 2016
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'padua_points_set_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  padua_points_set looks the Padua points up in a table.' )
 
  for l in range ( 3, 5 ):
    n = padua_order ( l )
    xy1 = padua_points ( l )
    xy2 = padua_points_set ( l )
    print ( '' )
    print ( '  Level %d  Padua points' % ( l ) )
    print ( '' )
    for j in range ( 0, n ):
      print ( '  %4d  %14.6g  %14.6g' % ( j, xy1[0,j], xy1[1,j] ) )
      print ( '        %14.6g  %14.6g' % (    xy2[0,j], xy2[1,j] ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'padua_points_set_test:' )
  print ( '  Normal end of execution.' )
  return

def padua_test ( ):

#*****************************************************************************80
#
## padua_test() tests padua().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    23 March 2021
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'padua_test:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test padua()' )

  padua_order_test ( )
  padua_plot_test ( )
  padua_points_test ( )
  padua_points_set_test ( )
  padua_weights_test ( )
  padua_weights_set_test ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'padua_test:' )
  print ( '  Normal end of execution.' )
  return

def padua_weights ( l ):

#*****************************************************************************80
#
## padua_weights() returns quadrature weights for Padua points.
#
#  Discussion:
#
#    The order of the weights corresponds to the ordering used
#    by the companion function padua_points().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 August 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Marco Caliari, Stefano de Marchi, Marco Vianello,
#    Bivariate interpolation on the square at new nodal sets,
#    Applied Mathematics and Computation,
#    Volume 165, Number 2, 2005, pages 261-274.
#
#  Input:
#
#    integer L, the level of the set.
#    0 <= L
#
#  Output:
#
#    real W((L+1)*(L+2)/2), the quadrature weights.
#
  import numpy as np

  n = ( ( l + 1 ) * ( l + 2 ) ) // 2
  w = np.zeros ( n )

  if ( l == 0 ):
    w[0] = 4.0
    return w
#
#  Relatives of L/2:
#
  lp1h = ( ( l + 1 ) // 2 )
  lp2h = ( ( l + 2 ) // 2 )
  lp3h = ( ( l + 3 ) // 2 )
#
#  TE1, TE2, TO1, TO2: 
#  Even and odd Chebyshev polynomials on subgrids 1 and 2.
#
  te1 = np.zeros ( [ lp2h, lp2h ] )
  
  for j in range ( 0, lp2h ):
    for i in range ( 0, lp2h ):
      te1[i,j] = \
        np.cos ( float ( 2 * i * 2 * j ) * np.pi / float ( l ) )

  te1[1:lp2h,0:lp2h] = te1[1:lp2h,0:lp2h] * np.sqrt ( 2.0 )

  to1 = np.zeros ( [ lp2h, lp1h ] )

  for j in range ( 0, lp1h ):
    for i in range ( 0, lp2h ):
      to1[i,j] = \
        np.cos ( float ( 2 * i * ( 2 * j + 1 ) ) * np.pi / float ( l ) )

  to1[1:lp2h,0:lp1h] = to1[1:lp2h,0:lp1h] * np.sqrt ( 2.0 )

  te2 = np.zeros ( [ lp2h, lp3h ] )

  for j in range ( 0, lp3h ):
    for i in range ( 0, lp2h ):
      te2[i,j] = \
        np.cos ( float ( 2 * i * 2 * j ) * np.pi / float ( l + 1 ) )

  te2[1:lp2h,0:lp3h] = te2[1:lp2h,0:lp3h] * np.sqrt ( 2.0 )

  to2 = np.zeros ( [ lp2h, lp2h ] )

  for j in range ( 0, lp2h ):
    for i in range ( 0, lp2h ):
      to2[i,j] = \
        np.cos ( float ( 2 * i * ( 2 * j + 1 ) ) * np.pi / float ( l + 1 ) )

  to2[1:lp2h,0:lp2h] = to2[1:lp2h,0:lp2h] * np.sqrt ( 2.0 )
#
#  MOM: Moments matrix for even * even pairs.
#
  mom = np.zeros ( [ lp2h, lp2h ] )

  for j in range ( 0, lp2h ):
    mj = 2.0 * np.sqrt ( 2.0 ) / float ( 1.0 - ( 2 * j ) ** 2 )
    for i in range ( 0, lp2h - j ):
      mi = 2.0 * np.sqrt ( 2.0 ) / float ( 1.0 - ( 2 * i ) ** 2 )
      mom[i,j] = mi * mj

  mom[0,0:lp2h] = mom[0,0:lp2h] / np.sqrt ( 2.0 )
  mom[0:lp2h,0] = mom[0:lp2h,0] / np.sqrt ( 2.0 )

  if ( ( l % 2 ) == 0 ):
    mom[lp2h-1,0] = mom[lp2h-1,0] / 2.0
#
#  TMTOE and TMTEO: matrix products.
#
  tmtoe = np.zeros ( [ lp2h, lp2h ] )

  for j2 in range ( 0, lp2h ):
    for i2 in range ( 0, lp2h - j2 ):
      for j in range ( 0, lp2h ):
        for i in range ( 0, lp2h ):
          tmtoe[i,j] = tmtoe[i,j] + to2[i2,i] * mom[j2,i2] * te1[j2,j]

  tmteo = np.zeros ( [ lp3h, lp1h ] )

  for j2 in range ( 0, lp2h ):
    for i2 in range ( 0, lp2h - j2 ):
      for j in range ( 0, lp1h ):
        for i in range ( 0, lp3h ):
          tmteo[i,j] = tmteo[i,j] + te2[i2,i] * mom[j2,i2] * to1[j2,j]
#
#  W1 and W2: Interpolation weight matrices.
#
  w1 = 2.0 * np.ones ( [ lp2h, lp2h ] ) / float ( l * ( l + 1 ) )

  w1[0:lp2h,0] = w1[0:lp2h,0] / 2.0

  if ( ( l % 2 ) == 0 ):
    w1[0:lp2h,lp2h-1] = w1[0:lp2h,lp2h-1] / 2.0
    w1[lp2h-1,0:lp2h] = w1[lp2h-1,0:lp2h] / 2.0

  w2 = 2.0 * np.ones ( [ lp3h, lp1h ] ) / float ( l * ( l + 1 ) )
  w2[0,0:lp1h] = w2[0,0:lp1h] / 2.0

  if ( ( l % 2 ) == 1 ):
    w2[lp3h-1,0:lp1h] = w2[lp3h-1,0:lp1h] / 2.0
    w2[0:lp3h,lp1h-1] = w2[0:lp3h,lp1h-1] / 2.0
#
#  Cubature weights as matrices on the subgrids.
#
  for j in range ( 0, lp2h ):
    for i in range ( 0, lp2h ):
      w1[i,j] = w1[i,j] * tmtoe[i,j]

  for j in range ( 0, lp1h ):
    for i in range ( 0, lp3h ):
      w2[i,j] = w2[i,j] * tmteo[i,j]
#
#  Pack weight matrices W1 and W2 into the vector W.
#
  n = ( ( l + 1 ) * ( l + 2 ) ) // 2
  w = np.zeros ( n )

  if ( ( l % 2 ) == 0 ):

    for i in range ( 0, lp2h ):
      for j in range ( 0, lp2h ):
        w[i+2*j*lp2h] = w1[i,j]

    for i in range ( 0, lp3h ):
      for j in range ( 0, lp1h ):
        w[i+(2*j+1)*lp2h] = w2[i,j]

  else:

    for j in range ( 0, lp1h ):
      for i in range ( 0, lp2h ):
        w[i+j*(l+2)] = w1[i,j]

    for j in range ( 0, lp1h ):
      for i in range ( 0, lp3h ):
        w[i+lp2h+j*(l+2)] = w2[i,j]

  return w

def padua_weights_test ( ):

#*****************************************************************************80
#
## padua_weights_test() tests padua_weights().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 August 2016
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'padua_weights_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  padua_weights returns the weights of a Padua rule.' )

  for l in range ( 0, 11 ):
    n = padua_order ( l )
    w = padua_weights ( l )
    label = '  Level %d Padua weightss:' % ( l )
    r8vec_print ( n, w, label )
#
#  Terminate.
#
  print ( '' )
  print ( 'padua_weights_test' )
  print ( '  Normal end of execution' )
  return

def padua_weights_set ( l ):

#*****************************************************************************80
#
## padua_weights_set() sets quadrature weights for the Padua points.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    01 September 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Marco Caliari, Stefano de Marchi, Marco Vianello,
#    Bivariate interpolation on the square at new nodal sets,
#    Applied Mathematics and Computation,
#    Volume 165, Number 2, 2005, pages 261-274.
#
#  Input:
#
#    integer L, the level.
#    0 <= L <= 10.
#
#  Output:
#
#    real W(N), the quadrature weights.
#
  import numpy as np

  if ( l == 0 ):
    w = np.array \
    ( [ \
      4.000000000000000E+00 \
    ] )
  elif ( l == 1 ):
    w = np.array \
    ( [ \
        1.000000000000000E+00, \
        1.000000000000000E+00, \
        2.000000000000000E+00 \
    ] )
  elif ( l == 2 ):
    w = np.array \
    ( [ \
       0.0E+00, \
       0.6666666666666663E+00, \
       2.222222222222222E+00, \
       0.4444444444444444E+00, \
       0.0E+00, \
       0.6666666666666664E+00 \
    ] )
  elif ( l == 3 ):
    w = np.array \
    ( [ \
      -0.5555555555555480E-01, \
       0.3333333333333331E+00, \
      -0.5555555555555580E-01, \
       0.8888888888888886E+00, \
       0.8888888888888893E+00, \
       0.2222222222222224E+00, \
       1.333333333333333E+00, \
       0.2222222222222220E+00, \
       0.1111111111111109E+00, \
       0.1111111111111112E+00 \
    ] )
  elif ( l == 4 ):
    w = np.array \
    ( [ \
      -0.8888888888888932E-02, \
       0.8104919101110961E-01, \
       0.6117303121111219E-01, \
       0.3874097078666789E+00, \
       0.6259236254666545E+00, \
       0.5333333333333362E-01, \
       0.7111111111111067E-01, \
       0.9830822022444241E+00, \
       0.5458066866444642E+00, \
       0.3874097078666780E+00, \
       0.6259236254666568E+00, \
       0.5333333333333383E-01, \
      -0.8888888888888703E-02, \
       0.8104919101110968E-01, \
       0.6117303121111135E-01 \
    ] )
  elif ( l == 5 ):
    w = np.array \
    ( [ \
      -0.1037037037037093E-01, \
       0.5037037037036911E-01, \
       0.5037037037037081E-01, \
      -0.1037037037036947E-01, \
       0.1876963678740801E+00, \
       0.3460933466518654E+00, \
       0.1876963678740763E+00, \
       0.4514390511851724E-01, \
       0.5541130536814713E+00, \
       0.5541130536814728E+00, \
       0.4514390511851834E-01, \
       0.2804517802740705E+00, \
       0.6376103570518378E+00, \
       0.2804517802740683E+00, \
       0.3189313191851883E-01, \
       0.3288499092814910E+00, \
       0.3288499092814925E+00, \
       0.3189313191851956E-01, \
       0.2074074074074123E-01, \
       0.3851851851851849E-01, \
       0.2074074074074051E-01 \
    ] )
  elif ( l == 6 ):
    w = np.array \
    ( [ \
      -0.3023431594858565E-02, \
       0.1957267632451884E-01, \
       0.2633929313290840E-01, \
       0.1425431928029237E-01, \
       0.1006383046329639E+00, \
       0.2208900184526934E+00, \
       0.1743144584714012E+00, \
       0.1209372637943976E-01, \
       0.1934996220710680E-01, \
       0.3245064820875231E+00, \
       0.4027058473592984E+00, \
       0.1677234226317961E+00, \
       0.1953319357827178E+00, \
       0.4489633053035124E+00, \
       0.3721824611057551E+00, \
       0.2479213907785274E-01, \
       0.1934996220710561E-01, \
       0.3245064820875153E+00, \
       0.4027058473592959E+00, \
       0.1677234226317944E+00, \
       0.1006383046329745E+00, \
       0.2208900184526933E+00, \
       0.1743144584714027E+00, \
       0.1209372637944051E-01, \
      -0.3023431594861990E-02, \
       0.1957267632451757E-01, \
       0.2633929313290797E-01, \
       0.1425431928029198E-01 \
    ] )
  elif ( l == 7 ):
    w = np.array \
    ( [ \
      -0.3287981859413765E-02, \
       0.1337868480725671E-01, \
       0.2063492063491996E-01, \
       0.1337868480725546E-01, \
      -0.3287981859408898E-02, \
       0.5949324721885513E-01, \
       0.1306477599993571E+00, \
       0.1306477599993581E+00, \
       0.5949324721885061E-01, \
       0.1263869091685831E-01, \
       0.1979944935601103E+00, \
       0.2832184784823740E+00, \
       0.1979944935601143E+00, \
       0.1263869091685747E-01, \
       0.1221817987389771E+00, \
       0.3150266070593529E+00, \
       0.3150266070593440E+00, \
       0.1221817987389802E+00, \
       0.1771365352315134E-01, \
       0.2490926964598258E+00, \
       0.3408041116306980E+00, \
       0.2490926964598291E+00, \
       0.1771365352314976E-01, \
       0.9646986307476696E-01, \
       0.2557725606433917E+00, \
       0.2557725606433927E+00, \
       0.9646986307476431E-01, \
       0.8649923133686802E-02, \
       0.1062007918394705E+00, \
       0.1505805844901012E+00, \
       0.1062007918394705E+00, \
       0.8649923133690016E-02, \
       0.6355881462931014E-02, \
       0.1405228180237514E-01, \
       0.1405228180237651E-01, \
       0.6355881462928496E-02 \
    ] )
  elif ( l == 8 ):
    w = np.array \
    ( [ \
      -0.1269841269835311E-02, \
       0.6706089639041270E-02, \
       0.1111455441352989E-01, \
       0.1026455026455282E-01, \
       0.4930678698742625E-02, \
       0.3633146869162523E-01, \
       0.8838322767333079E-01, \
       0.9965911758463214E-01, \
       0.6400185533755555E-01, \
       0.4061629144893127E-02, \
       0.6772486772485166E-02, \
       0.1258344472781388E+00, \
       0.1927501398511116E+00, \
       0.1699470899470907E+00, \
       0.6342599488133535E-01, \
       0.8376332474107638E-01, \
       0.2170841444607031E+00, \
       0.2477307250801775E+00, \
       0.1648098048612226E+00, \
       0.1004771829779292E-01, \
       0.1015873015872910E-01, \
       0.1784328991205164E+00, \
       0.2729409493576765E+00, \
       0.2364021164021134E+00, \
       0.8936689226256009E-01, \
       0.8376332474107701E-01, \
       0.2170841444607054E+00, \
       0.2477307250801761E+00, \
       0.1648098048612200E+00, \
       0.1004771829779330E-01, \
       0.6772486772485237E-02, \
       0.1258344472781358E+00, \
       0.1927501398511135E+00, \
       0.1699470899470926E+00, \
       0.6342599488133838E-01, \
       0.3633146869162453E-01, \
       0.8838322767332588E-01, \
       0.9965911758463601E-01, \
       0.6400185533755502E-01, \
       0.4061629144888279E-02, \
      -0.1269841269836355E-02, \
       0.6706089639046927E-02, \
       0.1111455441352761E-01, \
       0.1026455026454956E-01, \
       0.4930678698747173E-02 \
    ] )
  elif ( l == 9 ):
    w = np.array \
    ( [ \
      -0.1368606701945113E-02, \
       0.4837977417140975E-02, \
       0.8876308297144902E-02, \
       0.8876308297143068E-02, \
       0.4837977417150492E-02, \
      -0.1368606701935084E-02, \
       0.2425285860992349E-01, \
       0.5727330842923516E-01, \
       0.7008257906578071E-01, \
       0.5727330842922034E-01, \
       0.2425285860989794E-01, \
       0.4659404339099723E-02, \
       0.8354521980498550E-01, \
       0.1370796991940044E+00, \
       0.1370796991940248E+00, \
       0.8354521980500107E-01, \
       0.4659404339109654E-02, \
       0.5564545640233619E-01, \
       0.1524391996823315E+00, \
       0.1877107583774149E+00, \
       0.1524391996823176E+00, \
       0.5564545640232402E-01, \
       0.8186176158691754E-02, \
       0.1295355639606716E+00, \
       0.2061407656847711E+00, \
       0.2061407656847630E+00, \
       0.1295355639606894E+00, \
       0.8186176158692687E-02, \
       0.6234969028097752E-01, \
       0.1730419031522391E+00, \
       0.2169418247419051E+00, \
       0.1730419031522361E+00, \
       0.6234969028097048E-01, \
       0.7506172839505762E-02, \
       0.1142161960569350E+00, \
       0.1802176663769002E+00, \
       0.1802176663769038E+00, \
       0.1142161960569279E+00, \
       0.7506172839512260E-02, \
       0.4031900987631698E-01, \
       0.1142976211857364E+00, \
       0.1413353845521477E+00, \
       0.1142976211857414E+00, \
       0.4031900987631700E-01, \
       0.3239075586856897E-02, \
       0.4317587564913915E-01, \
       0.7015250533601934E-01, \
       0.7015250533601930E-01, \
       0.4317587564913908E-01, \
       0.3239075586852207E-02, \
       0.2550690557469151E-02, \
       0.6084230077461027E-02, \
       0.7421516754852508E-02, \
       0.6084230077458821E-02, \
       0.2550690557473353E-02 \
    ] )
  elif ( l == 10 ):
    w = np.array \
    ( [ \
       -0.6240762604463766E-03, \
        0.2843227149025789E-02, \
        0.5250031948150784E-02, \
        0.5891746241568810E-02, \
        0.4705736485964679E-02, \
        0.2135354637732944E-02, \
        0.1610939653924566E-01, \
        0.4099595211758227E-01, \
        0.5326500934654063E-01, \
        0.4863338516658277E-01, \
        0.2843474741781434E-01, \
        0.1719619179693151E-02, \
        0.2883769745121509E-02, \
        0.5724711668876453E-01, \
        0.9659872841640438E-01, \
        0.1053210323353631E+00, \
        0.8066212502628711E-01, \
        0.2855765663647366E-01, \
        0.3981286043310814E-01, \
        0.1090390674981577E+00, \
        0.1430169021081585E+00, \
        0.1313686303763064E+00, \
        0.7932850918298831E-01, \
        0.4610696968783255E-02, \
        0.5086495679684716E-02, \
        0.9311356395361167E-01, \
        0.1562320334111262E+00, \
        0.1696057154254139E+00, \
        0.1283581371975154E+00, \
        0.4603059518094556E-01, \
        0.4894888812994630E-01, \
        0.1347281473526573E+00, \
        0.1764193542601264E+00, \
        0.1635037456303485E+00, \
        0.9822749154565460E-01, \
        0.5704840613923174E-02, \
        0.5086495679679268E-02, \
        0.9311356395362781E-01, \
        0.1562320334111511E+00, \
        0.1696057154253968E+00, \
        0.1283581371975113E+00, \
        0.4603059518094044E-01, \
        0.3981286043311782E-01, \
        0.1090390674981293E+00, \
        0.1430169021081508E+00, \
        0.1313686303763217E+00, \
        0.7932850918299997E-01, \
        0.4610696968790496E-02, \
        0.2883769745110260E-02, \
        0.5724711668875122E-01, \
        0.9659872841642343E-01, \
        0.1053210323353932E+00, \
        0.8066212502626474E-01, \
        0.2855765663644533E-01, \
        0.1610939653928420E-01, \
        0.4099595211758404E-01, \
        0.5326500934649123E-01, \
        0.4863338516656233E-01, \
        0.2843474741784810E-01, \
        0.1719619179720036E-02, \
       -0.6240762604606350E-03, \
        0.2843227149011163E-02, \
        0.5250031948172295E-02, \
        0.5891746241587802E-02, \
        0.4705736485965663E-02, \
        0.2135354637703863E-02 \
    ] )
  else:
    print ( '' )
    print ( 'padua_weights_set - Fatal error!' )
    print ( '  Illegal value of L = %d' % ( l ) )
    print ( '  Legal values are 0 through 10.' )
    raise Exception ( 'padua_weights_set - Fatal error!' )
#
#  Reverse the indexing to match the published data.
#
  n = ( ( l + 1 ) * ( l + 2 ) ) // 2
  w = r8vec_reverse ( n, w )

  return w

def padua_weights_set_test ( ):

#*****************************************************************************80
#
## padua_weights_set_test() tests padua_weights_set().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    01 September 2016
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'padua_weights_set_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  padua_weights_set looks up Padua weights in a table.' )
 
  for l in range ( 3, 5 ):
    n = padua_order ( l )
    w1 = padua_weights ( l )
    w2 = padua_weights_set ( l )
    print ( '' )
    print ( '  Level %d  Padua weights' % ( l ) )
    print ( '' )
    diff = 0.0
    for j in range ( 0, n ):
      diff = max ( diff, abs ( w1[j] - w2[j] ) )
      print ( '  %4d  %14.6g  %14.6g' % ( j, w1[j], w2[j] ) )
    print ( '' )
    print ( '  Maximum difference = %g' % ( diff ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'padua_weights_set_test' )
  print ( '  Normal end of execution.' )
  return

def r8mat_transpose_print ( m, n, a, title ):

#*****************************************************************************80
#
## r8mat_transpose_print() prints an R8MAT, transposed.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 August 2014
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer M, the number of rows in A.
#
#    integer N, the number of columns in A.
#
#    real A(M,N), the matrix.
#
#    string TITLE, a title.
#
  r8mat_transpose_print_some ( m, n, a, 0, 0, m - 1, n - 1, title )

  return

def r8mat_transpose_print_test ( ):

#*****************************************************************************80
#
## r8mat_transpose_print_test() tests r8mat_transpose_print().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 October 2014
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'r8mat_transpose_print_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  r8mat_transpose_print prints an R8MAT.' )

  m = 4
  n = 3
  v = np.array ( [ \
    [ 11.0, 12.0, 13.0 ], 
    [ 21.0, 22.0, 23.0 ], 
    [ 31.0, 32.0, 33.0 ], 
    [ 41.0, 42.0, 43.0 ] ], dtype = np.float64 )
  r8mat_transpose_print ( m, n, v, '  Here is an R8MAT, transposed:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'r8mat_transpose_print_test:' )
  print ( '  Normal end of execution.' )
  return

def r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ):

#*****************************************************************************80
#
## r8mat_transpose_print_some() prints a portion of an R8MAT, transposed.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    13 November 2014
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer M, N, the number of rows and columns of the matrix.
#
#    real A(M,N), an M by N matrix to be printed.
#
#    integer ILO, JLO, the first row and column to print.
#
#    integer IHI, JHI, the last row and column to print.
#
#    string TITLE, a title.
#
  incx = 5

  print ( '' )
  print ( title )

  if ( m <= 0 or n <= 0 ):
    print ( '' )
    print ( '  (None)' )
    return

  for i2lo in range ( max ( ilo, 0 ), min ( ihi, m - 1 ), incx ):

    i2hi = i2lo + incx - 1
    i2hi = min ( i2hi, m - 1 )
    i2hi = min ( i2hi, ihi )
    
    print ( '' )
    print ( '  Row: ' ),

    for i in range ( i2lo, i2hi + 1 ):
      print ( '%7d       ' % ( i ) ),

    print ( '' )
    print ( '  Col' )

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

    for j in range ( j2lo, j2hi + 1 ):

      print ( '%7d :' % ( j ) ),
      
      for i in range ( i2lo, i2hi + 1 ):
        print ( '%12g  ' % ( a[i,j] ) ),

      print ( '' )

  return

def r8mat_transpose_print_some_test ( ):

#*****************************************************************************80
#
## r8mat_transpose_print_some_test() tests r8mat_transpose_print_some().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 October 2014
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'r8mat_transpose_print_some_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  r8mat_transpose_print_some prints some of an R8MAT, transposed.' )

  m = 4
  n = 6
  v = np.array ( [ \
    [ 11.0, 12.0, 13.0, 14.0, 15.0, 16.0 ], 
    [ 21.0, 22.0, 23.0, 24.0, 25.0, 26.0 ], 
    [ 31.0, 32.0, 33.0, 34.0, 35.0, 36.0 ], 
    [ 41.0, 42.0, 43.0, 44.0, 45.0, 46.0 ] ], dtype = np.float64 )
  r8mat_transpose_print_some ( m, n, v, 0, 3, 2, 5, '  R8MAT, rows 0:2, cols 3:5:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'r8mat_transpose_print_some_test:' )
  print ( '  Normal end of execution.' )
  return

def r8vec_print ( n, a, title ):

#*****************************************************************************80
#
## r8vec_print() prints an R8VEC.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    31 August 2014
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the dimension of the vector.
#
#    real A(N), the vector to be printed.
#
#    string TITLE, a title.
#
  print ( '' )
  print ( title )
  print ( '' )
  for i in range ( 0, n ):
    print ( '%6d:  %12g' % ( i, a[i] ) )

def r8vec_print_test ( ):

#*****************************************************************************80
#
## r8vec_print_test() tests r8vec_print().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    29 October 2014
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'r8vec_print_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  r8vec_print prints an R8VEC.' )

  n = 4
  v = np.array ( [ 123.456, 0.000005, -1.0E+06, 3.14159265 ], dtype = np.float64 )
  r8vec_print ( n, v, '  Here is an R8VEC:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'r8vec_print_test:' )
  print ( '  Normal end of execution.' )
  return

def r8vec_reverse ( n, a1 ):

#*****************************************************************************80
#
## r8vec_reverse() reverses the elements of an R8VEC.
#
#  Example:
#
#    Input:
#
#      N = 5, A = ( 11.0, 12.0, 13.0, 14.0, 15.0 ).
#
#    Output:
#
#      A = ( 15.0, 14.0, 13.0, 12.0, 11.0 ).
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    09 June 2015
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of entries in the array.
#
#    real A1(N), the array to be reversed.
#
#  Output:
#
#    real A2(N), the reversed array.
#
  import numpy as np

  a2 = np.zeros ( n )
  for i in range ( 0, n ):
    a2[i] = a1[n-1-i]

  return a2

def r8vec_reverse_test ( ):

#*****************************************************************************80
#
## r8vec_reverse_test() tests r8vec_reverse().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    30 August 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  n = 5

  print ( '' )
  print ( 'r8vec_reverse_test' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  r8vec_reverse reverses an R8VEC.' )
 
  a1 = np.zeros ( n )
  for i in range ( 0, n ):
    a1[i] = float ( i + 1 )
 
  r8vec_print ( n, a1, '  Original array:' )

  a2 = r8vec_reverse ( n, a1 )

  r8vec_print ( n, a2, '  Reversed array:' )
#
#  Terminate.
#
  print ( '' )
  print ( 'r8vec_reverse_test:' )
  print ( '  Normal end of execution.' )
  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    06 April 2013
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return None

if ( __name__ == '__main__' ):
  timestamp( )
  padua_test ( )
  timestamp ( )