subroutine filename_inc ( filename )

!*****************************************************************************80
!
!! filename_inc() increments a partially numeric filename.
!
!  Discussion:
!
!    It is assumed that the digits in the name, whether scattered or
!    connected, represent a number that is to be increased by 1 on
!    each call.  If this number is all 9's on input, the output number
!    is all 0's.  Non-numeric letters of the name are unaffected.
!
!    If the name is empty, then the routine stops.
!
!    If the name contains no digits, the empty string is returned.
!
!  Example:
!
!      Input            Output
!      -----            ------
!      'a7to11.txt'     'a7to12.txt'
!      'a7to99.txt'     'a8to00.txt'
!      'a9to99.txt'     'a0to00.txt'
!      'cat.txt'        ' '
!      ' '              STOP!
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    19 September 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character ( len = * ) FILENAME.
!    On input, a character string to be incremented.
!    On output, the incremented string.
!
  implicit none

  character c
  integer change
  integer digit
  character ( len = * ) filename
  integer i
  integer lens

  lens = len_trim ( filename )

  if ( lens <= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILENAME_INC - Fatal error!'
    write ( *, '(a)' ) '  The input string is empty.'
    stop 1
  end if

  change = 0

  do i = lens, 1, -1

    c = filename(i:i)

    if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then

      change = change + 1

      digit = ichar ( c ) - 48
      digit = digit + 1

      if ( digit == 10 ) then
        digit = 0
      end if

      c = char ( digit + 48 )

      filename(i:i) = c

      if ( c /= '0' ) then
        return
      end if

    end if

  end do
!
!  No digits were found.  Return blank.
!
  if ( change == 0 ) then
    filename = ' '
    return
  end if

  return
end
subroutine get_unit ( iunit )

!*****************************************************************************80
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!  Discussion:
!
!    A "free" FORTRAN unit number is a value between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5, 6 and 9, which
!    are commonly reserved for console I/O).
!
!    Otherwise, IUNIT is a value between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    26 October 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT, the free unit number.
!
  implicit none

  integer i
  integer ios
  integer iunit
  logical lopen
 
  iunit = 0
 
  do i = 1, 99
 
    if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )
 
      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if
 
  end do

  return
end
subroutine padua_order ( l, n )

!*****************************************************************************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:
!
!    09 February 2014
!
!  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.
!
!  Parameters:
!
!    Input, integer L, the level of the set.
!    0 <= L
!
!    Output, integer N, the order (number of points) in the set.
!
  implicit none

  integer i
  integer l
  integer n

  n = 0
  do i = 0, l
    n = n + ( l / 2 ) + 1
    if ( mod ( l, 2 ) == 1 .and. mod ( i, 2 ) == 1 ) then
      n = n + 1
    end if
  end do

  return
end
subroutine padua_plot ( l, filename )

!*****************************************************************************80
!
!! PADUA_PLOT plots the Padua points of given level.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    29 May 2014
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer L, the level of the set.
!    0 <= L
!
!    Input, character ( len = * ) FILENAME, a common filename prefix for 
!    the files to be created.
!
  implicit none

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

  character ( len = 255 ) command_filename
  integer command_unit
  character ( len = 255 ) data_filename
  integer data_unit
  character ( len = * ) filename
  integer j
  integer l
  integer n
  character ( len = 255 ) plot_filename
  real ( kind = rk ), allocatable :: xy(:,:)

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

  allocate ( xy(1:2,1:n) )

  call padua_points ( l, xy )
!
!  Create graphics data file.
!
  call get_unit ( data_unit )
  data_filename = trim ( filename ) // '_data.txt'
  open ( unit = data_unit, file = data_filename, status = 'replace' )
  do j = 1, n
    write ( data_unit, '(2x,g24.16,2x,g24.16)' ) xy(1:2,j)
  end do
  close ( unit = data_unit )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Created data file "' // trim ( data_filename ) // '".'
!
!  Create graphics command file.
!
  call get_unit ( command_unit )
  command_filename = trim ( filename ) // '_commands.txt'
  open ( unit = command_unit, file = command_filename, status = 'replace' )
  write ( command_unit, '(a)' ) '# ' // trim ( command_filename )
  write ( command_unit, '(a)' ) '#'
  write ( command_unit, '(a)' ) '# Usage:'
  write ( command_unit, '(a)' ) '#  gnuplot < ' // trim ( command_filename )
  write ( command_unit, '(a)' ) '#'
  write ( command_unit, '(a)' ) 'set term png'
  plot_filename = trim ( filename ) // '.png'
  write ( command_unit, '(a)' ) 'set output "' // trim ( plot_filename ) // '"'
  write ( command_unit, '(a)' ) 'set xlabel "<--- X --->"'
  write ( command_unit, '(a)' ) 'set ylabel "<--- Y --->"'
  write ( command_unit, '(a,i2,a)' ) 'set title "Padua Points, Level ', l, '"'
  write ( command_unit, '(a)' ) 'set grid'
  write ( command_unit, '(a)' ) 'set key off'
  write ( command_unit, '(a)' ) 'set size ratio -1'
  write ( command_unit, '(a)' ) 'set style data lines'
  write ( command_unit, '(a)' ) 'set timestamp'
  write ( command_unit, '(a)' ) 'plot [-1:+1] [-1:+1] "' // &
    trim ( data_filename ) // &
    '" using 1:2 with points lt 3 pt 3'
  close ( unit = command_unit )

  write ( *, '(a)' ) &
    '  Created command file "' // trim ( command_filename ) // '".'
