# include <cstdlib>
# include <iostream>
# include <iomanip>

using namespace std;

# include "spline.hpp"
# include "test_interp_fun.hpp"

int main ( );
void test01 ( );
void test02 ( );
void test03 ( );
void test04 ( );
void test05 ( );
void test06 ( );
void test07 ( );
void test08 ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN tests TEST_INTERP_FUN.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    25 August 2011
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );

  cout << "\n";
  cout << "TEST_INTERP_FUN_TEST\n";
  cout << "  C++ version\n";
  cout << "  Test the TEST_INTERP_FUN library.\n";

  test01 ( );
  test02 ( );
  test03 ( );
  test04 ( );
  test05 ( );
  test06 ( );
  test07 ( );
  test08 ( );
//
//  Terminate.
//
  cout << "\n";
  cout << "TEST_INTERP_FUN_TEST\n";
  cout << "  Normal end of execution.\n";

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

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

void test01 ( )

//****************************************************************************80
//
//  Purpose:
//
//    TEST01 shows how P00_TITLE can be called.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    25 August 2011
//
//  Author:
//
//    John Burkardt
//
{
  double a;
  double b;
  int prob;
  int prob_num;
  string title;

  cout << "\n";
  cout << "TEST01\n";
  cout << "  Demonstrate some of the bookkeeping routines.\n";
  cout << "  P00_PROB_NUM returns the number of problems.\n";
  cout << "  P00_TITLE returns the problem title.\n";
  cout << "  P00_LIMIT returns the problem limits.\n";

  prob_num = p00_prob_num ( );

  cout << "\n";
  cout << "  Number of problems = " << prob_num << "\n";

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    cout << "\n";
    cout << "  Problem " << prob << "\n";
    title = p00_title ( prob );
    cout << "  Problem TITLE = \"" << title << "\".\n";
    p00_lim ( prob, a, b );
    cout << "  Problem lower limit A = " << a << "\n";
    cout << "  Problem upper limit B = " << b << "\n";
  }

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

void test02 ( )

//****************************************************************************80
//
//  Purpose:
//
//    TEST02 shows how P00_STORY can be called.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    25 August 2011
//
//  Author:
//
//    John Burkardt
//
{
  int prob;
  int prob_num;

  cout << "\n";
  cout << "TEST02\n";
  cout << "  P00_STORY prints the problem \"story\".\n";

  prob_num = p00_prob_num ( );

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    cout << "\n";
    cout << "  Problem " << prob << "\n";
    p00_story ( prob );
  }

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

void test03 ( )

//****************************************************************************80
//
//  Purpose:
//
//    TEST03 uses equally spaced polynomial interpolation.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    25 August 2011
//
//  Author:
//
//    John Burkardt
//
{
  double a;
  double b;
  double *diftab;
  double error_max;
  int i;
  int imax;
  int max_tab = 21;
  int n;
  int prob;
  int prob_num;
  string title;
  int type;
  double x;
  double *xtab;
  double yapprox;
  double *ytab;
  double ytrue;

  cout << "\n";
  cout << "TEST03\n";
  cout << "  Equally spaced polynomial interpolation.\n";
  cout << "  Evaluate the function at N equally spaced points.\n";
  cout << "  Determine the N-1 degre polynomial interpolant.\n";
  cout << "  Estimate the maximum difference between the function\n";
  cout << "  and the interpolant.\n";

  prob_num = p00_prob_num ( );

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    title = p00_title ( prob );
    p00_lim ( prob, a, b );

    cout << "\n";
    cout << "  Problem " << prob << "\n";
    cout << "  Problem TITLE = \"" << title << "\".\n";
    cout << "\n";
    cout << "     N   Max ||Error||\n";
    cout << "\n";

    for ( n = 1; n <= max_tab; n = n + 4 )
    {
      xtab = new double[n];
      ytab = new double[n];
      diftab = new double[n];
//
//  Evaluate the function at N equally spaced points.
//
      for ( i = 0; i < n; i++ )
      {
        if ( n == 1 )
        {
          xtab[i] = 0.5 * ( a + b );
        }
        else
        {
          xtab[i] = ( ( double ) ( n - i - 1 ) * a 
                    + ( double ) (     i     ) * b 
                  ) / ( double ) ( n     - 1 );
        }
        ytab[i] = p00_fun ( prob, xtab[i] );
      }
//
//  Construct the interpolating polynomial via finite differences.
//
      data_to_dif ( n, xtab, ytab, diftab );
//
//  Now examine the approximation error.
//
      error_max = 0.0;
      imax = 100;

      for ( i = 0; i <= imax; i++ )
      {
        if ( imax == 0 )
        {
          x = 0.5 * ( a + b );
        }
        else
        {
          x = ( ( double ) ( imax - i ) * a    
              + ( double ) (        i ) * b  ) 
              / ( double ) ( imax     );
        }

        ytrue = p00_fun ( prob, x );
        yapprox = dif_val ( n, xtab, diftab, x );
        error_max = r8_max ( error_max, r8_abs ( ytrue - yapprox ) );
      }
      cout << "  " << setw(4) << n
           << "  " << setw(14) << error_max << "\n";
    }
    delete [] diftab;
    delete [] xtab;
    delete [] ytab;
  }
  return;
}
//****************************************************************************80

