program main

!*****************************************************************************80
!
!! MAIN is the main program for SQUARE_EXACTNESS_TEST.
!
!  Discussion:
!
!    SQUARE_EXACTNESS_TEST tests the SQUARE_EXACTNESS library.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    29 May 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'SQUARE_EXACTNESS_TEST'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Test the SQUARE_EXACTNESS library.'

  call test01 ( )
  call test02 ( )
!
!  Terminate.
!
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'SQUARE_EXACTNESS_TEST'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ''
  call timestamp ( )

  stop
end
subroutine test01 ( )

!*****************************************************************************80
!
!! TEST01 tests product Gauss-Legendre rules for the Legendre 2D integral.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    29 May 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real ( kind = 8 ) a(2)
  real ( kind = 8 ) b(2)
  integer ( kind = 4 ) l
  integer ( kind = 4 ) n
  integer ( kind = 4 ) n_1d
  integer ( kind = 4 ) p_max
  integer ( kind = 4 ) t
  real ( kind = 8 ), allocatable :: w(:)
  real ( kind = 8 ), allocatable :: x(:)
  real ( kind = 8 ), allocatable :: y(:)

  a(1:2) = -1.0D+00
  b(1:2) = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  Product Gauss-Legendre rules for the 2D Legendre integral.'
  write ( *, '(a)' ) '  Density function rho(x) = 1.'
  write ( *, '(a)' ) '  Region: -1 <= x <= +1.'
  write ( *, '(a)' ) '  Region: -1 <= y <= +1.'
  write ( *, '(a)' ) '  Level: L'
  write ( *, '(a)' ) '  Exactness: 2*L+1'
  write ( *, '(a)' ) '  Order: N = (L+1)*(L+1)'

  do l = 0, 5

    n_1d = l + 1
    n = n_1d * n_1d
    t = 2 * l + 1

    allocate ( x(1:n) )
    allocate ( y(1:n) )
    allocate ( w(1:n) )

    call legendre_2d_set ( a, b, n_1d, n_1d, x, y, w )

    p_max = t + 1
    call legendre_2d_exactness ( a, b, n, x, y, w, p_max )

    deallocate ( x )
    deallocate ( y )
    deallocate ( w )

  end do

  return
end
subroutine test02 ( )

!*****************************************************************************80
!
!! TEST02 tests Padua rules for the Legendre 2D integral.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    30 May 2014
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real ( kind = 8 ) a(2)
  real ( kind = 8 ) b(2)
  integer ( kind = 4 ) l
  integer ( kind = 4 ) n
  integer ( kind = 4 ) p_max
  real ( kind = 8 ), allocatable :: w(:)
  real ( kind = 8 ), allocatable :: x(:)
  real ( kind = 8 ), allocatable :: y(:)

  a(1:2) = -1.0D+00
  b(1:2) = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  Padua rule for the 2D Legendre integral.'
  write ( *, '(a)' ) '  Density function rho(x) = 1.'
  write ( *, '(a)' ) '  Region: -1 <= x <= +1.'
  write ( *, '(a)' ) '  Region: -1 <= y <= +1.'
  write ( *, '(a)' ) '  Level: L'
  write ( *, '(a)' ) '  Exactness: L+1 when L is 0,'
  write ( *, '(a)' ) '             L   otherwise.'
  write ( *, '(a)' ) '  Order: N = (L+1)*(L+2)/2'

  do l = 0, 5

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

    allocate ( w(1:n) )
    allocate ( x(1:n) )
    allocate ( y(1:n) )

    call padua_point_set ( l, x, y )
    call padua_weight_set ( l, w )
 
    if ( l == 0 ) then
      p_max = l + 2
    else
      p_max = l + 1
    end if

    call legendre_2d_exactness ( a, b, n, x, y, w, p_max )

    deallocate ( w )
    deallocate ( x )
    deallocate ( y )

  end do

  return
end
subroutine legendre_2d_set ( a, b, nx, ny, x, y, w )

!*****************************************************************************80
!
!! LEGENDRE_2D_SET: set a 2D Gauss-Legendre quadrature rule.
!
!  Discussion:
!
!    The integral:
!
!      integral ( a(2) <= y <= b(2) ) ( a(1) <= x <= b(1) ) f(x,y) dx dy
!
!    The quadrature rule:
!
!      sum ( 1 <= i <= n ) w(i) * f ( x(i),y(i) )
!
!    where n = nx * ny, the orders of the rule in the X and Y directions.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    29 May 2014
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A(2), B(2), the lower and upper integration
!    limits.
!
!    Input, integer ( kind = 4 ) NX, NY, the orders in the X and Y directions.
!    NX and NY must be between 1 and 10.
!
!    Output, real ( kind = 8 ) X(N), Y(N), the abscissas.
!
!    Output, real ( kind = 8 ) W(N), the weights.
!
  implicit none

  integer ( kind = 4 ) nx
  integer ( kind = 4 ) ny

  real ( kind = 8 ) a(2)
  real ( kind = 8 ) b(2)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  real ( kind = 8 ) w(nx*ny)
  real ( kind = 8 ) wx(nx)
  real ( kind = 8 ) wy(ny)
  real ( kind = 8 ) x(nx*ny)
  real ( kind = 8 ) xx(nx)
  real ( kind = 8 ) y(nx*ny)
  real ( kind = 8 ) yy(ny)
!
!  Get the rules for [-1,+1].
!
  call legendre_set ( nx, xx, wx )
  call legendre_set ( ny, yy, wy )