!
!  Free memory.
!
  deallocate ( xy )

  return
end
subroutine padua_points ( l, xy )

!*****************************************************************************80
!
!! PADUA_POINTS returns the Padua points of order N.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    09 June 2014
!
!  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.
!
!  Parameters:
!
!    Input, integer L, the level of the set.
!    0 <= L
!
!    Output, real ( kind = rk ) XY(2,((L+1)*(L+2))/2), the Padua points.
!
  implicit none

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

  integer l

  real ( kind = rk ) angle1
  real ( kind = rk ) angle2
  integer i
  integer j
  integer j_hi
  integer k
  real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00
  real ( kind = rk ) xy(2,((l+1)*(l+2))/2)

  if ( l == 0 ) then
    xy(1,1) = 0.0D+00
    xy(2,1) = 0.0D+00
    return
  end if

  k = 0

  do i = 0, l

    j_hi = ( l / 2 ) + 1
    if ( mod ( l, 2 ) == 1 .and. mod ( i, 2 ) == 1 ) then
      j_hi = j_hi + 1
    end if

    do j = 1, j_hi

     k = k + 1

      if ( i * 2 == l ) then
        xy(1,k) = 0.0D+00
      else
        angle1 = real ( i, kind = rk ) * r8_pi / real ( l, kind = rk )
        xy(1,k) = cos ( angle1 )
      end if

      if ( mod ( i, 2 ) == 0 ) then
        if ( 2 * ( 2 * j - 1 ) == l + 1 ) then
          xy(2,k) = 0.0D+00
        else
          angle2 = real ( 2 * j - 1, kind = rk ) * r8_pi &
                 / real ( l + 1, kind = rk )
          xy(2,k) = cos ( angle2 )
        end if
      else
        if ( 2 * ( 2 * j - 2 ) == l + 1 ) then
          xy(2,k) = 0.0D+00
        else
          angle2 = real ( 2 * j - 2, kind = rk ) * r8_pi &
                 / real ( l + 1, kind = rk )
          xy(2,k) = cos ( angle2 )
        end if
      end if

    end do

  end do

  return
end
subroutine padua_points_set ( l, xy )

