# include <cstdlib>
# include <iostream>
# include <fstream>
# include <sstream>
# include <iomanip>
# include <cmath>
# include <cstring>
# include <ctime>

using namespace std;

# include "padua.hpp"

//****************************************************************************80

void filename_inc ( string *filename )

//****************************************************************************80
//
//  Purpose:
//
//    FILENAME_INC increments a partially numeric file name.
//
//  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"  (typical case.  Last digit incremented)
//      "a7to99.txt"     "a8to00.txt"  (last digit incremented, with carry.)
//      "a9to99.txt"     "a0to00.txt"  (wrap around)
//      "cat.txt"        " "           (no digits to increment)
//      " "              STOP!         (error)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 November 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input/output, string *FILENAME, the filename to be incremented.
//
{
  char c;
  int change;
  int i;
  int lens;

  lens = (*filename).length ( );

  if ( lens <= 0 )
  {
    cerr << "\n";
    cerr << "FILENAME_INC - Fatal error!\n";
    cerr << "  The input string is empty.\n";
    exit ( 1 );
  }

  change = 0;

  for ( i = lens - 1; 0 <= i; i-- )
  {
    c = (*filename)[i];

    if ( '0' <= c && c <= '9' )
    {
      change = change + 1;

      if ( c == '9' )
      {
        c = '0';
        (*filename)[i] = c;
      }
      else
      {
        c = c + 1;
        (*filename)[i] = c;
        return;
      }
    }
  }
//
//  No digits were found.  Return blank.
//
  if ( change == 0 )
  {
    for ( i = lens - 1; 0 <= i; i-- )
    {
      (*filename)[i] = ' ';
    }
  }

  return;
}
//****************************************************************************80

string i4_to_string ( int i4 )

//****************************************************************************80
//
//  Purpose:
//
//    I4_TO_STRING converts an I4 to a C++ string.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    16 January 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I4, an integer.
//
//    Input, string FORMAT, the format string.
//
//    Output, string I4_TO_STRING, the string.
//
{
  ostringstream fred;
  string value;

  fred << i4;

  value = fred.str ( );

  return value;
}
//****************************************************************************80

int padua_order ( int l )

//****************************************************************************80
//
//  Purpose:
//
//    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 GNU LGPL license.
//
//  Modified:
//
//    07 April 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, int L, the level of the set.
//    0 <= L
//
//    Output, int PADUA_ORDER, the order (number of points) in the set.
//
{
  int i;
  int n;

  n = 0;
  for ( i = 0; i <= l; i++ )
  {
    n = n + ( l / 2 ) + 1;
    if ( ( l % 2 ) == 1 && ( i % 2 ) == 1 )
    {
      n = n + 1;
    }
  }

  return n;
}
//****************************************************************************80

void padua_plot ( int l, string filename )

//****************************************************************************80
//
//  Purpose:
//
//    PADUA_PLOT plots the Padua points of given level.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 June 2014
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int L, the level of the set.
//    0 <= L
//
//    Input, string FILENAME, a common filename prefix for 
//    the files to be created.
//
{
  string command_filename;
  ofstream command_unit;
  string data_filename;
  ofstream data_unit;
  int j;
  int n;
  string plot_filename;
  double *xy;

  n = padua_order ( l );

  xy = padua_points ( l );
//
//  Create graphics data file.
//
  data_filename = filename + "_data.txt";
  data_unit.open ( data_filename.c_str ( ) );
  for ( j = 0; j < n; j++ )
  {
    data_unit << "  " << setprecision(16) << xy[0+j*2]
              << "  " << setprecision(16) << xy[1+j*2] << "\n";
  }
  data_unit.close ( );
  cout << "\n";
  cout << "  Created data file '" << data_filename << "'.\n";
//
//  Create graphics command file.
//
  command_filename = filename + "_commands.txt";
  command_unit.open ( command_filename.c_str ( ) );
  command_unit << "# " << command_filename << "\n";
  command_unit << "#\n";
  command_unit << "# Usage:\n";
  command_unit << "#  gnuplot < " << command_filename << "\n";
  command_unit << "#\n";
  command_unit << "set term png\n";
  plot_filename = filename + ".png";
  command_unit << "set output '" << plot_filename << "'\n";
  command_unit << "set xlabel '<--- X --->'\n";
  command_unit << "set ylabel '<--- Y --->'\n";
  command_unit << "set title 'Padua Points, Level " << l << "\n";
  command_unit << "set grid\n";
  command_unit << "set key off\n";
  command_unit << "set size ratio -1\n";
  command_unit << "set style data lines\n";
  command_unit << "set timestamp\n";
  command_unit << "plot [-1:+1] [-1:+1] '" << data_filename 
               << "' using 1:2 with points lt 3 pt 3\n";
  command_unit.close ( );

  cout << "  Created command file '" << command_filename << "'.\n";
//
//  Free memory.
//
  delete [] xy;

  return;
}
//****************************************************************************80