!
!  Adjust from [-1,+1] to [A,B].
!
  do i = 1, nx
    xx(i) = ( ( 1.0D+00 - xx(i) ) * a(1)   &
            + ( 1.0D+00 + xx(i) ) * b(1) ) &
            /   2.0D+00
    wx(i) = wx(i) * ( b(1) - a(1) ) / 2.0D+00
  end do

  do j = 1, ny
    yy(j) = ( ( 1.0D+00 - yy(j) ) * a(2)   &
            + ( 1.0D+00 + yy(j) ) * b(2) ) &
            /   2.0D+00
    wy(j) = wy(j) * ( b(2) - a(2) ) / 2.0D+00
  end do
!
!  Compute the product rule.
!
  k = 0
  do j = 1, ny
    do i = 1, nx
      k = k + 1
      x(k) = xx(i)
      y(k) = yy(j)
      w(k) = wx(i) * wy(j)
    end do
  end do

  return
end
subroutine legendre_set ( n, x, w )

!*****************************************************************************80
!
!! LEGENDRE_SET sets abscissas and weights for Gauss-Legendre quadrature.
!
!  Discussion:
!
!    The integral:
!
!      integral ( -1 <= x <= 1 ) f(x) dx
!
!    The quadrature rule:
!
!      sum ( 1 <= i <= n ) w(i) * f ( x(i) )
!
!    The quadrature rule is exact for polynomials through degree 2*N-1.
!
!    The abscissas are the zeroes of the Legendre polynomial P(N)(X).
!
!    Mathematica can compute the abscissas and weights of a Gauss-Legendre
!    rule of order N for the interval [A,B] with P digits of precision
!    by the commands:
!
!    Needs["NumericalDifferentialEquationAnalysis`"]
!    GaussianQuadratureWeights [ n, a, b, p ]
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    20 April 2010
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    National Bureau of Standards, 1964,
!    ISBN: 0-486-61272-4,
!    LC: QA47.A34.
!
!    Vladimir Krylov,
!    Approximate Calculation of Integrals,
!    Dover, 2006,
!    ISBN: 0486445798,
!    LC: QA311.K713.
!
!    Arthur Stroud, Don Secrest,
!    Gaussian Quadrature Formulas,
!    Prentice Hall, 1966,
!    LC: QA299.4G3S7.
!
!    Stephen Wolfram,
!    The Mathematica Book,
!    Fourth Edition,
!    Cambridge University Press, 1999,
!    ISBN: 0-521-64314-7,
!    LC: QA76.95.W65.
!
!    Daniel Zwillinger, editor,
!    CRC Standard Mathematical Tables and Formulae,
!    30th Edition,
!    CRC Press, 1996,
!    ISBN: 0-8493-2479-3,
!    LC: QA47.M315.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the order.
!    N must be between 1 and 10.
!
!    Output, real ( kind = 8 ) X(N), the abscissas.
!
!    Output, real ( kind = 8 ) W(N), the weights.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) w(n)
  real ( kind = 8 ) x(n)

  if ( n == 1 ) then

    x(1) = 0.000000000000000000000000000000D+00

    w(1) = 2.000000000000000000000000000000D+00

  else if ( n == 2 ) then

    x(1) = -0.577350269189625764509148780502D+00
    x(2) = 0.577350269189625764509148780502D+00

    w(1) = 1.000000000000000000000000000000D+00
    w(2) = 1.000000000000000000000000000000D+00

  else if ( n == 3 ) then

    x(1) = -0.774596669241483377035853079956D+00
    x(2) = 0.000000000000000000000000000000D+00
    x(3) = 0.774596669241483377035853079956D+00

    w(1) = 0.555555555555555555555555555556D+00
    w(2) = 0.888888888888888888888888888889D+00
    w(3) = 0.555555555555555555555555555556D+00

  else if ( n == 4 ) then

    x(1) = -0.861136311594052575223946488893D+00
    x(2) = -0.339981043584856264802665759103D+00
    x(3) = 0.339981043584856264802665759103D+00
    x(4) = 0.861136311594052575223946488893D+00

    w(1) = 0.347854845137453857373063949222D+00
    w(2) = 0.652145154862546142626936050778D+00
    w(3) = 0.652145154862546142626936050778D+00
    w(4) = 0.347854845137453857373063949222D+00

  else if ( n == 5 ) then

    x(1) = -0.906179845938663992797626878299D+00
    x(2) = -0.538469310105683091036314420700D+00
    x(3) = 0.000000000000000000000000000000D+00
    x(4) = 0.538469310105683091036314420700D+00
    x(5) = 0.906179845938663992797626878299D+00

    w(1) = 0.236926885056189087514264040720D+00
    w(2) = 0.478628670499366468041291514836D+00
    w(3) = 0.568888888888888888888888888889D+00
    w(4) = 0.478628670499366468041291514836D+00
    w(5) = 0.236926885056189087514264040720D+00

  else if ( n == 6 ) then

    x(1) = -0.932469514203152027812301554494D+00
    x(2) = -0.661209386466264513661399595020D+00
    x(3) = -0.238619186083196908630501721681D+00
    x(4) = 0.238619186083196908630501721681D+00
    x(5) = 0.661209386466264513661399595020D+00
    x(6) = 0.932469514203152027812301554494D+00

    w(1) = 0.171324492379170345040296142173D+00
    w(2) = 0.360761573048138607569833513838D+00
    w(3) = 0.467913934572691047389870343990D+00
    w(4) = 0.467913934572691047389870343990D+00
    w(5) = 0.360761573048138607569833513838D+00
    w(6) = 0.171324492379170345040296142173D+00

  else if ( n == 7 ) then

    x(1) = -0.949107912342758524526189684048D+00
    x(2) = -0.741531185599394439863864773281D+00
    x(3) = -0.405845151377397166906606412077D+00
    x(4) = 0.000000000000000000000000000000D+00
    x(5) = 0.405845151377397166906606412077D+00
    x(6) = 0.741531185599394439863864773281D+00
    x(7) = 0.949107912342758524526189684048D+00

    w(1) = 0.129484966168869693270611432679D+00
    w(2) = 0.279705391489276667901467771424D+00
    w(3) = 0.381830050505118944950369775489D+00
    w(4) = 0.417959183673469387755102040816D+00
    w(5) = 0.381830050505118944950369775489D+00
    w(6) = 0.279705391489276667901467771424D+00
    w(7) = 0.129484966168869693270611432679D+00

  else if ( n == 8 ) then

    x(1) = -0.960289856497536231683560868569D+00
    x(2) = -0.796666477413626739591553936476D+00
    x(3) = -0.525532409916328985817739049189D+00
    x(4) = -0.183434642495649804939476142360D+00
    x(5) = 0.183434642495649804939476142360D+00
    x(6) = 0.525532409916328985817739049189D+00
    x(7) = 0.796666477413626739591553936476D+00
    x(8) = 0.960289856497536231683560868569D+00

    w(1) = 0.101228536290376259152531354310D+00
    w(2) = 0.222381034453374470544355994426D+00
    w(3) = 0.313706645877887287337962201987D+00
    w(4) = 0.362683783378361982965150449277D+00
    w(5) = 0.362683783378361982965150449277D+00
    w(6) = 0.313706645877887287337962201987D+00
    w(7) = 0.222381034453374470544355994426D+00
    w(8) = 0.101228536290376259152531354310D+00

  else if ( n == 9 ) then

    x(1) = -0.968160239507626089835576203D+00
    x(2) = -0.836031107326635794299429788D+00
    x(3) = -0.613371432700590397308702039D+00
    x(4) = -0.324253423403808929038538015D+00
    x(5) = 0.000000000000000000000000000D+00
    x(6) = 0.324253423403808929038538015D+00
    x(7) = 0.613371432700590397308702039D+00
    x(8) = 0.836031107326635794299429788D+00
    x(9) = 0.968160239507626089835576203D+00

    w(1) = 0.081274388361574411971892158111D+00
    w(2) = 0.18064816069485740405847203124D+00
    w(3) = 0.26061069640293546231874286942D+00
    w(4) = 0.31234707704000284006863040658D+00
    w(5) = 0.33023935500125976316452506929D+00
    w(6) = 0.31234707704000284006863040658D+00
    w(7) = 0.26061069640293546231874286942D+00
    w(8) = 0.18064816069485740405847203124D+00
    w(9) = 0.081274388361574411971892158111D+00

  else if ( n == 10 ) then

    x(1) = -0.973906528517171720077964012D+00
    x(2) = -0.865063366688984510732096688D+00
    x(3) = -0.679409568299024406234327365D+00
    x(4) = -0.433395394129247190799265943D+00
    x(5) = -0.148874338981631210884826001D+00
    x(6) = 0.148874338981631210884826001D+00
    x(7) = 0.433395394129247190799265943D+00
    x(8) = 0.679409568299024406234327365D+00
    x(9) = 0.865063366688984510732096688D+00
    x(10) = 0.973906528517171720077964012D+00

    w(1) = 0.066671344308688137593568809893D+00
    w(2) = 0.14945134915058059314577633966D+00
    w(3) = 0.21908636251598204399553493423D+00
    w(4) = 0.26926671930999635509122692157D+00
    w(5) = 0.29552422471475287017389299465D+00
    w(6) = 0.29552422471475287017389299465D+00
    w(7) = 0.26926671930999635509122692157D+00
    w(8) = 0.21908636251598204399553493423D+00
    w(9) = 0.14945134915058059314577633966D+00
    w(10) = 0.066671344308688137593568809893D+00

  else

    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'LEGENDRE_SET - Fatal error!'
    write ( *, '(a)' ) '  Illegal value of N.'
    stop 1

  end if

  return