void test04 ( )

//****************************************************************************80
//
//  Purpose:
//
//    TEST04 uses Bernstein polynomial approximation on functional problems.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    25 August 2011
//
//  Author:
//
//    John Burkardt
//
{
  double a;
  double b;
  double error_max;
  int i;
  int imax;
  int max_tab = 21;
  int ndata;
  int prob;
  int prob_num;
  string title;
  double *xdata;
  double xval;
  double *ydata;
  double ytrue;
  double yval;

  cout << "\n";
  cout << "TEST04\n";
  cout << "  Bernstein polynomial approximation.\n";
  cout << "  Evaluate the function at N equally spaced points.\n";
  cout << "  Determine the N-1 degree Bernstein polynomial approximant.\n";
  cout << "  Estimate the maximum difference between the function\n";
  cout << "  and the approximant.\n";

  prob_num = p00_prob_num ( );

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    title = p00_title ( prob );
    p00_lim ( prob, a, b );

    cout << "\n";
    cout << "  Problem " << prob << "\n";
    cout << "  Problem TITLE = \"" << title << "\".\n";
    cout << "\n";
    cout << "     N   Max ||Error||\n";
    cout << "\n";

    for ( ndata = 1; ndata <= max_tab; ndata = ndata + 4 )
    {
//
//  Evaluate the function at NDATA equally spaced points.
//
      xdata = new double[ndata];
      ydata = new double[ndata];

      for ( i = 0; i < ndata; i++ )
      {
        if ( ndata == 1 )
        {
          xdata[i] = 0.5 * ( a + b );
        }
        else
        {
          xdata[i] = ( ( double ) ( ndata - i - 1 ) * a 
                     + ( double ) (         i     ) * b 
                   ) / ( double ) ( ndata     - 1 );
        }
        ydata[i] = p00_fun ( prob, xdata[i] );
      }
//
//  Now examine the approximation error.
//
      error_max = 0.0;
      imax = 100;

      for ( i = 0; i <= imax; i++ )
      {
        if ( imax == 0 )
        {
          xval = 0.5 * ( a + b );
        }
        else
        {
          xval = ( ( double ) ( imax - i ) * a    
                 + ( double ) (        i ) * b  ) 
                 / ( double ) ( imax     );
        }

        ytrue = p00_fun ( prob, xval );

        yval = bpab_approx ( ndata - 1, a, b, ydata, xval );

        error_max = r8_max ( error_max, r8_abs ( ytrue - yval ) );
      }
      cout << "  " << setw(4) << ndata
           << "  " << setw(14) << error_max << "\n";

      delete [] xdata;
      delete [] ydata;
    }
  }
  return;
}
//****************************************************************************80

void test05 ( )