void padua_points_set ( int l, double x[], double y[] )

//****************************************************************************80
//
//  Purpose:
//
//    PADUA_POINTS_SET sets the Padua points.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL 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, int L, the level.
//    0 <= L <= 10.
//
//    Output, double X[N], Y[N], the Padua points.
//
{
  int n;

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

  return;
}
//****************************************************************************80

double *padua_points ( int l )

//****************************************************************************80
//
//  Purpose:
//
//    PADUA_POINTS returns the Padua points of level L.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL 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, int L, the level of the set.
//    0 <= L
//
//    Output, double PADUA_POINTS[2*((L+1)*(L+2))/2)], the Padua points.
//
{
  double angle1;
  double angle2;
  int i;
  int j;
  int j_hi;
  int k;
  int n;
  const double r8_pi = 3.141592653589793;
  double *xy;

  n = ( ( l + 1 ) * ( l + 2 ) ) / 2;
  xy = new double[2*n];

  if ( l == 0 )
  {
    xy[0+0*2] = 0.0;
    xy[1+0*2] = 0.0;
    return xy;
  }

  k = 0;

  for ( i = 0; i <= l; i++ )
  {
    j_hi = ( l / 2 ) + 1;
    if ( ( l % 2 ) == 1 && ( i % 2 ) == 1 )
    {
      j_hi = j_hi + 1;
    }

    for ( j = 1; j <= j_hi; j++ )
    {
      if ( i * 2 == l )
      {
        xy[0+k*2] = 0.0;
      }
      else
      {
        angle1 = ( double ) ( i ) * r8_pi / ( double ) ( l );
        xy[0+k*2] = cos ( angle1 );
      }

      if ( ( i % 2 ) == 0 )
      {
        if ( 2 * ( 2 * j - 1 ) == l + 1 )
        {
          xy[1+k*2] = 0.0;
        }
        else
        {
          angle2 = ( double ) ( 2 * j - 1 ) * r8_pi / ( double ) ( l + 1 );
          xy[1+k*2] = cos ( angle2 );
        }
      }
      else
      {
        if ( 2 * ( 2 * j - 2 ) == l + 1 )
        {
          xy[1+k*2] = 0.0;
        }
        else
        {
          angle2 = ( double ) ( 2 * j - 2 ) * r8_pi / ( double ) ( l + 1 );
          xy[1+k*2] = cos ( angle2 );
        }
      }
      k = k + 1;
    }
  }

  return xy;
}
//****************************************************************************80

double *padua_weights_set ( int l )

//****************************************************************************80
//
//  Purpose:
//
//    PADUA_WEIGHTS_SET sets quadrature weights for the Padua points.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 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, int L, the level.
//    0 <= L <= 10.
//
//    Output, double PADUA_WEIGHTS_SET[N], the quadrature weights.
//
{
  int n;
  double *w;

  n = ( ( l + 1 ) * ( l + 2 ) ) / 2;
  w = new double[n];

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

  return w;
}
//****************************************************************************80

double *padua_weights ( int l )