end
subroutine padua_point_set ( l, x, y )

!*****************************************************************************80
!
!! PADUA_POINT_SET sets the Padua points.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    29 May 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 ( kind = 4 ) L, the level.
!    0 <= L <= 10.
!
!    Output, real ( kind = 8 ) X(N), Y(N), the Padua points.
!
  implicit none

  integer ( kind = 4 ) l

  real ( kind = 8 ) x(((l+1)*(l+2))/2)
  real ( kind = 8 ) y(((l+1)*(l+2))/2)

  if ( l == 0 ) then
    x( 1) =  0.000000000000000D+00
    y( 1) =  0.000000000000000D+00    
  else if ( l == 1 ) then
    x( 1) = -1.000000000000000D+00
    y( 1) = -1.000000000000000D+00    
    x( 2) = -1.000000000000000D+00
    y( 2) =  1.000000000000000D+00    
    x( 3) =  1.000000000000000D+00
    y( 3) =  0.000000000000000D+00    
  else if ( l == 2 ) then
    x( 1) = -1.000000000000000D+00
    y( 1) = -1.000000000000000D+00    
    x( 2) = -1.000000000000000D+00
    y( 2) =  0.5000000000000001D+00
    x( 3) =  0.000000000000000D+00
    y( 3) = -0.4999999999999998D+00    
    x( 4) =  0.000000000000000D+00
    y( 4) =  1.000000000000000D+00    
    x( 5) =  1.000000000000000D+00
    y( 5) = -1.000000000000000D+00
    x( 6) =  1.000000000000000D+00
    y( 6) =  0.5000000000000001D+00    
  else if ( l == 3 ) then
    x( 1) = -1.000000000000000D+00
    y( 1) = -1.000000000000000D+00    
    x( 2) = -1.000000000000000D+00
    y( 2) =  0.000000000000000D+00
    x( 3) = -1.000000000000000D+00
    y( 3) =  1.000000000000000D+00    
    x( 4) = -0.4999999999999998D+00
    y( 4) = -0.7071067811865475D+00    
    x( 5) = -0.4999999999999998D+00
    y( 5) =  0.7071067811865476D+00    
    x( 6) =  0.5000000000000001D+00
    y( 6) = -1.000000000000000D+00    
    x( 7) =  0.5000000000000001D+00
    y( 7) =  0.000000000000000D+00    
    x( 8) =  0.5000000000000001D+00
    y( 8) =  1.000000000000000D+00    
    x( 9) =  1.000000000000000D+00
    y( 9) = -0.7071067811865475D+00    
    x(10) =  1.000000000000000D+00
    y(10) =  0.7071067811865476D+00    
  else if ( l == 4 ) then
    x( 1) = -1.000000000000000D+00
    y( 1) = -1.000000000000000D+00   
    x( 2) = -1.000000000000000D+00
    y( 2) = -0.3090169943749473D+00    
    x( 3) = -1.000000000000000D+00
    y( 3) =  0.8090169943749475D+00    
    x( 4) = -0.7071067811865475D+00
    y( 4) = -0.8090169943749473D+00    
    x( 5) = -0.7071067811865475D+00
    y( 5) =  0.3090169943749475D+00
    x( 6) = -0.7071067811865475D+00
    y( 6) =  1.000000000000000D+00    
    x( 7) =  0.000000000000000D+00
    y( 7) = -1.000000000000000D+00   
    x( 8) =  0.000000000000000D+00
    y( 8) = -0.3090169943749473D+00    
    x( 9) =  0.000000000000000D+00
    y( 9) =  0.8090169943749475D+00    
    x(10) =  0.7071067811865476D+00
    y(10) = -0.8090169943749473D+00    
    x(11) =  0.7071067811865476D+00
    y(11) =  0.3090169943749475D+00 
    x(12) =  0.7071067811865476D+00
    y(12) =  1.000000000000000D+00    
    x(13) =  1.000000000000000D+00
    y(13) = -1.000000000000000D+00   
    x(14) =  1.000000000000000D+00
    y(14) = -0.3090169943749473D+00    
    x(15) =  1.000000000000000D+00
    y(15) =  0.8090169943749475D+00    
  else if ( l == 5 ) then
    x( 1) = -1.000000000000000D+00
    y( 1) = -1.000000000000000D+00
    x( 2) = -1.000000000000000D+00
    y( 2) = -0.4999999999999998D+00    
    x( 3) = -1.000000000000000D+00
    y( 3) =  0.5000000000000001D+00    
    x( 4) = -1.000000000000000D+00
    y( 4) =  1.000000000000000D+00    
    x( 5) = -0.8090169943749473D+00
    y( 5) = -0.8660254037844387D+00
    x( 6) = -0.8090169943749473D+00
    y( 6) =  0.000000000000000D+00   
    x( 7) = -0.8090169943749473D+00
    y( 7) =  0.8660254037844387D+00
    x( 8) = -0.3090169943749473D+00
    y( 8) = -1.000000000000000D+00   
    x( 9) = -0.3090169943749473D+00
    y( 9) = -0.4999999999999998D+00    
    x(10) = -0.3090169943749473D+00
    y(10) =  0.5000000000000001D+00
    x(11) = -0.3090169943749473D+00
    y(11) =  1.000000000000000D+00   
    x(12) =  0.3090169943749475D+00
    y(12) = -0.8660254037844387D+00
    x(13) =  0.3090169943749475D+00
    y(13) =  0.000000000000000D+00    
    x(14) =  0.3090169943749475D+00
    y(14) =  0.8660254037844387D+00    
    x(15) =  0.8090169943749475D+00
    y(15) = -1.000000000000000D+00    
    x(16) =  0.8090169943749475D+00
    y(16) = -0.4999999999999998D+00    
    x(17) =  0.8090169943749475D+00
    y(17) =  0.5000000000000001D+00
    x(18) =  0.8090169943749475D+00
    y(18) =  1.000000000000000D+00
    x(19) =  1.000000000000000D+00
    y(19) = -0.8660254037844387D+00    
    x(20) =  1.000000000000000D+00
    y(20) =  0.000000000000000D+00
    x(21) =  1.000000000000000D+00
    y(21) =  0.8660254037844387D+00 
  else if ( l == 6 ) then   
    x( 1) = -1.000000000000000D+00
    y( 1) = -1.000000000000000D+00   
    x( 2) = -1.000000000000000D+00
    y( 2) = -0.6234898018587335D+00    
    x( 3) = -1.000000000000000D+00
    y( 3) =  0.2225209339563144D+00    
    x( 4) = -1.000000000000000D+00
    y( 4) =  0.9009688679024191D+00    
    x( 5) = -0.8660254037844387D+00
    y( 5) = -0.9009688679024190D+00    
    x( 6) = -0.8660254037844387D+00
    y( 6) = -0.2225209339563143D+00    
    x( 7) = -0.8660254037844387D+00
    y( 7) =  0.6234898018587336D+00
    x( 8) = -0.8660254037844387D+00
    y( 8) =  1.000000000000000D+00
    x( 9) = -0.4999999999999998D+00
    y( 9) = -1.000000000000000D+00   
    x(10) = -0.4999999999999998D+00
    y(10) = -0.6234898018587335D+00    
    x(11) = -0.4999999999999998D+00
    y(11) =  0.2225209339563144D+00    
    x(12) = -0.4999999999999998D+00
    y(12) =  0.9009688679024191D+00    
    x(13) =  0.000000000000000D+00
    y(13) = -0.9009688679024190D+00    
    x(14) =  0.000000000000000D+00
    y(14) = -0.2225209339563143D+00
    x(15) =  0.000000000000000D+00
    y(15) =  0.6234898018587336D+00    
    x(16) =  0.000000000000000D+00
    y(16) =  1.000000000000000D+00
    x(17) =  0.5000000000000001D+00
    y(17) = -1.000000000000000D+00    
    x(18) =  0.5000000000000001D+00
    y(18) = -0.6234898018587335D+00    
    x(19) =  0.5000000000000001D+00
    y(19) =  0.2225209339563144D+00    
    x(20) =  0.5000000000000001D+00
    y(20) =  0.9009688679024191D+00    
    x(21) =  0.8660254037844387D+00
    y(21) = -0.9009688679024190D+00    
    x(22) =  0.8660254037844387D+00
    y(22) = -0.2225209339563143D+00    
    x(23) =  0.8660254037844387D+00
    y(23) =  0.6234898018587336D+00
    x(24) =  0.8660254037844387D+00
    y(24) =  1.000000000000000D+00    
    x(25) =  1.000000000000000D+00
    y(25) = -1.000000000000000D+00   
    x(26) =  1.000000000000000D+00
    y(26) = -0.6234898018587335D+00    
    x(27) =  1.000000000000000D+00
    y(27) =  0.2225209339563144D+00
    x(28) =  1.000000000000000D+00
    y(28) =  0.9009688679024191D+00
  else if ( l == 7 ) then    
    x( 1) = -1.000000000000000D+00
    y( 1) = -1.000000000000000D+00
    x( 2) = -1.000000000000000D+00
    y( 2) = -0.7071067811865475D+00    
    x( 3) = -1.000000000000000D+00
    y( 3) =  0.000000000000000D+00   
    x( 4) = -1.000000000000000D+00
    y( 4) =  0.7071067811865476D+00 
    x( 5) = -1.000000000000000D+00
    y( 5) =  1.000000000000000D+00    
    x( 6) = -0.9009688679024190D+00
    y( 6) = -0.9238795325112867D+00    
    x( 7) = -0.9009688679024190D+00
    y( 7) = -0.3826834323650897D+00    
    x( 8) = -0.9009688679024190D+00
    y( 8) =  0.3826834323650898D+00    
    x( 9) = -0.9009688679024190D+00
    y( 9) =  0.9238795325112867D+00
    x(10) = -0.6234898018587335D+00
    y(10) = -1.000000000000000D+00   
    x(11) = -0.6234898018587335D+00
    y(11) = -0.7071067811865475D+00
    x(12) = -0.6234898018587335D+00
    y(12) =  0.000000000000000D+00   
    x(13) = -0.6234898018587335D+00
    y(13) =  0.7071067811865476D+00
    x(14) = -0.6234898018587335D+00
    y(14) =  1.000000000000000D+00   
    x(15) = -0.2225209339563143D+00
    y(15) = -0.9238795325112867D+00    
    x(16) = -0.2225209339563143D+00
    y(16) = -0.3826834323650897D+00    
    x(17) = -0.2225209339563143D+00
    y(17) =  0.3826834323650898D+00    
    x(18) = -0.2225209339563143D+00
    y(18) =  0.9238795325112867D+00
    x(19) =  0.2225209339563144D+00
    y(19) = -1.000000000000000D+00   
    x(20) =  0.2225209339563144D+00
    y(20) = -0.7071067811865475D+00
    x(21) =  0.2225209339563144D+00
    y(21) =  0.000000000000000D+00
    x(22) =  0.2225209339563144D+00
    y(22) =  0.7071067811865476D+00
    x(23) =  0.2225209339563144D+00
    y(23) =  1.000000000000000D+00   
    x(24) =  0.6234898018587336D+00
    y(24) = -0.9238795325112867D+00    
    x(25) =  0.6234898018587336D+00
    y(25) = -0.3826834323650897D+00    
    x(26) =  0.6234898018587336D+00
    y(26) =  0.3826834323650898D+00    
    x(27) =  0.6234898018587336D+00
    y(27) =  0.9238795325112867D+00
    x(28) =  0.9009688679024191D+00
    y(28) = -1.000000000000000D+00   
    x(29) =  0.9009688679024191D+00
    y(29) = -0.7071067811865475D+00
    x(30) =  0.9009688679024191D+00
    y(30) =  0.000000000000000D+00   
    x(31) =  0.9009688679024191D+00
    y(31) =  0.7071067811865476D+00
    x(32) =  0.9009688679024191D+00
    y(32) =  1.000000000000000D+00   
    x(33) =  1.000000000000000D+00
    y(33) = -0.9238795325112867D+00    
    x(34) =  1.000000000000000D+00
    y(34) = -0.3826834323650897D+00    
    x(35) =  1.000000000000000D+00
    y(35) =  0.3826834323650898D+00    
    x(36) =  1.000000000000000D+00
    y(36) =  0.9238795325112867D+00 
  else if ( l == 8 ) then 
    x( 1) = -1.000000000000000D+00
    y( 1) = -1.000000000000000D+00   
    x( 2) = -1.000000000000000D+00
    y( 2) = -0.7660444431189779D+00    
    x( 3) = -1.000000000000000D+00
    y( 3) = -0.1736481776669303D+00    
    x( 4) = -1.000000000000000D+00
    y( 4) =  0.5000000000000001D+00    
    x( 5) = -1.000000000000000D+00
    y( 5) =  0.9396926207859084D+00    
    x( 6) = -0.9238795325112867D+00
    y( 6) = -0.9396926207859083D+00    
    x( 7) = -0.9238795325112867D+00
    y( 7) = -0.4999999999999998D+00    
    x( 8) = -0.9238795325112867D+00
    y( 8) =  0.1736481776669304D+00    
    x( 9) = -0.9238795325112867D+00
    y( 9) =  0.7660444431189780D+00
    x(10) = -0.9238795325112867D+00
    y(10) =  1.000000000000000D+00
    x(11) = -0.7071067811865475D+00
    y(11) = -1.000000000000000D+00   
    x(12) = -0.7071067811865475D+00
    y(12) = -0.7660444431189779D+00    
    x(13) = -0.7071067811865475D+00
    y(13) = -0.1736481776669303D+00    
    x(14) = -0.7071067811865475D+00
    y(14) =  0.5000000000000001D+00    
    x(15) = -0.7071067811865475D+00
    y(15) =  0.9396926207859084D+00    
    x(16) = -0.3826834323650897D+00
    y(16) = -0.9396926207859083D+00    
    x(17) = -0.3826834323650897D+00
    y(17) = -0.4999999999999998D+00    
    x(18) = -0.3826834323650897D+00
    y(18) =  0.1736481776669304D+00    
    x(19) = -0.3826834323650897D+00
    y(19) =  0.7660444431189780D+00
    x(20) = -0.3826834323650897D+00
    y(20) =  1.000000000000000D+00    
    x(21) =  0.000000000000000D+00
    y(21) = -1.000000000000000D+00   
    x(22) =  0.000000000000000D+00
    y(22) = -0.7660444431189779D+00    
    x(23) =  0.000000000000000D+00
    y(23) = -0.1736481776669303D+00    
    x(24) =  0.000000000000000D+00
    y(24) =  0.5000000000000001D+00    
    x(25) =  0.000000000000000D+00
    y(25) =  0.9396926207859084D+00    
    x(26) =  0.3826834323650898D+00
    y(26) = -0.9396926207859083D+00    
    x(27) =  0.3826834323650898D+00
    y(27) = -0.4999999999999998D+00    
    x(28) =  0.3826834323650898D+00
    y(28) =  0.1736481776669304D+00    
    x(29) =  0.3826834323650898D+00
    y(29) =  0.7660444431189780D+00
    x(30) =  0.3826834323650898D+00
    y(30) =  1.000000000000000D+00
    x(31) =  0.7071067811865476D+00
    y(31) = -1.000000000000000D+00   
    x(32) =  0.7071067811865476D+00
    y(32) = -0.7660444431189779D+00    
    x(33) =  0.7071067811865476D+00
    y(33) = -0.1736481776669303D+00    
    x(34) =  0.7071067811865476D+00
    y(34) =  0.5000000000000001D+00    
    x(35) =  0.7071067811865476D+00
    y(35) =  0.9396926207859084D+00    
    x(36) =  0.9238795325112867D+00
    y(36) = -0.9396926207859083D+00    
    x(37) =  0.9238795325112867D+00
    y(37) = -0.4999999999999998D+00    
    x(38) =  0.9238795325112867D+00
    y(38) =  0.1736481776669304D+00    
    x(39) =  0.9238795325112867D+00
    y(39) =  0.7660444431189780D+00
    x(40) =  0.9238795325112867D+00
    y(40) =  1.000000000000000D+00    
    x(41) =  1.000000000000000D+00
    y(41) = -1.000000000000000D+00   
    x(42) =  1.000000000000000D+00
    y(42) = -0.7660444431189779D+00    
    x(43) =  1.000000000000000D+00
    y(43) = -0.1736481776669303D+00    
    x(44) =  1.000000000000000D+00
    y(44) =  0.5000000000000001D+00
    x(45) =  1.000000000000000D+00
    y(45) =  0.9396926207859084D+00   
  else if ( l == 9 ) then 
    x( 1) = -1.000000000000000D+00
    y( 1) = -1.000000000000000D+00   
    x( 2) = -1.000000000000000D+00
    y( 2) = -0.8090169943749473D+00    
    x( 3) = -1.000000000000000D+00
    y( 3) = -0.3090169943749473D+00    
    x( 4) = -1.000000000000000D+00
    y( 4) =  0.3090169943749475D+00
    x( 5) = -1.000000000000000D+00
    y( 5) =  0.8090169943749475D+00    
    x( 6) = -1.000000000000000D+00
    y( 6) =  1.000000000000000D+00   
    x( 7) = -0.9396926207859083D+00
    y( 7) = -0.9510565162951535D+00    
    x( 8) = -0.9396926207859083D+00
    y( 8) = -0.5877852522924730D+00
    x( 9) = -0.9396926207859083D+00
    y( 9) =  0.000000000000000D+00   
    x(10) = -0.9396926207859083D+00
    y(10) =  0.5877852522924731D+00    
    x(11) = -0.9396926207859083D+00
    y(11) =  0.9510565162951535D+00
    x(12) = -0.7660444431189779D+00
    y(12) = -1.000000000000000D+00   
    x(13) = -0.7660444431189779D+00
    y(13) = -0.8090169943749473D+00    
    x(14) = -0.7660444431189779D+00
    y(14) = -0.3090169943749473D+00    
    x(15) = -0.7660444431189779D+00
    y(15) =  0.3090169943749475D+00    
    x(16) = -0.7660444431189779D+00
    y(16) =  0.8090169943749475D+00 
    x(17) = -0.7660444431189779D+00
    y(17) =  1.000000000000000D+00   
    x(18) = -0.4999999999999998D+00
    y(18) = -0.9510565162951535D+00    
    x(19) = -0.4999999999999998D+00
    y(19) = -0.5877852522924730D+00
    x(20) = -0.4999999999999998D+00
    y(20) =  0.000000000000000D+00   
    x(21) = -0.4999999999999998D+00
    y(21) =  0.5877852522924731D+00    
    x(22) = -0.4999999999999998D+00
    y(22) =  0.9510565162951535D+00
    x(23) = -0.1736481776669303D+00
    y(23) = -1.000000000000000D+00   
    x(24) = -0.1736481776669303D+00
    y(24) = -0.8090169943749473D+00    
    x(25) = -0.1736481776669303D+00
    y(25) = -0.3090169943749473D+00    
    x(26) = -0.1736481776669303D+00
    y(26) =  0.3090169943749475D+00    
    x(27) = -0.1736481776669303D+00
    y(27) =  0.8090169943749475D+00
    x(28) = -0.1736481776669303D+00
    y(28) =  1.000000000000000D+00   
    x(29) =  0.1736481776669304D+00
    y(29) = -0.9510565162951535D+00    
    x(30) =  0.1736481776669304D+00
    y(30) = -0.5877852522924730D+00
    x(31) =  0.1736481776669304D+00
    y(31) =  0.000000000000000D+00   
    x(32) =  0.1736481776669304D+00
    y(32) =  0.5877852522924731D+00    
    x(33) =  0.1736481776669304D+00
    y(33) =  0.9510565162951535D+00
    x(34) =  0.5000000000000001D+00
    y(34) = -1.000000000000000D+00   
    x(35) =  0.5000000000000001D+00
    y(35) = -0.8090169943749473D+00    
    x(36) =  0.5000000000000001D+00
    y(36) = -0.3090169943749473D+00    
    x(37) =  0.5000000000000001D+00
    y(37) =  0.3090169943749475D+00    
    x(38) =  0.5000000000000001D+00
    y(38) =  0.8090169943749475D+00
    x(39) =  0.5000000000000001D+00
    y(39) =  1.000000000000000D+00   
    x(40) =  0.7660444431189780D+00
    y(40) = -0.9510565162951535D+00    
    x(41) =  0.7660444431189780D+00
    y(41) = -0.5877852522924730D+00
    x(42) =  0.7660444431189780D+00
    y(42) =  0.000000000000000D+00   
    x(43) =  0.7660444431189780D+00
    y(43) =  0.5877852522924731D+00    
    x(44) =  0.7660444431189780D+00
    y(44) =  0.9510565162951535D+00
    x(45) =  0.9396926207859084D+00
    y(45) = -1.000000000000000D+00   
    x(46) =  0.9396926207859084D+00
    y(46) = -0.8090169943749473D+00    
    x(47) =  0.9396926207859084D+00
    y(47) = -0.3090169943749473D+00    
    x(48) =  0.9396926207859084D+00
    y(48) =  0.3090169943749475D+00    
    x(49) =  0.9396926207859084D+00
    y(49) =  0.8090169943749475D+00
    x(50) =  0.9396926207859084D+00
    y(50) =  1.000000000000000D+00   
    x(51) =  1.000000000000000D+00
    y(51) = -0.9510565162951535D+00
    x(52) =  1.000000000000000D+00
    y(52) = -0.5877852522924730D+00    
    x(53) =  1.000000000000000D+00
    y(53) =  0.000000000000000D+00   
    x(54) =  1.000000000000000D+00
    y(54) =  0.5877852522924731D+00
    x(55) =  1.000000000000000D+00
    y(55) =  0.9510565162951535D+00    
  else if ( l == 10 ) then
    x( 1) = -1.000000000000000D+00
    y( 1) = -1.000000000000000D+00   
    x( 2) = -1.000000000000000D+00
    y( 2) = -0.8412535328311811D+00    
    x( 3) = -1.000000000000000D+00
    y( 3) = -0.4154150130018863D+00    
    x( 4) = -1.000000000000000D+00
    y( 4) =  0.1423148382732851D+00    
    x( 5) = -1.000000000000000D+00
    y( 5) =  0.6548607339452851D+00    
    x( 6) = -1.000000000000000D+00
    y( 6) =  0.9594929736144974D+00    
    x( 7) = -0.9510565162951535D+00
    y( 7) = -0.9594929736144974D+00    
    x( 8) = -0.9510565162951535D+00
    y( 8) = -0.6548607339452850D+00    
    x( 9) = -0.9510565162951535D+00
    y( 9) = -0.1423148382732850D+00    
    x(10) = -0.9510565162951535D+00
    y(10) =  0.4154150130018864D+00    
    x(11) = -0.9510565162951535D+00
    y(11) =  0.8412535328311812D+00
    x(12) = -0.9510565162951535D+00
    y(12) =  1.000000000000000D+00
    x(13) = -0.8090169943749473D+00
    y(13) = -1.000000000000000D+00   
    x(14) = -0.8090169943749473D+00
    y(14) = -0.8412535328311811D+00    
    x(15) = -0.8090169943749473D+00
    y(15) = -0.4154150130018863D+00    
    x(16) = -0.8090169943749473D+00
    y(16) =  0.1423148382732851D+00    
    x(17) = -0.8090169943749473D+00
    y(17) =  0.6548607339452851D+00    
    x(18) = -0.8090169943749473D+00
    y(18) =  0.9594929736144974D+00    
    x(19) = -0.5877852522924730D+00
    y(19) = -0.9594929736144974D+00    
    x(20) = -0.5877852522924730D+00
    y(20) = -0.6548607339452850D+00    
    x(21) = -0.5877852522924730D+00
    y(21) = -0.1423148382732850D+00    
    x(22) = -0.5877852522924730D+00
    y(22) =  0.4154150130018864D+00    
    x(23) = -0.5877852522924730D+00
    y(23) =  0.8412535328311812D+00
    x(24) = -0.5877852522924730D+00
    y(24) =  1.000000000000000D+00
    x(25) = -0.3090169943749473D+00
    y(25) = -1.000000000000000D+00   
    x(26) = -0.3090169943749473D+00
    y(26) = -0.8412535328311811D+00    
    x(27) = -0.3090169943749473D+00
    y(27) = -0.4154150130018863D+00    
    x(28) = -0.3090169943749473D+00
    y(28) =  0.1423148382732851D+00    
    x(29) = -0.3090169943749473D+00
    y(29) =  0.6548607339452851D+00    
    x(30) = -0.3090169943749473D+00
    y(30) =  0.9594929736144974D+00    
    x(31) =  0.000000000000000D+00
    y(31) = -0.9594929736144974D+00    
    x(32) =  0.000000000000000D+00
    y(32) = -0.6548607339452850D+00    
    x(33) =  0.000000000000000D+00
    y(33) = -0.1423148382732850D+00    
    x(34) =  0.000000000000000D+00
    y(34) =  0.4154150130018864D+00
    x(35) =  0.000000000000000D+00
    y(35) =  0.8412535328311812D+00    
    x(36) =  0.000000000000000D+00
    y(36) =  1.000000000000000D+00 
    x(37) =  0.3090169943749475D+00
    y(37) = -1.000000000000000D+00    
    x(38) =  0.3090169943749475D+00
    y(38) = -0.8412535328311811D+00    
    x(39) =  0.3090169943749475D+00
    y(39) = -0.4154150130018863D+00    
    x(40) =  0.3090169943749475D+00
    y(40) =  0.1423148382732851D+00    
    x(41) =  0.3090169943749475D+00
    y(41) =  0.6548607339452851D+00    
    x(42) =  0.3090169943749475D+00
    y(42) =  0.9594929736144974D+00    
    x(43) =  0.5877852522924731D+00
    y(43) = -0.9594929736144974D+00    
    x(44) =  0.5877852522924731D+00
    y(44) = -0.6548607339452850D+00    
    x(45) =  0.5877852522924731D+00
    y(45) = -0.1423148382732850D+00    
    x(46) =  0.5877852522924731D+00
    y(46) =  0.4154150130018864D+00    
    x(47) =  0.5877852522924731D+00
    y(47) =  0.8412535328311812D+00
    x(48) =  0.5877852522924731D+00
    y(48) =  1.000000000000000D+00
    x(49) =  0.8090169943749475D+00
    y(49) = -1.000000000000000D+00   
    x(50) =  0.8090169943749475D+00
    y(50) = -0.8412535328311811D+00    
    x(51) =  0.8090169943749475D+00
    y(51) = -0.4154150130018863D+00    
    x(52) =  0.8090169943749475D+00
    y(52) =  0.1423148382732851D+00    
    x(53) =  0.8090169943749475D+00
    y(53) =  0.6548607339452851D+00    
    x(54) =  0.8090169943749475D+00
    y(54) =  0.9594929736144974D+00    
    x(55) =  0.9510565162951535D+00
    y(55) = -0.9594929736144974D+00    
    x(56) =  0.9510565162951535D+00
    y(56) = -0.6548607339452850D+00    
    x(57) =  0.9510565162951535D+00
    y(57) = -0.1423148382732850D+00    
    x(58) =  0.9510565162951535D+00
    y(58) =  0.4154150130018864D+00    
    x(59) =  0.9510565162951535D+00
    y(59) =  0.8412535328311812D+00
    x(60) =  0.9510565162951535D+00
    y(60) =  1.000000000000000D+00    
    x(61) =  1.000000000000000D+00
    y(61) = -1.000000000000000D+00   
    x(62) =  1.000000000000000D+00
    y(62) = -0.8412535328311811D+00    
    x(63) =  1.000000000000000D+00
    y(63) = -0.4154150130018863D+00    
    x(64) =  1.000000000000000D+00
    y(64) =  0.1423148382732851D+00    
    x(65) =  1.000000000000000D+00
    y(65) =  0.6548607339452851D+00 
    x(66) =  1.000000000000000D+00
    y(66) =  0.9594929736144974D+00
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PADUA_POINT_SET - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal value of L = ', l
    write ( *, '(a)' ) '  Legal values are 0 through 10.'
    stop 1
  end if

  return
end
subroutine padua_weight_set ( l, w )

!*****************************************************************************80
!
!! PADUA_WEIGHT_SET sets quadrature weights for the Padua points.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    29 May 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 ( kind = 4 ) L, the level.
!    0 <= L <= 10.
!
!    Output, real ( kind = 8 ) W(N), the quadrature weights.
!
  implicit none

  integer ( kind = 4 ) l

  real ( kind = 8 ) 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_WEIGHT_SET - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal value of L = ', l
    write ( *, '(a)' ) '  Legal values are 0 through 10.'
    stop 1
  end if

  return
end