# include <cmath>
# include <cstdlib>
# include <ctime>
# include <iomanip>
# include <iostream>

using namespace std;

# include "triangle_symq_rule.hpp"

int main ( );
void triangle_symq_rule_test01 ( int p );
void triangle_symq_rule_test02 ( int p );
void triangle_symq_rule_test03 ( int p_lo, int p_hi );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    triangle_symq_rule_test() tests triangle_symq_rule().
//
//  Licensing:
//
//    This code is distributed under the GNU GPL license.
//
//  Modified:
//
//    10 July 2023
//
//  Author:
//
//    John Burkardt.
//
{
  int p;
  int p_hi;
  int p_lo;

  timestamp ( );
  cout << "\n";
  cout << "triangle_symq_rule_test():\n";
  cout << "  C++ version\n";
  cout << "  Test triangle_symq_rule().\n";

  p = 5;
  triangle_symq_rule_test01 ( p );

  p = 5;
  triangle_symq_rule_test02 ( p );

  p_lo = 0;
  p_hi = 50;
  triangle_symq_rule_test03 ( p_lo, p_hi );
//
//  Terminate.
//
  cout << "\n";
  cout << "triangle_symq_rule_test()\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void triangle_symq_rule_test01 ( int p )

//****************************************************************************80
//
//  Purpose:
//
//    triangle_symq_rule_test01() prints the rule of precision P.
//
//  Licensing:
//
//    This code is distributed under the GNU GPL license.
//
//  Modified:
//
//    11 June 2023
//
//  Author:
//
//    John Burkardt.
//
{
  double *a;
  double *b;
  double *c;
  int i;
  int n;
  double *w;
  double w_sum;

  cout << "\n";
  cout << "triangle_symq_rule_test01():\n";
  cout << "  Quadrature rule for the triangle,\n";
  cout << "  given in barycentric coordinates.\n";
  cout << "  Precision p = " << p << "\n";
//
//  Retrieve the quadrature rule.
//
  n = rule_order ( p );
  cout << "\n";
  cout << "  Number of nodes N = " << n << "\n";

  a = new double[n];
  b = new double[n];
  c = new double[n];
  w = new double[n];

  triangle_symq_rule ( p, n, a, b, c, w );
//
//  Print the rule.
//
  cout << "\n";
  cout << "     I      W          A           B           C\n";
  cout << "\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  " << setw(4) << i
         << "  " << setw(10) << w[i]
         << "  " << setw(10) << a[i]
         << "  " << setw(10) << b[i]
         << "  " << setw(10) << c[i] << "\n";
  }
//
//  Verify weights sum to 1.
//
  w_sum = 0.0;
  for ( i = 0; i < n; i++ )
  {
    w_sum = w_sum + w[i];
  }

  cout << "\n";
  cout << "  Weight Sum = " << w_sum << "\n";
//
//  Free memory.
//
  delete [] a;
  delete [] b;
  delete [] c;
  delete [] w;

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

void triangle_symq_rule_test02 ( int p )

//****************************************************************************80
//
//  Purpose:
//
//    triangle_symq_rule_test02() tests a rule of precision P.
//
//  Licensing:
//
//    This code is distributed under the GNU GPL license.
//
//  Modified:
//
//    11 June 2023
//
//  Author:
//
//    John Burkardt.
//
{
  double *a;
  double *b;
  double *c;
  int degree;
  const int dim_num = 2;
  double exact;
  int expon[2];
  int h;
  int i;
  int ij;
  double max_error;
  bool more;
  double quad;
  double quad_error;
  int n;
  int t;
  double *v;
  double *w;
  double *xy;

  cout << "\n";
  cout << "triangle_symq_rule_test02():\n";
  cout << "  Test the precision of a quadrature rule for the unit triangle.\n";
//
//  Retrieve the quadrature rule.
//
  n = rule_order ( p );
  cout << "\n";
  cout << "  Number of nodes N = " << n << "\n";

  a = new double[n];
  b = new double[n];
  c = new double[n];
  w = new double[n];

  triangle_symq_rule ( p, n, a, b, c, w );
//
//  Pack the x, y, z vectors as columns of an array.
//
  xy = new double[n*2];
  ij = 0;
  for ( i = 0; i < n; i++ )
  {
    xy[ij] = a[i];
    ij = ij + 1;
    xy[ij] = b[i];
    ij = ij + 1;
  }

  cout << "\n";
  cout << "  Stated precision of rule    = " << p << "\n";
  cout << "  Number of quadrature points = " << n << "\n";
  cout << "\n";
  cout << "  Degree  Maximum error\n";
  cout << "\n";

  v = new double[n];

  for ( degree = 0; degree <= p + 2; degree++ )
  {
    more = 0;
    h = 0;
    t = 0;
    max_error = 0.0;

    while ( true )
    {
      comp_next ( degree, dim_num, expon, more, h, t );

      v = monomial_value ( dim_num, n, expon, xy );

      quad = 0.0;
      for ( i = 0; i < n; i++ )
      {
        quad = quad + w[i] * v[i];
      }
      quad = triangle_unit_area ( ) * quad;

      exact = triangle_unit_monomial_integral ( expon );

      quad_error = fabs ( quad - exact );
      max_error = fmax ( max_error, quad_error );

      if ( ! more )
      {
        break;
      }
    }
    cout << "  " << setw(2) << degree
         << "  " << setw(24) << max_error << "\n";
  }
//
//  Free memory.
//
  delete [] a;
  delete [] b;
  delete [] c;
  delete [] v;
  delete [] w;
  delete [] xy;

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

void triangle_symq_rule_test03 ( int p_lo, int p_hi )

//****************************************************************************80
//
//  Purpose:
//
//    triangle_symq_rule_test03() tests absolute and relative precision.
//
//  Licensing:
//
//    This code is distributed under the GNU GPL license.
//
//  Modified:
//
//    10 July 2023
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    integer p_lo, p_hi: the lowest and highest rules to check.
//
{
  double *a;
  double *b;
  double *c;
  int degree;
  const int dim_num = 2;
  double exact;
  int expon[2];
  int h;
  int i;
  int ij;
  double max_abs;
  double max_rel;
  bool more;
  int p;
  double quad;
  double quad_error;
  int n;
  int t;
  double *v;
  double *w;
  double *xy;

  cout << "\n";
  cout << "triangle_symq_rule_test03():\n";
  cout << "  Test the precision of a quadrature rule for the unit triangle.\n";
  cout << "  Check rules of precision p = " << p_lo << " through " << p_hi << "\n";
  cout << "  for error in approximating integrals of monomials.\n";

  cout << "\n";
  cout << "              maximum                   maximum\n";
  cout << "   p          absolute                  relative\n";
  cout << "              error                     error\n";
  cout << "\n";

  for ( p = p_lo; p <= p_hi; p++ )
  {
    n = rule_order ( p );

    a = new double[n];
    b = new double[n];
    c = new double[n];
    w = new double[n];

    triangle_symq_rule ( p, n, a, b, c, w );

    xy = new double[n*2];
    ij = 0;
    for ( i = 0; i < n; i++ )
    {
      xy[ij] = a[i];
      ij = ij + 1;
      xy[ij] = b[i];
      ij = ij + 1;
    }

    v = new double[n];

    max_abs = 0.0;
    max_rel = 0.0;

    for ( degree = 0; degree <= p; degree++ )
    {
      more = 0;
      h = 0;
      t = 0;

      while ( true )
      {
        comp_next ( degree, dim_num, expon, more, h, t );

        v = monomial_value ( dim_num, n, expon, xy );

        quad = 0.0;
        for ( i = 0; i < n; i++ )
        {
          quad = quad + w[i] * v[i];
        }
        quad = triangle_unit_area ( ) * quad;

        exact = triangle_unit_monomial_integral ( expon );

        quad_error = fabs ( quad - exact );
        max_abs = fmax ( max_abs, quad_error );
        max_rel = fmax ( max_rel, quad_error / exact );

        if ( ! more )
        {
          break;
        }

      }
    }
    cout << "  " << setw(2) << p
         << "  " << setw(24) << max_abs
         << "  " << setw(24) << max_rel << "\n";
//
//  Free memory.
//
    delete [] a;
    delete [] b;
    delete [] c;
    delete [] v;
    delete [] w;
    delete [] xy;
  }

  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 MIT license. 
//
//  Modified:
//
//    19 March 2018
//
//  Author:
//
//    John Burkardt
//
{
# 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
}