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

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

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_points_set ( int l )

//****************************************************************************80
//
//  Purpose:
//
//    PADUA_POINTS_SET sets the Padua points.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    30 August 2016
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Marco Caliari, Stefano de Marchi, Marco Vianello,
//    Bivariate interpolation on the square at new nodal sets,
//    Applied Mathematics and Computation,
//    Volume 165, Number 2, 2005, pages 261-274.
//
//  Parameters:
//
//    Input, int L, the level.
//    0 <= L <= 10.
//
//    Output, double PADUA_POINTS_SET[2*N], the Padua points.
//
{
  int j1;
  int j1_hi;
  int j2;
  int n;
  double t;
  double *xy;
  double xy00[2*1] = {
      0.000000000000000, 0.000000000000000 };
  double xy01[2*3] = {
      -1.000000000000000, -1.000000000000000,  
      -1.000000000000000,  1.000000000000000,
       1.000000000000000,  0.000000000000000 };
  double xy02[2*6] = {
      -1.000000000000000, -1.000000000000000,
      -1.000000000000000,  0.5000000000000001,
       0.000000000000000, -0.4999999999999998,  
       0.000000000000000,  1.000000000000000,    
       1.000000000000000, -1.000000000000000,
       1.000000000000000,  0.5000000000000001 };
  double xy03[2*10] = {
      -1.000000000000000,  -1.000000000000000,    
      -1.000000000000000,   0.000000000000000,
      -1.000000000000000,   1.000000000000000,    
      -0.4999999999999998, -0.7071067811865475,    
      -0.4999999999999998,  0.7071067811865476,    
       0.5000000000000001, -1.000000000000000,    
       0.5000000000000001,  0.000000000000000,    
       0.5000000000000001,  1.000000000000000,    
       1.000000000000000,  -0.7071067811865475,    
       1.000000000000000,   0.7071067811865476 };
  double xy04[2*15] = {
      -1.000000000000000,  -1.000000000000000,
      -1.000000000000000,  -0.3090169943749473,
      -1.000000000000000,   0.8090169943749475,
      -0.7071067811865475, -0.8090169943749473,
      -0.7071067811865475,  0.3090169943749475,
      -0.7071067811865475,  1.000000000000000,
       0.000000000000000,  -1.000000000000000,
       0.000000000000000,  -0.3090169943749473,
       0.000000000000000,   0.8090169943749475,
       0.7071067811865476, -0.8090169943749473,
       0.7071067811865476,  0.3090169943749475,
       0.7071067811865476,  1.000000000000000,
       1.000000000000000,  -1.000000000000000,
       1.000000000000000,  -0.3090169943749473,
       1.000000000000000,   0.8090169943749475 };
  double xy05[2*21] = {
      -1.000000000000000,  -1.000000000000000,
      -1.000000000000000,  -0.4999999999999998,   
      -1.000000000000000,   0.5000000000000001,    
      -1.000000000000000,   1.000000000000000,   
      -0.8090169943749473, -0.8660254037844387,
      -0.8090169943749473,  0.000000000000000, 
      -0.8090169943749473,  0.8660254037844387,
      -0.3090169943749473, -1.000000000000000, 
      -0.3090169943749473, -0.4999999999999998,   
      -0.3090169943749473,  0.5000000000000001,
      -0.3090169943749473,  1.000000000000000, 
       0.3090169943749475, -0.8660254037844387,
       0.3090169943749475,  0.000000000000000,  
       0.3090169943749475,  0.8660254037844387,   
       0.8090169943749475, -1.000000000000000,  
       0.8090169943749475, -0.4999999999999998,  
       0.8090169943749475,  0.5000000000000001,
       0.8090169943749475,  1.000000000000000,
       1.000000000000000,  -0.8660254037844387,   
       1.000000000000000,   0.000000000000000,
       1.000000000000000,   0.8660254037844387 };
  double xy06[2*28] = {
      -1.000000000000000,  -1.000000000000000,  
      -1.000000000000000,  -0.6234898018587335,    
      -1.000000000000000,   0.2225209339563144,    
      -1.000000000000000,   0.9009688679024191,    
      -0.8660254037844387, -0.9009688679024190,    
      -0.8660254037844387, -0.2225209339563143, 
      -0.8660254037844387,  0.6234898018587336,
      -0.8660254037844387,  1.000000000000000,
      -0.4999999999999998, -1.000000000000000,  
      -0.4999999999999998, -0.6234898018587335,    
      -0.4999999999999998,  0.2225209339563144,    
      -0.4999999999999998,  0.9009688679024191,    
       0.000000000000000,  -0.9009688679024190, 
       0.000000000000000,  -0.2225209339563143,
       0.000000000000000,   0.6234898018587336,    
       0.000000000000000,   1.000000000000000,
       0.5000000000000001, -1.000000000000000,   
       0.5000000000000001, -0.6234898018587335,    
       0.5000000000000001,  0.2225209339563144,    
       0.5000000000000001,  0.9009688679024191,    
       0.8660254037844387, -0.9009688679024190,    
       0.8660254037844387, -0.2225209339563143, 
       0.8660254037844387,  0.6234898018587336,
       0.8660254037844387,  1.000000000000000,    
       1.000000000000000,  -1.000000000000000,  
       1.000000000000000,  -0.6234898018587335, 
       1.000000000000000,   0.2225209339563144,
       1.000000000000000,   0.9009688679024191 };
  double xy07[2*36] = {
      -1.000000000000000, -1.000000000000000,
      -1.000000000000000, -0.7071067811865475,    
      -1.000000000000000,  0.000000000000000,
      -1.000000000000000,  0.7071067811865476, 
      -1.000000000000000,  1.000000000000000,  
      -0.9009688679024190, -0.9238795325112867,    
      -0.9009688679024190, -0.3826834323650897,    
      -0.9009688679024190,  0.3826834323650898, 
      -0.9009688679024190,  0.9238795325112867,
      -0.6234898018587335, -1.000000000000000,
      -0.6234898018587335, -0.7071067811865475,
      -0.6234898018587335,  0.000000000000000,
      -0.6234898018587335,  0.7071067811865476,
      -0.6234898018587335, 1.000000000000000, 
      -0.2225209339563143, -0.9238795325112867,    
      -0.2225209339563143, -0.3826834323650897,    
      -0.2225209339563143,  0.3826834323650898, 
      -0.2225209339563143,  0.9238795325112867,
       0.2225209339563144, -1.000000000000000,
       0.2225209339563144, -0.7071067811865475,
       0.2225209339563144,  0.000000000000000,
       0.2225209339563144,  0.7071067811865476,
       0.2225209339563144,  1.000000000000000,  
       0.6234898018587336, -0.9238795325112867,    
       0.6234898018587336, -0.3826834323650897,    
       0.6234898018587336,  0.3826834323650898, 
       0.6234898018587336,  0.9238795325112867,
       0.9009688679024191, -1.000000000000000,
       0.9009688679024191, -0.7071067811865475,
       0.9009688679024191,  0.000000000000000,
       0.9009688679024191,  0.7071067811865476,
       0.9009688679024191,  1.000000000000000,   
       1.000000000000000, -0.9238795325112867,    
       1.000000000000000, -0.3826834323650897,    
       1.000000000000000,  0.3826834323650898,
       1.000000000000000,  0.9238795325112867 };
  double xy08[2*45] = {
      -1.000000000000000, -1.000000000000000,   
      -1.000000000000000, -0.7660444431189779,    
      -1.000000000000000, -0.1736481776669303,    
      -1.000000000000000,  0.5000000000000001,    
      -1.000000000000000,  0.9396926207859084,    
      -0.9238795325112867, -0.9396926207859083,    
      -0.9238795325112867, -0.4999999999999998,    
      -0.9238795325112867,  0.1736481776669304,    
      -0.9238795325112867,  0.7660444431189780,
      -0.9238795325112867,  1.000000000000000,
      -0.7071067811865475, -1.000000000000000,   
      -0.7071067811865475, -0.7660444431189779,    
      -0.7071067811865475, -0.1736481776669303,    
      -0.7071067811865475,  0.5000000000000001,    
      -0.7071067811865475,  0.9396926207859084,    
      -0.3826834323650897, -0.9396926207859083,    
      -0.3826834323650897, -0.4999999999999998,    
      -0.3826834323650897,  0.1736481776669304,    
      -0.3826834323650897,  0.7660444431189780,
      -0.3826834323650897,  1.000000000000000,    
       0.000000000000000, -1.000000000000000,   
       0.000000000000000, -0.7660444431189779,    
       0.000000000000000, -0.1736481776669303,    
       0.000000000000000,  0.5000000000000001,    
       0.000000000000000,  0.9396926207859084,    
       0.3826834323650898, -0.9396926207859083,    
       0.3826834323650898, -0.4999999999999998,    
       0.3826834323650898,  0.1736481776669304,    
       0.3826834323650898,  0.7660444431189780,
       0.3826834323650898,  1.000000000000000,
       0.7071067811865476, -1.000000000000000,   
       0.7071067811865476, -0.7660444431189779,    
       0.7071067811865476, -0.1736481776669303,    
       0.7071067811865476,  0.5000000000000001,    
       0.7071067811865476,  0.9396926207859084,    
       0.9238795325112867, -0.9396926207859083,    
       0.9238795325112867, -0.4999999999999998,    
       0.9238795325112867,  0.1736481776669304,    
       0.9238795325112867,  0.7660444431189780,
       0.9238795325112867,  1.000000000000000,    
       1.000000000000000, -1.000000000000000,   
       1.000000000000000, -0.7660444431189779,    
       1.000000000000000, -0.1736481776669303,  
       1.000000000000000,  0.5000000000000001,
       1.000000000000000,  0.9396926207859084 };
  double xy09[2*55] = {
      -1.000000000000000, -1.000000000000000,   
      -1.000000000000000, -0.8090169943749473,    
      -1.000000000000000, -0.3090169943749473,    
      -1.000000000000000,  0.3090169943749475,
      -1.000000000000000,  0.8090169943749475,    
      -1.000000000000000,  1.000000000000000,   
      -0.9396926207859083, -0.9510565162951535,    
      -0.9396926207859083, -0.5877852522924730,
      -0.9396926207859083,  0.000000000000000,   
      -0.9396926207859083,  0.5877852522924731,    
      -0.9396926207859083,  0.9510565162951535,
      -0.7660444431189779, -1.000000000000000,   
      -0.7660444431189779, -0.8090169943749473,    
      -0.7660444431189779, -0.3090169943749473,    
      -0.7660444431189779,  0.3090169943749475,    
      -0.7660444431189779,  0.8090169943749475, 
      -0.7660444431189779,  1.000000000000000,   
      -0.4999999999999998, -0.9510565162951535,    
      -0.4999999999999998, -0.5877852522924730,
      -0.4999999999999998,  0.000000000000000,   
      -0.4999999999999998,  0.5877852522924731,    
      -0.4999999999999998,  0.9510565162951535,
      -0.1736481776669303, -1.000000000000000,   
      -0.1736481776669303, -0.8090169943749473,    
      -0.1736481776669303, -0.3090169943749473,    
      -0.1736481776669303,  0.3090169943749475,    
      -0.1736481776669303,  0.8090169943749475,
      -0.1736481776669303,  1.000000000000000,   
       0.1736481776669304, -0.9510565162951535,    
       0.1736481776669304, -0.5877852522924730,
       0.1736481776669304,  0.000000000000000,   
       0.1736481776669304,  0.5877852522924731,    
       0.1736481776669304,  0.9510565162951535,
       0.5000000000000001, -1.000000000000000,   
       0.5000000000000001, -0.8090169943749473,    
       0.5000000000000001, -0.3090169943749473,    
       0.5000000000000001,  0.3090169943749475,    
       0.5000000000000001,  0.8090169943749475,
       0.5000000000000001,  1.000000000000000,   
       0.7660444431189780, -0.9510565162951535,    
       0.7660444431189780, -0.5877852522924730,
       0.7660444431189780,  0.000000000000000,   
       0.7660444431189780,  0.5877852522924731,    
       0.7660444431189780,  0.9510565162951535,
       0.9396926207859084, -1.000000000000000,   
       0.9396926207859084, -0.8090169943749473,    
       0.9396926207859084, -0.3090169943749473,    
       0.9396926207859084,  0.3090169943749475,    
       0.9396926207859084,  0.8090169943749475,
       0.9396926207859084,  1.000000000000000,   
       1.000000000000000, -0.9510565162951535,
       1.000000000000000, -0.5877852522924730,    
       1.000000000000000,  0.000000000000000,
       1.000000000000000,  0.5877852522924731,
       1.000000000000000,  0.9510565162951535 };
  double xy10[2*66] = {
      -1.000000000000000, -1.000000000000000,   
      -1.000000000000000, -0.8412535328311811,    
      -1.000000000000000, -0.4154150130018863,    
      -1.000000000000000,  0.1423148382732851,    
      -1.000000000000000,  0.6548607339452851,    
      -1.000000000000000,  0.9594929736144974,    
      -0.9510565162951535, -0.9594929736144974,    
      -0.9510565162951535, -0.6548607339452850,    
      -0.9510565162951535, -0.1423148382732850,    
      -0.9510565162951535,  0.4154150130018864,    
      -0.9510565162951535,  0.8412535328311812,
      -0.9510565162951535,  1.000000000000000,
      -0.8090169943749473, -1.000000000000000,   
      -0.8090169943749473, -0.8412535328311811,    
      -0.8090169943749473, -0.4154150130018863,    
      -0.8090169943749473,  0.1423148382732851,    
      -0.8090169943749473,  0.6548607339452851,    
      -0.8090169943749473,  0.9594929736144974,    
      -0.5877852522924730, -0.9594929736144974,    
      -0.5877852522924730, -0.6548607339452850,    
      -0.5877852522924730, -0.1423148382732850,    
      -0.5877852522924730,  0.4154150130018864,    
      -0.5877852522924730,  0.8412535328311812,
      -0.5877852522924730,  1.000000000000000,
      -0.3090169943749473, -1.000000000000000,   
      -0.3090169943749473, -0.8412535328311811,    
      -0.3090169943749473, -0.4154150130018863,    
      -0.3090169943749473,  0.1423148382732851,    
      -0.3090169943749473,  0.6548607339452851,   
      -0.3090169943749473, 0.9594929736144974,    
       0.000000000000000, -0.9594929736144974,    
       0.000000000000000, -0.6548607339452850,    
       0.000000000000000, -0.1423148382732850,    
       0.000000000000000,  0.4154150130018864,
       0.000000000000000,  0.8412535328311812,    
       0.000000000000000,  1.000000000000000, 
       0.3090169943749475, -1.000000000000000,    
       0.3090169943749475, -0.8412535328311811,    
       0.3090169943749475, -0.4154150130018863,    
       0.3090169943749475,  0.1423148382732851,    
       0.3090169943749475,  0.6548607339452851,    
       0.3090169943749475,  0.9594929736144974,    
       0.5877852522924731, -0.9594929736144974,    
       0.5877852522924731, -0.6548607339452850,    
       0.5877852522924731, -0.1423148382732850,    
       0.5877852522924731,  0.4154150130018864,    
       0.5877852522924731,  0.8412535328311812,
       0.5877852522924731, 1.000000000000000,
       0.8090169943749475, -1.000000000000000,   
       0.8090169943749475, -0.8412535328311811,    
       0.8090169943749475, -0.4154150130018863,    
       0.8090169943749475,  0.1423148382732851,    
       0.8090169943749475,  0.6548607339452851,    
       0.8090169943749475,  0.9594929736144974,    
       0.9510565162951535, -0.9594929736144974,    
       0.9510565162951535, -0.6548607339452850,    
       0.9510565162951535, -0.1423148382732850,    
       0.9510565162951535,  0.4154150130018864,    
       0.9510565162951535,  0.8412535328311812,
       0.9510565162951535,  1.000000000000000,    
       1.000000000000000, -1.000000000000000,   
       1.000000000000000, -0.8412535328311811,    
       1.000000000000000, -0.4154150130018863,    
       1.000000000000000,  0.1423148382732851,    
       1.000000000000000,  0.6548607339452851,
       1.000000000000000,  0.9594929736144974 };

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

  if ( l == 0 )
  {
    xy = r8vec_copy_new ( 2 * n, xy00 );
  }
  else if ( l == 1 )
  {
    xy = r8vec_copy_new ( 2 * n, xy01 ); 
  }
  else if ( l == 2 )
  {
    xy = r8vec_copy_new ( 2 * n, xy02 );  
  }
  else if ( l == 3 )
  {
    xy = r8vec_copy_new ( 2 * n, xy03 );   
  }
  else if ( l == 4 )
  {
    xy = r8vec_copy_new ( 2 * n, xy04 );   
  }
  else if ( l == 5 )
  {
    xy = r8vec_copy_new ( 2 * n, xy05 );
  }
  else if ( l == 6 )
  {   
    xy = r8vec_copy_new ( 2 * n, xy06 );
  }
  else if ( l == 7 )
  {    
    xy = r8vec_copy_new ( 2 * n, xy07 );
  }
  else if ( l == 8 )
  { 
    xy = r8vec_copy_new ( 2 * n, xy08 );
  }
  else if ( l == 9 )
  { 
    xy = r8vec_copy_new ( 2 * n, xy09 ); 
  }
  else if ( l == 10 )
  {
    xy = r8vec_copy_new ( 2 * n, xy10 );
  }
  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 data to match published information.
//
  j1_hi = ( n - 1 ) / 2;

  for ( j1 = 0; j1 < j1_hi; j1++ )
  {
    j2 = n - 1 - j1;
    t          = xy[0+2*j1];
    xy[0+2*j1] = xy[0+2*j2];
    xy[0+2*j2] = t;
    t          = xy[1+2*j1];
    xy[1+2*j1] = xy[1+2*j2];
    xy[1+2*j2] = t;
  }

  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

double *r8vec_copy_new ( int n, double a1[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_COPY_NEW copies an R8VEC to a new R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 July 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vectors.
//
//    Input, double A1[N], the vector to be copied.
//
//    Output, double R8VEC_COPY_NEW[N], the copy of A1.
//
{
  double *a2;
  int i;

  a2 = new double[n];

  for ( i = 0; i < n; i++ )
  {
    a2[i] = a1[i];
  }
  return a2;
}
//****************************************************************************80

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

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_PRINT prints an R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    16 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of components of the vector.
//
//    Input, double A[N], the vector to be printed.
//
//    Input, string TITLE, a title.
//
{
  int i;

  cout << "\n";
  cout << title << "\n";
  cout << "\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  " << setw(8)  << i
         << ": " << setw(14) << a[i]  << "\n";
  }

  return;
}
//****************************************************************************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;
  std::time_t now;

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

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

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

  return;
# undef TIME_SIZE
}