!*****************************************************************************80
!
!! PADUA_POINTS_SET sets the Padua points.
!
!  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.
!
!  Parameters:
!
!    Input, integer L, the level.
!    0 <= L <= 10.
!
!    Output, real ( kind = rk ) XY(2,N), the Padua points.
!
  implicit none

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

  integer l

  integer j1
  integer j1_hi
  integer j2
  integer n
  real ( kind = rk ) t
  real ( kind = rk ) xy(2,((l+1)*(l+2))/2)

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

  if ( l == 0 ) then

    xy = reshape ( (/ &
      0.000000000000000D+00, 0.000000000000000D+00 /), (/ 2, n /) )
 
  else if ( l == 1 ) then

    xy = reshape ( (/ &
      -1.000000000000000D+00, -1.000000000000000D+00, &  
      -1.000000000000000D+00,  1.000000000000000D+00, &
       1.000000000000000D+00,  0.000000000000000D+00 /), (/ 2, n /) )

  else if ( l == 2 ) then

    xy = reshape ( (/ &
      -1.000000000000000D+00, -1.000000000000000D+00, &
      -1.000000000000000D+00,  0.5000000000000001D+00, &
       0.000000000000000D+00, -0.4999999999999998D+00, &  
       0.000000000000000D+00,  1.000000000000000D+00, &    
       1.000000000000000D+00, -1.000000000000000D+00, &
       1.000000000000000D+00,  0.5000000000000001D+00 /), (/ 2, n /) )
 
  else if ( l == 3 ) then

    xy = reshape ( (/ &
      -1.000000000000000D+00,  -1.000000000000000D+00, &    
      -1.000000000000000D+00,   0.000000000000000D+00, &
      -1.000000000000000D+00,   1.000000000000000D+00, &    
      -0.4999999999999998D+00, -0.7071067811865475D+00, &    
      -0.4999999999999998D+00,  0.7071067811865476D+00, &    
       0.5000000000000001D+00, -1.000000000000000D+00, &    
       0.5000000000000001D+00,  0.000000000000000D+00, &    
       0.5000000000000001D+00,  1.000000000000000D+00, &    
       1.000000000000000D+00,  -0.7071067811865475D+00, &    
       1.000000000000000D+00,   0.7071067811865476D+00 /), (/ 2, n /) )
  
  else if ( l == 4 ) then

    xy = reshape ( (/ &
      -1.000000000000000D+00,  -1.000000000000000D+00, &
      -1.000000000000000D+00,  -0.3090169943749473D+00, &
      -1.000000000000000D+00,   0.8090169943749475D+00, &
      -0.7071067811865475D+00, -0.8090169943749473D+00, &
      -0.7071067811865475D+00,  0.3090169943749475D+00, &
      -0.7071067811865475D+00,  1.000000000000000D+00, &
       0.000000000000000D+00,  -1.000000000000000D+00, &
       0.000000000000000D+00,  -0.3090169943749473D+00, &
       0.000000000000000D+00,   0.8090169943749475D+00, &
       0.7071067811865476D+00, -0.8090169943749473D+00, &
       0.7071067811865476D+00,  0.3090169943749475D+00, &
       0.7071067811865476D+00,  1.000000000000000D+00, &
       1.000000000000000D+00,  -1.000000000000000D+00, &
       1.000000000000000D+00,  -0.3090169943749473D+00, &
       1.000000000000000D+00,   0.8090169943749475D+00 /), (/ 2, n /) )

  else if ( l == 5 ) then

    xy = reshape ( (/ &
      -1.000000000000000D+00,  -1.000000000000000D+00, &
      -1.000000000000000D+00,  -0.4999999999999998D+00, &   
      -1.000000000000000D+00,   0.5000000000000001D+00, &    
      -1.000000000000000D+00,   1.000000000000000D+00, &   
      -0.8090169943749473D+00, -0.8660254037844387D+00, &
      -0.8090169943749473D+00,  0.000000000000000D+00, & 
      -0.8090169943749473D+00,  0.8660254037844387D+00, &
      -0.3090169943749473D+00, -1.000000000000000D+00, & 
      -0.3090169943749473D+00, -0.4999999999999998D+00, &   
      -0.3090169943749473D+00,  0.5000000000000001D+00, &
      -0.3090169943749473D+00,  1.000000000000000D+00, & 
       0.3090169943749475D+00, -0.8660254037844387D+00, &
       0.3090169943749475D+00,  0.000000000000000D+00, &  
       0.3090169943749475D+00,  0.8660254037844387D+00, &   
       0.8090169943749475D+00, -1.000000000000000D+00, &  
       0.8090169943749475D+00, -0.4999999999999998D+00, &  
       0.8090169943749475D+00,  0.5000000000000001D+00, &
       0.8090169943749475D+00,  1.000000000000000D+00, &
       1.000000000000000D+00,  -0.8660254037844387D+00, &   
       1.000000000000000D+00,   0.000000000000000D+00, &
       1.000000000000000D+00,   0.8660254037844387D+00 /), (/ 2, n /) )

  else if ( l == 6 ) then   

    xy = reshape ( (/ &
      -1.000000000000000D+00,  -1.000000000000000D+00, &  
      -1.000000000000000D+00,  -0.6234898018587335D+00, &    
      -1.000000000000000D+00,   0.2225209339563144D+00, &    
      -1.000000000000000D+00,   0.9009688679024191D+00, &    
      -0.8660254037844387D+00, -0.9009688679024190D+00, &    
      -0.8660254037844387D+00, -0.2225209339563143D+00, & 
      -0.8660254037844387D+00,  0.6234898018587336D+00, &
      -0.8660254037844387D+00,  1.000000000000000D+00, &
      -0.4999999999999998D+00, -1.000000000000000D+00, &  
      -0.4999999999999998D+00, -0.6234898018587335D+00, &    
      -0.4999999999999998D+00,  0.2225209339563144D+00, &    
      -0.4999999999999998D+00,  0.9009688679024191D+00, &    
       0.000000000000000D+00,  -0.9009688679024190D+00, & 
       0.000000000000000D+00,  -0.2225209339563143D+00, &
       0.000000000000000D+00,   0.6234898018587336D+00, &    
       0.000000000000000D+00,   1.000000000000000D+00, &
       0.5000000000000001D+00, -1.000000000000000D+00, &   
       0.5000000000000001D+00, -0.6234898018587335D+00, &    
       0.5000000000000001D+00,  0.2225209339563144D+00, &    
       0.5000000000000001D+00,  0.9009688679024191D+00, &    
       0.8660254037844387D+00, -0.9009688679024190D+00, &    
       0.8660254037844387D+00, -0.2225209339563143D+00, & 
       0.8660254037844387D+00,  0.6234898018587336D+00, &
       0.8660254037844387D+00,  1.000000000000000D+00, &    
       1.000000000000000D+00,  -1.000000000000000D+00, &  
       1.000000000000000D+00,  -0.6234898018587335D+00, & 
       1.000000000000000D+00,   0.2225209339563144D+00, &
       1.000000000000000D+00,   0.9009688679024191D+00 /), (/ 2, n /) )

  else if ( l == 7 ) then    

    xy = reshape ( (/ &
      -1.000000000000000D+00, -1.000000000000000D+00, &
      -1.000000000000000D+00, -0.7071067811865475D+00, &    
      -1.000000000000000D+00,  0.000000000000000D+00, &
      -1.000000000000000D+00,  0.7071067811865476D+00, & 
      -1.000000000000000D+00,  1.000000000000000D+00, &  
      -0.9009688679024190D+00, -0.9238795325112867D+00, &    
      -0.9009688679024190D+00, -0.3826834323650897D+00, &    
      -0.9009688679024190D+00,  0.3826834323650898D+00, & 
      -0.9009688679024190D+00,  0.9238795325112867D+00, &
      -0.6234898018587335D+00, -1.000000000000000D+00, &
      -0.6234898018587335D+00, -0.7071067811865475D+00, &
      -0.6234898018587335D+00,  0.000000000000000D+00, &
      -0.6234898018587335D+00,  0.7071067811865476D+00, &
      -0.6234898018587335D+00, 1.000000000000000D+00, & 
      -0.2225209339563143D+00, -0.9238795325112867D+00, &    
      -0.2225209339563143D+00, -0.3826834323650897D+00, &    
      -0.2225209339563143D+00,  0.3826834323650898D+00, & 
      -0.2225209339563143D+00,  0.9238795325112867D+00, &
       0.2225209339563144D+00, -1.000000000000000D+00, &
       0.2225209339563144D+00, -0.7071067811865475D+00, &
       0.2225209339563144D+00,  0.000000000000000D+00, &
       0.2225209339563144D+00,  0.7071067811865476D+00, &
       0.2225209339563144D+00,  1.000000000000000D+00, &  
       0.6234898018587336D+00, -0.9238795325112867D+00, &    
       0.6234898018587336D+00, -0.3826834323650897D+00, &    
       0.6234898018587336D+00,  0.3826834323650898D+00, & 
       0.6234898018587336D+00,  0.9238795325112867D+00, &
       0.9009688679024191D+00, -1.000000000000000D+00, &
       0.9009688679024191D+00, -0.7071067811865475D+00, &
       0.9009688679024191D+00,  0.000000000000000D+00, &
       0.9009688679024191D+00,  0.7071067811865476D+00, &
       0.9009688679024191D+00,  1.000000000000000D+00, &   
       1.000000000000000D+00, -0.9238795325112867D+00, &    
       1.000000000000000D+00, -0.3826834323650897D+00, &    
       1.000000000000000D+00,  0.3826834323650898D+00, &
       1.000000000000000D+00,  0.9238795325112867D+00 /), (/ 2, n /) )

  else if ( l == 8 ) then 

    xy = reshape ( (/ &
      -1.000000000000000D+00, -1.000000000000000D+00, &   
      -1.000000000000000D+00, -0.7660444431189779D+00, &    
      -1.000000000000000D+00, -0.1736481776669303D+00, &    
      -1.000000000000000D+00,  0.5000000000000001D+00, &    
      -1.000000000000000D+00,  0.9396926207859084D+00, &    
      -0.9238795325112867D+00, -0.9396926207859083D+00, &    
      -0.9238795325112867D+00, -0.4999999999999998D+00, &    
      -0.9238795325112867D+00,  0.1736481776669304D+00, &    
      -0.9238795325112867D+00,  0.7660444431189780D+00, &
      -0.9238795325112867D+00,  1.000000000000000D+00, &
      -0.7071067811865475D+00, -1.000000000000000D+00, &   
      -0.7071067811865475D+00, -0.7660444431189779D+00, &    
      -0.7071067811865475D+00, -0.1736481776669303D+00, &    
      -0.7071067811865475D+00,  0.5000000000000001D+00, &    
      -0.7071067811865475D+00,  0.9396926207859084D+00, &    
      -0.3826834323650897D+00, -0.9396926207859083D+00, &    
      -0.3826834323650897D+00, -0.4999999999999998D+00, &    
      -0.3826834323650897D+00,  0.1736481776669304D+00, &    
      -0.3826834323650897D+00,  0.7660444431189780D+00, &
      -0.3826834323650897D+00,  1.000000000000000D+00, &    
       0.000000000000000D+00, -1.000000000000000D+00, &   
       0.000000000000000D+00, -0.7660444431189779D+00, &    
       0.000000000000000D+00, -0.1736481776669303D+00, &    
       0.000000000000000D+00,  0.5000000000000001D+00, &    
       0.000000000000000D+00,  0.9396926207859084D+00, &    
       0.3826834323650898D+00, -0.9396926207859083D+00, &    
       0.3826834323650898D+00, -0.4999999999999998D+00, &    
       0.3826834323650898D+00,  0.1736481776669304D+00, &    
       0.3826834323650898D+00,  0.7660444431189780D+00, &
       0.3826834323650898D+00,  1.000000000000000D+00, &
       0.7071067811865476D+00, -1.000000000000000D+00, &   
       0.7071067811865476D+00, -0.7660444431189779D+00, &    
       0.7071067811865476D+00, -0.1736481776669303D+00, &    
       0.7071067811865476D+00,  0.5000000000000001D+00, &    
       0.7071067811865476D+00,  0.9396926207859084D+00, &    
       0.9238795325112867D+00, -0.9396926207859083D+00, &    
       0.9238795325112867D+00, -0.4999999999999998D+00, &    
       0.9238795325112867D+00,  0.1736481776669304D+00, &    
       0.9238795325112867D+00,  0.7660444431189780D+00, &
       0.9238795325112867D+00,  1.000000000000000D+00, &    
       1.000000000000000D+00, -1.000000000000000D+00, &   
       1.000000000000000D+00, -0.7660444431189779D+00, &    
       1.000000000000000D+00, -0.1736481776669303D+00, &  
       1.000000000000000D+00,  0.5000000000000001D+00, &
       1.000000000000000D+00,  0.9396926207859084D+00 /), (/ 2, n /) )

  else if ( l == 9 ) then 

    xy = reshape ( (/ &
      -1.000000000000000D+00, -1.000000000000000D+00, &   
      -1.000000000000000D+00, -0.8090169943749473D+00, &    
      -1.000000000000000D+00, -0.3090169943749473D+00, &    
      -1.000000000000000D+00,  0.3090169943749475D+00, &
      -1.000000000000000D+00,  0.8090169943749475D+00, &    
      -1.000000000000000D+00,  1.000000000000000D+00, &   
      -0.9396926207859083D+00, -0.9510565162951535D+00, &    
      -0.9396926207859083D+00, -0.5877852522924730D+00, &
      -0.9396926207859083D+00,  0.000000000000000D+00, &   
      -0.9396926207859083D+00,  0.5877852522924731D+00, &    
      -0.9396926207859083D+00,  0.9510565162951535D+00, &
      -0.7660444431189779D+00, -1.000000000000000D+00, &   
      -0.7660444431189779D+00, -0.8090169943749473D+00, &    
      -0.7660444431189779D+00, -0.3090169943749473D+00, &    
      -0.7660444431189779D+00,  0.3090169943749475D+00, &    
      -0.7660444431189779D+00,  0.8090169943749475D+00, & 
      -0.7660444431189779D+00,  1.000000000000000D+00, &   
      -0.4999999999999998D+00, -0.9510565162951535D+00, &    
      -0.4999999999999998D+00, -0.5877852522924730D+00, &
      -0.4999999999999998D+00,  0.000000000000000D+00, &   
      -0.4999999999999998D+00,  0.5877852522924731D+00, &    
      -0.4999999999999998D+00,  0.9510565162951535D+00, &
      -0.1736481776669303D+00, -1.000000000000000D+00, &   
      -0.1736481776669303D+00, -0.8090169943749473D+00, &    
      -0.1736481776669303D+00, -0.3090169943749473D+00, &    
      -0.1736481776669303D+00,  0.3090169943749475D+00, &    
      -0.1736481776669303D+00,  0.8090169943749475D+00, &
      -0.1736481776669303D+00,  1.000000000000000D+00, &   
       0.1736481776669304D+00, -0.9510565162951535D+00, &    
       0.1736481776669304D+00, -0.5877852522924730D+00, &
       0.1736481776669304D+00,  0.000000000000000D+00, &   
       0.1736481776669304D+00,  0.5877852522924731D+00, &    
       0.1736481776669304D+00,  0.9510565162951535D+00, &
       0.5000000000000001D+00, -1.000000000000000D+00, &   
       0.5000000000000001D+00, -0.8090169943749473D+00, &    
       0.5000000000000001D+00, -0.3090169943749473D+00, &    
       0.5000000000000001D+00,  0.3090169943749475D+00, &    
       0.5000000000000001D+00,  0.8090169943749475D+00, &
       0.5000000000000001D+00,  1.000000000000000D+00, &   
       0.7660444431189780D+00, -0.9510565162951535D+00, &    
       0.7660444431189780D+00, -0.5877852522924730D+00, &
       0.7660444431189780D+00,  0.000000000000000D+00, &   
       0.7660444431189780D+00,  0.5877852522924731D+00, &    
       0.7660444431189780D+00,  0.9510565162951535D+00, &
       0.9396926207859084D+00, -1.000000000000000D+00, &   
       0.9396926207859084D+00, -0.8090169943749473D+00, &    
       0.9396926207859084D+00, -0.3090169943749473D+00, &    
       0.9396926207859084D+00,  0.3090169943749475D+00, &    
       0.9396926207859084D+00,  0.8090169943749475D+00, &
       0.9396926207859084D+00,  1.000000000000000D+00, &   
       1.000000000000000D+00, -0.9510565162951535D+00, &
       1.000000000000000D+00, -0.5877852522924730D+00, &    
       1.000000000000000D+00,  0.000000000000000D+00, &
       1.000000000000000D+00,  0.5877852522924731D+00, &
       1.000000000000000D+00,  0.9510565162951535D+00 /), (/ 2, n /) ) 

  else if ( l == 10 ) then

    xy = reshape ( (/ &
      -1.000000000000000D+00, -1.000000000000000D+00, &   
      -1.000000000000000D+00, -0.8412535328311811D+00, &    
      -1.000000000000000D+00, -0.4154150130018863D+00, &    
      -1.000000000000000D+00,  0.1423148382732851D+00, &    
      -1.000000000000000D+00,  0.6548607339452851D+00, &    
      -1.000000000000000D+00,  0.9594929736144974D+00, &    
      -0.9510565162951535D+00, -0.9594929736144974D+00, &    
      -0.9510565162951535D+00, -0.6548607339452850D+00, &    
      -0.9510565162951535D+00, -0.1423148382732850D+00, &    
      -0.9510565162951535D+00,  0.4154150130018864D+00, &    
      -0.9510565162951535D+00,  0.8412535328311812D+00, &
      -0.9510565162951535D+00,  1.000000000000000D+00, &
      -0.8090169943749473D+00, -1.000000000000000D+00, &   
      -0.8090169943749473D+00, -0.8412535328311811D+00, &    
      -0.8090169943749473D+00, -0.4154150130018863D+00, &    
      -0.8090169943749473D+00,  0.1423148382732851D+00, &    
      -0.8090169943749473D+00,  0.6548607339452851D+00, &    
      -0.8090169943749473D+00,  0.9594929736144974D+00, &    
      -0.5877852522924730D+00, -0.9594929736144974D+00, &    
      -0.5877852522924730D+00, -0.6548607339452850D+00, &    
      -0.5877852522924730D+00, -0.1423148382732850D+00, &    
      -0.5877852522924730D+00,  0.4154150130018864D+00, &    
      -0.5877852522924730D+00,  0.8412535328311812D+00, &
      -0.5877852522924730D+00,  1.000000000000000D+00, &
      -0.3090169943749473D+00, -1.000000000000000D+00, &   
      -0.3090169943749473D+00, -0.8412535328311811D+00, &    
      -0.3090169943749473D+00, -0.4154150130018863D+00, &    
      -0.3090169943749473D+00,  0.1423148382732851D+00, &    
      -0.3090169943749473D+00,  0.6548607339452851D+00, &   
      -0.3090169943749473D+00, 0.9594929736144974D+00, &    
       0.000000000000000D+00, -0.9594929736144974D+00, &    
       0.000000000000000D+00, -0.6548607339452850D+00, &    
       0.000000000000000D+00, -0.1423148382732850D+00, &    
       0.000000000000000D+00,  0.4154150130018864D+00, &
       0.000000000000000D+00,  0.8412535328311812D+00, &    
       0.000000000000000D+00,  1.000000000000000D+00, & 
       0.3090169943749475D+00, -1.000000000000000D+00, &    
       0.3090169943749475D+00, -0.8412535328311811D+00, &    
       0.3090169943749475D+00, -0.4154150130018863D+00, &    
       0.3090169943749475D+00,  0.1423148382732851D+00, &    
       0.3090169943749475D+00,  0.6548607339452851D+00, &    
       0.3090169943749475D+00,  0.9594929736144974D+00, &    
       0.5877852522924731D+00, -0.9594929736144974D+00, &    
       0.5877852522924731D+00, -0.6548607339452850D+00, &    
       0.5877852522924731D+00, -0.1423148382732850D+00, &    
       0.5877852522924731D+00,  0.4154150130018864D+00, &    
       0.5877852522924731D+00,  0.8412535328311812D+00, &
       0.5877852522924731D+00, 1.000000000000000D+00, &
       0.8090169943749475D+00, -1.000000000000000D+00, &   
       0.8090169943749475D+00, -0.8412535328311811D+00, &    
       0.8090169943749475D+00, -0.4154150130018863D+00, &    
       0.8090169943749475D+00,  0.1423148382732851D+00, &    
       0.8090169943749475D+00,  0.6548607339452851D+00, &    
       0.8090169943749475D+00,  0.9594929736144974D+00, &    
       0.9510565162951535D+00, -0.9594929736144974D+00, &    
       0.9510565162951535D+00, -0.6548607339452850D+00, &    
       0.9510565162951535D+00, -0.1423148382732850D+00, &    
       0.9510565162951535D+00,  0.4154150130018864D+00, &    
       0.9510565162951535D+00,  0.8412535328311812D+00, &
       0.9510565162951535D+00,  1.000000000000000D+00, &    
       1.000000000000000D+00, -1.000000000000000D+00, &   
       1.000000000000000D+00, -0.8412535328311811D+00, &    
       1.000000000000000D+00, -0.4154150130018863D+00, &    
       1.000000000000000D+00,  0.1423148382732851D+00, &    
       1.000000000000000D+00,  0.6548607339452851D+00, &
       1.000000000000000D+00,  0.9594929736144974D+00 /), (/ 2, n /) ) 

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PADUA_POINTS_SET - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal value of L = ', l
    write ( *, '(a)' ) '  Legal values are 0 through 10.'
    stop 1

  end if