//****************************************************************************80
//
//  Purpose:
//
//    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 GNU LGPL license.
//
//  Modified:
//
//    11 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, int  L, the level of the set.
//    0 <= L
//
//    Output, double PADUA_WEIGHTS[(L+1)*(L+2)/2], the quadrature weights.
//
{
  double angle;
  int i;
  int i2;
  int j;
  int j2;
  int lp1h;
  int lp2h;
  int lp3h;
  double mi;
  double mj;
  double *mom;
  int n;
  const double r8_pi = 3.141592653589793;
  double *te1;
  double *te2;
  double *tmteo;
  double *tmtoe;
  double *to1;
  double *to2;
  double *w;
  double *w1;
  double *w2;

  n = ( ( l + 1 ) * ( l + 2 ) ) / 2;
  w = new double[n];

  if ( l == 0 )
  {
    w[0] = 4.0;
    return w;
  }
//
//  Relatives of L/2:
//
  lp1h = ( l + 1 ) / 2;
  lp2h = ( l + 2 ) / 2;
  lp3h = ( l + 3 ) / 2;
//
//  TE1, TE2, TO1, TO2: 
//  Even and odd Chebyshev polynomials on subgrids 1 and 2.
//
  te1 = new double[lp2h*lp2h];

  for ( j = 0; j < lp2h; j++ )
  {
    for ( i = 0; i < lp2h; i++ )
    {
      angle = r8_pi * ( double ) ( 2 * i * 2 * j ) / ( double ) ( l );
      te1[i+j*lp2h] = cos ( angle );
    }
  }

  for ( j = 0; j < lp2h; j++ )
  {
    for ( i = 1; i < lp2h; i++ )
    {
      te1[i+j*lp2h] = te1[i+j*lp2h] * sqrt ( 2.0 );
    }
  }

  to1 = new double[lp2h*lp1h];

  for ( j = 0; j < lp1h; j++ )
  {
    for ( i = 0; i < lp2h; i++ )
    {
      angle = r8_pi * ( double ) ( 2 * i * ( 2 * j + 1 ) ) 
        / ( double ) ( l );
      to1[i+j*lp2h] = cos ( angle );
    }
  }

  for ( j = 0; j < lp1h; j++ )
  {
    for ( i = 1; i < lp2h; i++ )
    {
      to1[i+j*lp2h] = to1[i+j*lp2h] * sqrt ( 2.0 );
    }
  }

  te2 = new double[lp2h*lp3h];

  for ( j = 0; j < lp3h; j++ )
  {
    for ( i = 0; i < lp2h; i++ )
    {
      angle = r8_pi * ( double ) ( 2 * i * 2 * j ) / ( double ) ( l + 1 );
      te2[i+j*lp2h] = cos ( angle );
    }
  }

  for ( j = 0; j < lp3h; j++ )
  {
    for ( i = 1; i < lp2h; i++ )
    {
      te2[i+j*lp2h] = te2[i+j*lp2h] * sqrt ( 2.0 );
    }
  }

  to2 = new double[lp2h*lp2h];

  for ( j = 0; j < lp2h; j++ )
  {
    for ( i = 0; i < lp2h; i++ )
    {
      angle = r8_pi * ( double ) ( 2 * i * ( 2 * j + 1 ) ) 
        / ( double ) ( l + 1 );
      to2[i+j*lp2h] = cos ( angle );
    }
  }

  for ( j = 0; j < lp2h; j++ )
  {
    for ( i = 1; i < lp2h; i++ )
    {
      to2[i+j*lp2h] = to2[i+j*lp2h] * sqrt ( 2.0 );
    }
  }
//
//  MOM: Moments matrix do even * even pairs.
//
  mom = new double[lp2h*lp2h];

  for ( j = 0; j < lp2h; j++ )
  {
    mj = 2.0 * sqrt ( 2.0 ) / ( double ) ( 1 - pow ( 2 * j, 2 ) );
    for ( i = 0; i < lp2h - j; i++ )
    {
      mi = 2.0 * sqrt ( 2.0 ) / ( double ) ( 1 - pow ( 2 * i, 2 ) );
      mom[i+j*lp2h] = mi * mj;
    }
  }

  i = 0;
  for ( j = 0; j < lp2h; j++ )
  {
    mom[i+j*lp2h] = mom[i+j*lp2h] / sqrt ( 2.0 );
  }

  j = 0;
  for ( i = 0; i < lp2h; i++ )
  {
    mom[i+j*lp2h] = mom[i+j*lp2h] / sqrt ( 2.0 );
  }
  if ( ( l % 2 ) == 0 )
  {
    i = lp2h-1;
    j = 0;
    mom[i+j*lp2h] = mom[i+j*lp2h] / 2.0;
  }
//
//  TMTOE and TMTEO: matrix products.
//
  tmtoe = new double[lp2h*lp2h];

  for ( j = 0; j < lp2h; j++ )
  {
    for ( i = 0; i < lp2h; i++ )
    {
      tmtoe[i+j*lp2h] = 0.0;
    }
  }

  for ( j2 = 0; j2 < lp2h; j2++ )
  {
    for ( i2 = 0; i2 < lp2h - j2; i2++ )
    {
      for ( j = 0; j < lp2h; j++ )
      {
        for ( i = 0; i < lp2h; i++ )
        {
          tmtoe[i+j*lp2h] = tmtoe[i+j*lp2h] 
            + to2[i2+i*lp2h] * mom[j2+i2*lp2h] * te1[j2+j*lp2h];
        }
      }
    }
  }

  tmteo = new double[lp3h*lp1h];

  for ( j = 0; j < lp1h; j++ )
  {
    for ( i = 0; i < lp3h; i++ )
    {
      tmteo[i+j*lp3h] = 0.0;
    }
  }

  for ( j2 = 0; j2 < lp2h; j2++ )
  {
    for ( i2 = 0; i2 < lp2h - j2; i2++ )
    {
      for ( j = 0; j < lp1h; j++ )
      {
        for ( i = 0; i < lp3h; i++ )
        {
          tmteo[i+j*lp3h] = tmteo[i+j*lp3h] 
            + te2[i2+i*lp2h] * mom[j2+i2*lp2h] * to1[j2+j*lp2h];
        }
      }
    }
  }
//
//  W1 and W2: Interpolation weight matrices.
//
  w1 = new double[lp2h*lp2h];

  for ( j = 0; j < lp2h; j++ )
  {
    for ( i = 0; i < lp2h; i++ )
    {
      w1[i+j*lp2h] = 2.0 / ( double ) ( l * ( l + 1 ) );
    }
  }

  j = 0;
  for ( i = 0; i < lp2h; i++ )
  {
    w1[i+j*lp2h] = w1[i+j*lp2h] / 2.0;
  }

  if ( ( l % 2 ) == 0 )
  {
    j = lp2h - 1;
    for ( i = 0; i < lp2h; i++ )
    {
      w1[i+j*lp2h] = w1[i+j*lp2h] / 2.0;
    }

    i = lp2h-1;
    for ( j = 0; j < lp2h; j++ )
    {
      w1[i+j*lp2h] = w1[i+j*lp2h] / 2.0;
    }
  }

  w2 = new double[lp3h*lp1h];

  for ( j = 0; j < lp1h; j++ )
  {
    for ( i = 0; i < lp3h; i++ )
    {
      w2[i+j*lp3h] = 2.0 / ( double ) ( l * ( l + 1 ) );
    }
  }

  i = 0;
  for ( j = 0; j < lp1h; j++ )
  {
    w2[i+j*lp3h] = w2[i+j*lp3h] / 2.0;
  }

  if ( ( l % 2 ) == 1 )
  {
    i = lp3h - 1;
    for ( j = 0; j < lp1h; j++ )
    {
      w2[i+j*lp3h] = w2[i+j*lp3h] / 2.0;
    }
    j = lp1h - 1;
    for ( i = 0; i < lp3h; i++ )
    {
      w2[i+j*lp3h] = w2[i+j*lp3h] / 2.0;
    }
  }
//
//  Cubature weights as matrices on the subgrids.
//
  for ( j = 0; j < lp2h; j++ )
  {
    for ( i = 0; i < lp2h; i++ )
    {
      w1[i+j*lp2h] = w1[i+j*lp2h] * tmtoe[i+j*lp2h];
    }
  }

  for ( j = 0; j < lp1h; j++ )
  {
    for ( i = 0; i < lp3h; i++ )
    {
      w2[i+j*lp3h] = w2[i+j*lp3h] * tmteo[i+j*lp3h];
    }
  }
//
//  Pack weight matrices W1 and W2 into the vector W.
//
  if ( ( l % 2 ) == 0 )
  {
    for ( j = 0; j < lp2h; j++ )
    {
      for ( i = 0; i < lp2h; i++ )
      {
        w[i+2*j*lp2h] = w1[i+j*lp2h];
      }
    }

    for ( j = 0; j < lp1h; j++ )
    {
      for ( i = 0; i < lp3h; i++ )
      {
        w[i+(2*j+1)*lp2h] = w2[i+j*lp3h];
      }
    }
  }
  else
  {
    for ( j = 0; j < lp1h; j++ )
    {
      for ( i = 0; i < lp2h; i++ )
      {
        w[i+j*(l+2)] = w1[i+j*lp2h];
      }
    }

    for ( j = 0; j < lp1h; j++ )
    {
      for ( i = 0; i < lp3h; i++ )
      {
        w[i+lp2h+j*(l+2)] = w2[i+j*lp3h];
      }
    }

  }
//
//  Free memory.
//
  delete [] te1;
  delete [] te2;
  delete [] tmteo;
  delete [] tmtoe;
  delete [] to1;
  delete [] to2;
  delete [] w1;
  delete [] w2;

  return w;
}
//****************************************************************************80

