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

using namespace std;

int main ( int argc, char *argv[] );
void clenshaw_curtis_compute ( int order, double xtab[], double weight[] );
void clenshaw_curtis_handle ( int order, string output );
void r8mat_write ( string output_filename, int m, int n, double table[] );
void timestamp ( );

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

int main ( int argc, char *argv[] )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for CLENSHAW_CURTIS_RULE.
//
//  Discussion:
//
//    This program computes a standard Clenshaw Curtis quadrature rule
//    and writes it to a file.
//
//    The user specifies:
//    * the ORDER (number of points) in the rule
//    * the OUTPUT option.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    27 June 2009
//
//  Author:
//
//    John Burkardt
//
{
  int order;
  string prefix;

  timestamp ( );
  cout << "\n";
  cout << "CLENSHAW_CURTIS_RULE\n";
  cout << "  C++ version\n";
  cout << "\n";
  cout << "  Compiled on " << __DATE__ << " at " << __TIME__ << ".\n";
  cout << "\n";
  cout << "  Compute a Clenshaw Curtis rule for approximating\n";
  cout << "\n";
  cout << "    Integral ( -1 <= x <= +1 ) f(x) dx\n";
  cout << "\n";
  cout << "  of order ORDER.\n";
  cout << "\n";
  cout << "  The user specifies ORDER and OUTPUT.\n";
  cout << "\n";
  cout << "  OUTPUT is:\n";
  cout << "\n";
  cout << "  \"C++\" for printed C++ output;\n";
  cout << "  \"F77\" for printed Fortran77 output;\n";
  cout << "  \"F90\" for printed Fortran90 output;\n";
  cout << "  \"MAT\" for printed MATLAB output;\n";
  cout << "\n";
  cout << "  or:\n";
  cout << "\n";
  cout << "  \"filename\" to generate 3 files:\n";
  cout << "\n";
  cout << "    filename_w.txt - the weight file\n";
  cout << "    filename_x.txt - the abscissa file.\n";
  cout << "    filename_r.txt - the region file.\n";
//
//  Get the order.
//
  if ( 1 < argc )
  {
    order = atoi ( argv[1] );
  }
  else
  {
    cout << "\n";
    cout << "  Enter the value of ORDER (1 or greater)\n";
    cin >> order;
  }

  cout << "\n";
  cout << "  The requested order of the rule is = " << order << "\n";
//
//  Get the output option or quadrature file root name:
//
  if ( 2 < argc )
  {
    prefix = argv[2];
  }
  else
  {
    cout << "\n";
    cout << "  Enter OUTPUT (one of C++, F77, F90, MAT\n";
    cout << "  or else the \"root name\" of the quadrature files).\n";
    cin >> prefix;
  }

  cout << "\n";
  cout << "  OUTPUT option is \"" << prefix << "\".\n";
//
//  Construct the rule and output it.
//
  clenshaw_curtis_handle ( order, prefix );

  cout << "\n";
  cout << "CLENSHAW_CURTIS_RULE:\n";
  cout << "  Normal end of execution.\n";

  cout << "\n";
  timestamp ( );

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

void clenshaw_curtis_compute ( int order, double x[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    CLENSHAW_CURTIS_COMPUTE computes a Clenshaw Curtis quadrature rule.
//
//  Discussion:
//
//    The integration interval is [ -1, 1 ].
//
//    The weight function is w(x) = 1.0.
//
//    The integral to approximate:
//
//      Integral ( -1 <= X <= 1 ) F(X) dX
//
//    The quadrature rule:
//
//      Sum ( 1 <= I <= ORDER ) W(I) * F ( X(I) )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    19 March 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int ORDER, the order of the rule.
//    1 <= ORDER.
//
//    Output, double X[ORDER], the abscissas.
//
//    Output, double W[ORDER], the weights.
//
{
  double b;
  int i;
  int j;
  double pi = 3.141592653589793;
  double theta;

  if ( order < 1 )
  {
    std::cerr << "\n";
    std::cerr << "CLENSHAW_CURTIS_COMPUTE - Fatal error!\n";
    std::cerr << "  Illegal value of ORDER = " << order << "\n";
    std::exit ( 1 );
  }
  else if ( order == 1 )
  {
    x[0] = 0.0;
    w[0] = 2.0;
  }
  else
  {
    for ( i = 0; i < order; i++ )
    {
      x[i] =  std::cos ( ( double ) ( order - 1 - i ) * pi 
                       / ( double ) ( order - 1     ) );
    }
    x[0] = -1.0;
    if ( ( order % 2 ) == 1 )
    {
      x[(order-1)/2] = 0.0;
    }
    x[order-1] = +1.0;

    for ( i = 0; i < order; i++ )
    {
      theta = ( double ) ( i ) * pi / ( double ) ( order - 1 );

      w[i] = 1.0;

      for ( j = 1; j <= ( order - 1 ) / 2; j++ )
      {
        if ( 2 * j == ( order - 1 ) )
        {
          b = 1.0;
        }
        else
        {
          b = 2.0;
        }

        w[i] = w[i] - b *  std::cos ( 2.0 * ( double ) ( j ) * theta ) 
          / ( double ) ( 4 * j * j - 1 );
      }
    }

    w[0] = w[0] / ( double ) ( order - 1 );
    for ( i = 1; i < order - 1; i++ )
    {
      w[i] = 2.0 * w[i] / ( double ) ( order - 1 );
    }
    w[order-1] = w[order-1] / ( double ) ( order - 1 );
  }

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

void clenshaw_curtis_handle ( int order, string prefix )

//****************************************************************************80
//
//  Purpose:
//
//    CLENSHAW_CURTIS_HANDLE computes a Clenshaw Curtis rule and outputs it.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    27 June 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int ORDER, the order of the rule.
//
//    Input, string PREFIX, specifies the output.
//    * "C++'", print as C++ code.
//    * "F77", print as FORTRAN77 code.
//    * "F90", print as FORTRAN90 code.
//    * "MAT", print as MATLAB code.
//    * file,  write files "file_w.txt", "file_x.txt", "file_r.txt" 
//      defining weights, abscissas, and region.
// 
{
  int i;
  string output_r;
  string output_w;
  string output_x;
  double *r;
  double *w;
  double *x;

  r = new double[2];
  w = new double[order];
  x = new double[order];

  r[0] = - 1.0;
  r[1] = + 1.0;

  clenshaw_curtis_compute ( order, x, w );

  if ( prefix == "C++" || prefix == "c++" )
  {
    cout << "//\n";
    cout << "//  Weights W, abscissas X and range R\n";
    cout << "//  for a Clenshaw Curtis quadrature rule\n";
    cout << "//  ORDER = " << order << "\n";
    cout << "//\n";
    cout << "//  Standard rule:\n";
	cout << "//    Integral ( -1 <= x <= +1 ) f(x) dx\n";
	cout << "//    is to be approximated by\n";
	cout << "//    sum ( 1 <= I <= ORDER ) w(i) * f(x(i)).\n";
    cout << "//\n";

    for ( i = 0; i < order; i++ )
    {
      cout << "  w[" << i << "] = " 
           << setprecision(16) << w[i] << ";\n";
    }
    cout << "\n";
    for ( i = 0; i < order; i++ )
    {
      cout << "  x[" << i << "] = " 
           << setprecision(16) << x[i] << ";\n";
    }
    cout << "\n";
    for ( i = 0; i < 2; i++ )
    {
      cout << "  r[" << i << "] = " << r[i] << ";\n";
    }
   }
   else if ( prefix == "F77" || prefix == "f77" )
   {
    cout << "c\n";
    cout << "c  Weights W, abscissas X and range R\n";
    cout << "c  for a Clenshaw Curtis quadrature rule\n";
    cout << "c  ORDER = " << order << "\n";
    cout << "c\n";
	cout << "c  Standard rule:\n";
	cout << "c    Integral ( -1 <= x <= +1 ) f(x) dx\n";
	cout << "c    is to be approximated by\n";
	cout << "c    sum ( 1 <= I <= ORDER ) w(i) * f(x(i)).\n";
    cout << "c\n";

    for ( i = 0; i < order; i++ )
    {
      cout << "      w(" << i + 1 << ") = " 
           << setprecision(16) << w[i] << "\n";
    }
    cout << "\n";
    for ( i = 0; i < order; i++ )
    {
      cout << "      x(" << i + 1 << ") = " 
           << setprecision(16) << x[i] << "\n";
    }
    cout << "\n";
    for ( i = 0; i < 2; i++ )
    {
      cout << "      r(" << i + 1 << ") = " << r[i] << "\n";
    }
  }
  else if ( prefix == "F90" || prefix == "f90" )
  {  
    cout << "!\n";
    cout << "!  Weights W, abscissas X and range R\n";
    cout << "!  for a Clenshaw Curtis quadrature rule\n";
    cout << "!  ORDER = " << order << "\n";
    cout << "!\n";
	cout << "!  Standard rule:\n";
	cout << "!    Integral ( -1 <= x <= +1 ) f(x) dx\n";
	cout << "!    is to be approximated by\n";
	cout << "!    sum ( 1 <= I <= ORDER ) w(i) * f(x(i)).\n";
    cout << "!\n";

    for ( i = 0; i < order; i++ )
    {
      cout << "  w(" << i + 1 << ") = " 
           << setprecision(16) << w[i] << "\n";
    }
    cout << "\n";
    for ( i = 0; i < order; i++ )
    {
      cout << "  x(" << i + 1 << ") = " 
           << setprecision(16) << x[i] << "\n";
    }
    cout << "\n";
    for ( i = 0; i < 2; i++ )
    {
      cout << "  r(" << i + 1 << ") = " << r[i] << "\n";
    }
  }
  else if ( prefix == "MAT" || prefix == "mat" )
  {
    cout << "%\n";
    cout << "%  Weights W, abscissas X and range R\n";
    cout << "%  for a Clenshaw Curtis quadrature rule\n";
    cout << "%  ORDER = " << order << "\n";
    cout << "%\n";
	cout << "%  Standard rule:\n";
	cout << "%    Integral ( -1 <= x <= +1 ) f(x) dx\n";
	cout << "%    is to be approximated by\n";
	cout << "%    sum ( 1 <= I <= ORDER ) w(i) * f(x(i)).\n";
    cout << "%\n";

    for ( i = 0; i < order; i++ )
    {
      cout << "  w(" << i + 1 << ") = " 
           << setprecision(16) << w[i] << ";\n";
    }
    cout << "\n";
    for ( i = 0; i < order; i++ )
    {
      cout << "  x(" << i + 1 << ") = " 
           << setprecision(16) << x[i] << ";\n";
    }
    cout << "\n";
    for ( i = 0; i < 2; i++ )
    {
      cout << "  r(" << i + 1 << ") = " << r[i] << ";\n";
    }
  } 
  else
  {
    output_w = prefix + "_w.txt";
    output_x = prefix + "_x.txt";
    output_r = prefix + "_r.txt";

    cout << "\n";
    cout << "  Creating quadrature files.\n";
    cout << "\n";
    cout << "  Root file name is     \"" << prefix   << "\".\n";
    cout << "\n";
    cout << "  Weight file will be   \"" << output_w << "\".\n";
    cout << "  Abscissa file will be \"" << output_x << "\".\n";
    cout << "  Region file will be   \"" << output_r << "\".\n";
            
    r8mat_write ( output_w, 1, order, w );
    r8mat_write ( output_x, 1, order, x );
    r8mat_write ( output_r, 1, 2,     r );
  }
  
  delete [] r;
  delete [] w;
  delete [] x;

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

void r8mat_write ( string output_filename, int m, int n, double table[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_WRITE writes an R8MAT file with no header.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 June 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string OUTPUT_FILENAME, the output filename.
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, double TABLE[M*N], the table data.
//
{
  int i;
  int j;
  ofstream output;
//
//  Open the file.
//
  output.open ( output_filename.c_str ( ) );

  if ( !output )
  {
    cerr << "\n";
    cerr << "R8MAT_WRITE - Fatal error!\n";
    cerr << "  Could not open the output file.\n";
    return;
  }
//
//  Write the data.
//
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      output << "  " << setw(24) << setprecision(16) << table[i+j*m];
    }
    output << "\n";
  }
//
//  Close the file.
//
  output.close ( );

  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
}