!
!  Reverse the order to match the published data.
!
  j1_hi = n / 2

  do j1 = 1, j1_hi
    j2 = n + 1 - j1
    t        = xy(1,j1)
    xy(1,j1) = xy(1,j2)
    xy(1,j2) = t
    t        = xy(2,j1)
    xy(2,j1) = xy(2,j2)
    xy(2,j2) = t
  end do

  return
end
subroutine padua_weights_set ( l, w )

!*****************************************************************************80
!
!! PADUA_WEIGHTS_SET sets quadrature weights for the Padua points.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    09 June 2014
!
!  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.
!
!  Parameters:
!
!    Input, integer L, the level.
!    0 <= L <= 10.
!
!    Output, real ( kind = rk ) W(N), the quadrature weights.
!
  implicit none

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

  integer l

  integer n
  real ( kind = rk ) w(((l+1)*(l+2))/2)

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

  return
end
subroutine padua_weights ( l, w )

!*****************************************************************************80
!
!! PADUA_WEIGHTS returns quadrature weights do Padua points.
!
!  Discussion:
!
!    The order of the weights corresponds to the ordering used
!    by the companion function padua_points().
!
!    Caliari, de Marchi and Vianello supplied a MATLAB code pdwtsMM
!    which carries out this same computation in a way that makes
!    more efficient use of MATLAB's vector and matrix capabilities.  
!    This version of the computation was painfully rewritten to display 
!    the individual scalar computations, so that it could be translated
!    into other languages.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    10 June 2014
!
!  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.
!
!  Parameters:
!
!    Input, integer L, the level of the set.
!    0 <= L
!
!    Output, real ( kind = rk ) W((L+1)*(L+2)/2), the quadrature weights.
!
  implicit none

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

  integer l

  real ( kind = rk ) angle
  integer i
  integer i2
  integer j
  integer j2
  integer lp1h
  integer lp2h
  integer lp3h
  real ( kind = rk ) mi
  real ( kind = rk ) mj
  real ( kind = rk ), allocatable :: mom(:,:)
  real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00
  real ( kind = rk ), allocatable :: te1(:,:)
  real ( kind = rk ), allocatable :: te2(:,:)
  real ( kind = rk ), allocatable :: tmteo(:,:)
  real ( kind = rk ), allocatable :: tmtoe(:,:)
  real ( kind = rk ), allocatable :: to1(:,:)
  real ( kind = rk ), allocatable :: to2(:,:)
  real ( kind = rk ) w(((l+1)*(l+2))/2)
  real ( kind = rk ), allocatable :: w1(:,:)
  real ( kind = rk ), allocatable :: w2(:,:)

  if ( l == 0 ) then
    w(1) = 4.0D+00
    return
  end if