double r8_max ( double x, double y )

//****************************************************************************80
//
//  Purpose:
//
//    R8_MAX returns the maximum of two R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    18 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, Y, the quantities to compare.
//
//    Output, double R8_MAX, the maximum of X and Y.
//
{
  double value;

  if ( y < x )
  {
    value = x;
  }
  else
  {
    value = y;
  }
  return value;
}
//****************************************************************************80

void r8mat_transpose_print ( int m, int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    10 September 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Input, double A[M*N], an M by N matrix to be printed.
//
//    Input, string TITLE, a title.
//
{
  r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title );

  return;
}
//****************************************************************************80

void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo,
  int ihi, int jhi, string title )

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 April 2014
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, N, the number of rows and columns.
//
//    Input, double A[M*N], an M by N matrix to be printed.
//
//    Input, int ILO, JLO, the first row and column to print.
//
//    Input, int IHI, JHI, the last row and column to print.
//
//    Input, string TITLE, a title.
//
{
# define INCX 5

  int i;
  int i2;
  int i2hi;
  int i2lo;
  int i2lo_hi;
  int i2lo_lo;
  int inc;
  int j;
  int j2hi;
  int j2lo;

  cout << "\n";
  cout << title << "\n";

  if ( m <= 0 || n <= 0 )
  {
    cout << "\n";
    cout << "  (None)\n";
    return;
  }

  if ( ilo < 1 )
  {
    i2lo_lo = 1;
  }
  else
  {
    i2lo_lo = ilo;
  }

  if ( ihi < m )
  {
    i2lo_hi = m;
  }
  else
  {
    i2lo_hi = ihi;
  }

  for ( i2lo = i2lo_lo; i2lo <= i2lo_hi; i2lo = i2lo + INCX )
  {
    i2hi = i2lo + INCX - 1;

    if ( m < i2hi )
    {
      i2hi = m;
    }
    if ( ihi < i2hi )
    {
      i2hi = ihi;
    }

    inc = i2hi + 1 - i2lo;

    cout << "\n";
    cout << "  Row: ";
    for ( i = i2lo; i <= i2hi; i++ )
    {
      cout << setw(7) << i - 1 << "       ";
    }
    cout << "\n";
    cout << "  Col\n";
    cout << "\n";

    if ( jlo < 1 )
    {
      j2lo = 1;
    }
    else
    {
      j2lo = jlo;
    }
    if ( n < jhi )
    {
      j2hi = n;
    }
    else
    {
      j2hi = jhi;
    }

    for ( j = j2lo; j <= j2hi; j++ )
    {
      cout << setw(5) << j - 1 << ":";
      for ( i2 = 1; i2 <= inc; i2++ )
      {
        i = i2lo - 1 + i2;
        cout << setw(14) << a[(i-1)+(j-1)*m];
      }
      cout << "\n";
    }
  }

  return;
# undef INCX
}
//****************************************************************************80

void r8vec_reverse ( int n, double a[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_REVERSE reverses the elements of an R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  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 GNU LGPL license.
//
//  Modified:
//
//    18 September 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the array.
//
//    Input/output, double A[N], the array to be reversed.
//
{
  int i;
  int i_hi;
  double temp;

  i_hi = n / 2;

  for ( i = 1; i <= i_hi; i++ )
  {
    temp   = a[i-1];
    a[i-1] = a[n-i];
    a[n-i] = temp;
  }

  return;
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    08 July 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  size_t len;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  len = std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}