//****************************************************************************80
//
//  Purpose:
//
//    TEST05 uses linear spline interpolation on all problems.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    25 August 2011
//
//  Author:
//
//    John Burkardt
//
{
  double a;
  double b;
  double error_max;
  int i;
  int imax;
  char mark;
  int max_data = 21;
  int ndata;
  int prob;
  int prob_num;
  string title;
  double *xdata;
  double xval;
  double *ydata;
  double ypval;
  double ytrue;
  double yval;

  cout << "\n";
  cout << "TEST05\n";
  cout << "  Linear spline interpolation.\n";

  prob_num = p00_prob_num ( );

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    title = p00_title ( prob );
    p00_lim ( prob, a, b );

    cout << "\n";
    cout << "  Problem " << prob << "\n";
    cout << "  Problem TITLE = \"" << title << "\".\n";
    cout << "\n";
    cout << "     N   Max ||Error||\n";
    cout << "\n";

    for ( ndata = 2; ndata <= max_data; ndata = ndata + 4 )
    {
      xdata = new double[ndata];
      ydata = new double[ndata];
//
//  Evaluate the function at NDATA equally spaced points.
//
      for ( i = 0; i < ndata; i++ )
      {
        if ( ndata == 1 )
        {
          xdata[i] = 0.5 * ( a + b );
        }
        else
        {
          xdata[i] = ( ( double ) ( ndata - i - 1 ) * a 
                     + ( double ) (         i     ) * b 
                   ) / ( double ) ( ndata     - 1 );
        }
        ydata[i] = p00_fun ( prob, xdata[i] );
      }
//
//  Evaluate the interpolation function.
//
      error_max = 0.0;
      imax = 100;

      for ( i = 0; i <= imax; i++ )
      {
        if ( imax == 0 )
        {
          xval = 0.5 * ( a + b );
        }
        else
        {
          xval = ( ( double ) ( imax - i ) * a    
                 + ( double ) (        i ) * b  ) 
                 / ( double ) ( imax     );
        }

        ytrue = p00_fun ( prob, xval );

        spline_linear_val ( ndata, xdata, ydata, xval, &yval, &ypval );

        ytrue = p00_fun ( prob, xval );

        error_max = r8_max ( error_max, r8_abs ( yval - ytrue ) );
      }

      cout << "  " << setw(4) << ndata
           << "  " << setw(14) << error_max << "\n";

      delete [] xdata;
      delete [] ydata;
    }
  }
  return;
}
//****************************************************************************80

void test06 ( )

//****************************************************************************80
//
//  Purpose:
//
//    TEST06 uses Overhauser spline interpolation on all problems.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    25 August 2011
//
//  Author:
//
//    John Burkardt
//
{
  double a;
  double b;
  int i;
  int j;
  int jhi;
  int jmax;
  char mark;
  int ndata;
  int num_dim = 1;
  int prob;
  int prob_num;
  string title;
  double *xdata;
  double xval;
  double *ydata;
  double yval;

  cout << "\n";
  cout << "TEST06\n";
  cout << "  Overhauser spline interpolation.\n";

  prob_num = p00_prob_num ( );

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    title = p00_title ( prob );

    p00_lim ( prob, a, b );

    ndata = 11;
    xdata = new double[ndata];
    ydata = new double[ndata];

    for ( i = 0; i < ndata; i++ )
    {
      if ( ndata == 1 )
      {
        xdata[i] = 0.5 * ( a + b );
      }
      else
      {
        xdata[i] = ( ( double ) ( ndata - i - 1 ) * a 
                   + ( double ) (         i     ) * b 
                 ) / ( double ) ( ndata     - 1 );
      }
      ydata[i] = p00_fun ( prob, xdata[i] );
    }
    cout << "\n";
    cout << "  Problem " << prob << "\n";
    cout << "  Problem TITLE = \"" << title << "\".\n";
    cout << "\n";
    cout << "  X   Y\n";
    cout << "\n";
//
//  Evaluate the interpolation function.
//
    for ( i = 0; i < ndata - 1; i++ )
    {
      jmax = 3;

      if ( i == ndata - 2 )
      {
        jhi = jmax;
      }
      else
      {
        jhi = jmax - 1;
      }

      for ( j = 1; j <= jhi; j++ )
      {
        xval = ( ( double ) ( jmax - j     ) * xdata[i]     
               + ( double ) (        j - 1 ) * xdata[i+1] )
               / ( double ) ( jmax     - 1 );

        spline_overhauser_val ( num_dim, ndata, xdata, ydata, xval, &yval );

        if ( j == 1 || j == 3 )
        {
          mark = '*';
        }
        else
        {
          mark = ' ';
        }
        cout << "  " << mark 
             << "  " << setw(14) << xval
             << "  " << setw(14) << yval << "\n";
      }
    }
    delete [] xdata;
    delete [] ydata;
  }

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

void test07 ( )