!
!  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.
!
  allocate ( te1(1:lp2h,1:lp2h) )
  do j = 1, lp2h
    do i = 1, lp2h
      angle = r8_pi * real ( 2 * ( i - 1 ) * ( 2 * j - 2 ), kind = rk ) &
        / real ( l, kind = rk )
      te1(i,j) = cos ( angle )
    end do
  end do
  te1(2:lp2h,1:lp2h) = te1(2:lp2h,1:lp2h) * sqrt ( 2.0D+00 )

  allocate ( to1(1:lp2h,1:lp1h) )
  do j = 1, lp1h
    do i = 1, lp2h
      angle = r8_pi * real ( 2 * ( i - 1 ) * ( 2 * j - 1 ), kind = rk ) &
        / real ( l, kind = rk )
      to1(i,j) = cos ( angle )
    end do
  end do
  to1(2:lp2h,1:lp1h) = to1(2:lp2h,1:lp1h) * sqrt ( 2.0D+00 )

  allocate ( te2(1:lp2h,1:lp3h) )
  do j = 1, lp3h
    do i = 1, lp2h
      angle = r8_pi * real ( 2 * ( i - 1 ) * ( 2 * j - 2 ), kind = rk ) &
        / real ( l + 1, kind = rk )
      te2(i,j) = cos ( angle )
    end do
  end do
  te2(2:lp2h,1:lp3h) = te2(2:lp2h,1:lp3h) * sqrt ( 2.0D+00 )

  allocate ( to2(1:lp2h,1:lp2h) )
  do j = 1, lp2h
    do i = 1, lp2h
      angle = r8_pi * real ( 2 * ( i - 1 ) * ( 2 * j - 1 ), kind = rk ) &
        / real ( l + 1, kind = rk )
      to2(i,j) = cos ( angle )
    end do
  end do
  to2(2:lp2h,1:lp2h) = to2(2:lp2h,1:lp2h) * sqrt ( 2.0D+00 )
!
!  MOM: Moments matrix do even * even pairs.
!
  allocate ( mom(1:lp2h,1:lp2h) )

  do j = 1, lp2h
    mj = 2.0D+00 * sqrt ( 2.0D+00 ) / real ( 1 - ( 2 * j - 2 ) ** 2, kind = rk )
    do i = 1, lp2h + 1 - j
      mi = 2.0D+00 * sqrt ( 2.0D+00 ) &
        / real ( 1 - ( 2 * i - 2 ) ** 2, kind = rk )
      mom(i,j) = mi * mj
    end do
  end do

  mom(1,1:lp2h) = mom(1,1:lp2h) / sqrt ( 2.0D+00 )
  mom(1:lp2h,1) = mom(1:lp2h,1) / sqrt ( 2.0D+00 )

  if ( mod ( l, 2 ) == 0 ) then
    mom(lp2h,1) = mom(lp2h,1) / 2.0D+00
  end if
!
!  TMTOE and TMTEO: matrix products.
!
  allocate ( tmtoe(1:lp2h,1:lp2h) )
  tmtoe(1:lp2h,1:lp2h) = 0.0D+00

  do j2 = 1, lp2h
    do i2 = 1, lp2h + 1 - j2
      do j = 1, lp2h
        do i = 1, lp2h
          tmtoe(i,j) = tmtoe(i,j) + to2(i2,i) * mom(j2,i2) * te1(j2,j)
        end do
      end do
    end do
  end do

  allocate ( tmteo(1:lp3h,1:lp1h) )
  tmteo(1:lp3h,1:lp1h) = 0.0D+00

  do j2 = 1, lp2h
    do i2 = 1, lp2h + 1 - j2
      do j = 1, lp1h
        do i = 1, lp3h
          tmteo(i,j) = tmteo(i,j) + te2(i2,i) * mom(j2,i2) * to1(j2,j)
        end do
      end do
    end do
  end do