//****************************************************************************80
//
//  Purpose:
//
//    TEST07 uses cubic spline interpolation on all problems.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    25 August 2011
//
//  Author:
//
//    John Burkardt
//
{
  double a;
  double b;
  int i;
  int ibcbeg;
  int ibcend;
  int j;
  int jhi;
  int jmax;
  char mark;
  int ndata;
  int prob;
  int prob_num;
  string title;
  double *xdata;
  double xval;
  double ybcbeg;
  double ybcend;
  double *ydata;
  double *ypp;
  double yppval;
  double ypval;
  double yval;

  ibcbeg = 0;
  ibcend = 0;
  ybcbeg = 0.0;
  ybcend = 0.0;

  cout << "\n";
  cout << "TEST07\n";
  cout << "  Cubic spline interpolation.\n";

  prob_num = p00_prob_num ( );

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    title = p00_title ( prob );
    p00_lim ( prob, a, b );

    ndata = 11;
    xdata = new double[ndata];
    ydata = new double[ndata];
    ypp = new double[ndata];

    for ( i = 0; i < ndata; i++ )
    {
      if ( ndata == 1 )
      {
        xdata[i] = 0.5 * ( a + b );
      }
      else
      {
        xdata[i] = ( ( double ) ( ndata - i - 1 ) * a 
                   + ( double ) (         i     ) * b 
                 ) / ( double ) ( ndata     - 1 );
      }
      ydata[i] = p00_fun ( prob, xdata[i] );
    }
//
//  Set up the interpolation function.
//
    ypp = spline_cubic_set ( ndata, xdata, ydata, ibcbeg, ybcbeg, ibcend, ybcend );

    cout << "\n";
    cout << "  Problem " << prob << "\n";
    cout << "  Problem TITLE = \"" << title << "\".\n";
    cout << "\n";
    cout << "    X   Y\n";
    cout << "\n";
//
//  Evaluate the interpolation function.
//
    for ( i = 0; i < ndata - 1; i++ )
    {
      jmax = 3;

      if ( i == ndata - 2 )
      {
        jhi = jmax;
      }
      else
      {
        jhi = jmax - 1;
      }

      for ( j = 1; j <= jhi; j++ )
      {
        xval = ( ( double ) ( jmax - j     ) * xdata[i]
               + ( double ) (        j - 1 ) * xdata[i+1] )
               / ( double ) ( jmax     - 1 );

        yval = spline_cubic_val ( ndata, xdata, ydata, ypp, xval, &ypval, &yppval );

        if ( j == 1 || j == 3 )
        {
          mark = '*';
        }
        else
        {
          mark = ' ';
        }
        cout << "  " << mark 
             << "  " << setw(14) << xval
             << "  " << setw(14) << yval << "\n";
      }

    }
    delete [] xdata;
    delete [] ydata;
    delete [] ypp;
  }
  return;
}
//****************************************************************************80

void test08 ( )

//****************************************************************************80
//
//  Purpose:
//
//    TEST08 uses B spline approximation on all problems.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    25 August 2011
//
//  Author:
//
//    John Burkardt
//
{
  double a;
  double b;
  int i;
  int j;
  int jhi;
  int jmax;
  char mark;
  int ndata;
  int prob;
  int prob_num;
  string title;
  double *xdata;
  double xval;
  double *ydata;
  double yval;

  cout << "\n";
  cout << "TEST08\n";
  cout << "  B spline approximation.\n";

  prob_num = p00_prob_num ( );

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    title = p00_title ( prob );
    p00_lim ( prob, a, b );

    ndata = 11;
    xdata = new double[ndata];
    ydata = new double[ndata];

    for ( i = 0; i < ndata; i++ )
    {
      if ( ndata == 1 )
      {
        xdata[i] = 0.5 * ( a + b );
      }
      else
      {
        xdata[i] = ( ( double ) ( ndata - i - 1 ) * a 
                   + ( double ) (         i     ) * b 
                 ) / ( double ) ( ndata     - 1 );
      }
      ydata[i] = p00_fun ( prob, xdata[i] );
    }

    cout << "\n";
    cout << "  Problem " << prob << "\n";
    cout << "  Problem TITLE = \"" << title << "\".\n";
    cout << "\n";
    cout << "       X        Y\n";
    cout << "\n";
//
//  Evaluate the interpolation function.
//
    for ( i = 0; i < ndata - 1; i++ )
    {
      jmax = 3;

      if ( i == ndata - 2 )
      {
        jhi = jmax;
      }
      else
      {
        jhi = jmax - 1;
      }

      for ( j = 1; j <= jhi; j++ )
      {
        xval = ( ( double ) ( jmax - j     ) * xdata[i]
               + ( double ) (        j - 1 ) * xdata[i+1] ) 
               / ( double ) ( jmax     - 1 );

        yval = spline_b_val ( ndata, xdata, ydata, xval );

        if ( j == 1 || j == 3 )
        {
          mark = '*';
        }
        else
        {
          mark = ' ';
        }
        cout << "  " << mark 
             << "  " << setw(14) << xval
             << "  " << setw(14) << yval << "\n";
      }

    }
    delete [] xdata;
    delete [] ydata;
  }
  return;
}