!
!  W1 and W2: Interpolation weight matrices.
!
  allocate ( w1(1:lp2h,1:lp2h) )
  w1(1:lp2h,1:lp2h) = 2.0D+00 / real ( l * ( l + 1 ), kind = rk )

  w1(1:lp2h,1) = w1(1:lp2h,1) / 2.0D+00

  if ( mod ( l, 2 ) == 0 ) then
    w1(1:lp2h,lp2h) = w1(1:lp2h,lp2h) / 2.0D+00
    w1(lp2h,1:lp2h) = w1(lp2h,1:lp2h) / 2.0D+00
  end if

  allocate ( w2(1:lp3h,1:lp1h) )
  w2(1:lp3h,1:lp1h) = 2.0D+00  / real ( l * ( l + 1 ), kind = rk )

  w2(1,1:lp1h) = w2(1,1:lp1h) / 2.0D+00

  if ( mod ( l, 2 ) == 1 ) then
    w2(lp3h,1:lp1h) = w2(lp3h,1:lp1h) / 2.0D+00
    w2(1:lp3h,lp1h) = w2(1:lp3h,lp1h) / 2.0D+00
  end if
!
!  Cubature weights as matrices on the subgrids.
!
  do j = 1, lp2h
    do i = 1, lp2h
      w1(i,j) = w1(i,j) * tmtoe(i,j)
    end do
  end do

  do j = 1, lp1h
    do i = 1, lp3h
      w2(i,j) = w2(i,j) * tmteo(i,j)
    end do
  end do
!
!  Pack weight matrices W1 and W2 into the vector W.
!
  if ( mod ( l, 2 ) == 0 ) then

    do j = 1, lp2h
      do i = 1, lp2h
        w(i+(2*j-2)*lp2h) = w1(i,j)
      end do
    end do

    do j = 1, lp1h
      do i = 1, lp3h
        w(i+(2*j-1)*lp2h) = w2(i,j)
      end do
    end do

  else

    do j = 1, lp1h
      do i = 1, lp2h
        w(i+(j-1)*(l+2)) = w1(i,j)
      end do
    end do

    do j = 1, lp1h
      do i = 1, lp3h
        w(i+lp2h+(j-1)*(l+2)) = w2(i,j)
      end do
    end do

  end if
!
!  Free memory.
!
  deallocate ( te1 )
  deallocate ( te2 )
  deallocate ( tmteo )
  deallocate ( tmtoe )
  deallocate ( to1 )
  deallocate ( to2 )
  deallocate ( w1 )
  deallocate ( w2 )

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

!*****************************************************************************80
!
!! R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed.
!
!  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:
!
!    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 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 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 )

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

  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,a,5a14)' ) j, ':', ( ctemp(i), i = 1, inc )

    end do

  end do

  return
end
subroutine r8vec_print ( n, a, title )

!*****************************************************************************80
!
!! R8VEC_PRINT prints an R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    22 August 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components of the vector.
!
!    Input, real ( kind = rk ) A(N), the vector to be printed.
!
!    Input, character ( len = * ) TITLE, a title.
!
  implicit none

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

  integer n

  real ( kind = rk ) a(n)
  integer i
  character ( len = * ) title

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

  do i = 1, n
    write ( *, '(2x,i8,a,1x,g16.8)' ) i, ':', a(i)
  end do

  return
end
subroutine r8vec_reverse ( n, a )

!*****************************************************************************80
!
!! R8VEC_REVERSE reverses the elements of an R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!    In FORTRAN90, calling R8VEC_REVERSE is equivalent to
!
!      A(1:N) = A(N:1:-1)
!
!  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:
!
!    30 September 2009
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the array.
!
!    Input/output, real ( kind = rk ) A(N), the array to be reversed.
!
  implicit none

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

  integer n

  real ( kind = rk ) a(n)

  a(1:n) = a(n:1:-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
!
  implicit none

  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.2,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