# include <stdlib.h>
# include <stdio.h>

# include "spline.h"
# include "test_interp_fun.h"

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

/******************************************************************************/

int main ( void )

/******************************************************************************/
/*
  Purpose:

    MAIN tests TEST_INTERP_FUN.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 February 2012

  Author:

    John Burkardt
*/
{
  timestamp ( );

  printf ( "\n" );
  printf ( "TEST_INTERP_FUN_TEST\n" );
  printf ( "  C version\n" );
  printf ( "  Test the TEST_INTERP_FUN library.\n" );

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

  printf ( "\n" );
  timestamp ( );

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

void test01 ( )

/******************************************************************************/
/*
  Purpose:

    TEST01 shows how P00_TITLE can be called.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 February 2012

  Author:

    John Burkardt
*/
{
  double a;
  double b;
  int prob;
  int prob_num;
  char title[100];

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

  prob_num = p00_prob_num ( );

  printf ( "\n" );
  printf ( "  Number of problems = %d\n", prob_num );

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

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

void test02 ( )

/******************************************************************************/
/*
  Purpose:

    TEST02 shows how P00_STORY can be called.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 February 2012

  Author:

    John Burkardt
*/
{
  int prob;
  int prob_num;

  printf ( "\n" );
  printf ( "TEST02\n" );
  printf ( "  P00_STORY prints the problem \"story\".\n" );

  prob_num = p00_prob_num ( );

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    printf ( "\n" );
    printf ( "  Problem %d\n", prob );
    p00_story ( prob );
  }

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

void test03 ( )

/******************************************************************************/
/*
  Purpose:

    TEST03 uses equally spaced polynomial interpolation.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 February 2012

  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;
  char title[100];
  int type;
  double x;
  double *xtab;
  double yapprox;
  double *ytab;
  double ytrue;

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

  prob_num = p00_prob_num ( );

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

    printf ( "\n" );
    printf ( "  Problem %d\n", prob );
    printf ( "  Problem TITLE = \"%s\".\n", title );
    printf ( "\n" );
    printf ( "     N   Max ||Error||\n" );
    printf ( "\n" );

    for ( n = 1; n <= max_tab; n = n + 4 )
    {
      xtab = ( double * ) malloc ( n * sizeof ( double ) );
      ytab = ( double * ) malloc ( n * sizeof ( double ) );
      diftab = ( double * ) malloc ( n * sizeof ( double ) );
/*
  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 ) );
      }
      printf ( "  %4d  %14g\n", n, error_max );
    }
    free ( diftab );
    free ( xtab );
    free ( ytab );
  }
  return;
}
/******************************************************************************/

void test04 ( )

/******************************************************************************/
/*
  Purpose:

    TEST04 uses Bernstein polynomial approximation on functional problems.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 February 2012

  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;
  char title[100];
  double *xdata;
  double xval;
  double *ydata;
  double ytrue;
  double yval;

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

  prob_num = p00_prob_num ( );

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

    printf ( "\n" );
    printf ( "  Problem %d\n", prob );
    printf ( "  Problem TITLE = \"%s\".\n", title );
    printf ( "\n" );
    printf ( "     N   Max ||Error||\n" );
    printf ( "\n" );

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

      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 ) );
      }
      printf ( "  %4d  %14g\n", ndata, error_max );

      free ( xdata );
      free ( ydata );
    }
  }
  return;
}
/******************************************************************************/

void test05 ( )

/******************************************************************************/
/*
  Purpose:

    TEST05 uses linear spline interpolation on all problems.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 February 2012

  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;
  char title[100];
  double *xdata;
  double xval;
  double *ydata;
  double ypval;
  double ytrue;
  double yval;

  printf ( "\n" );
  printf ( "TEST05\n" );
  printf ( "  Linear spline interpolation.\n" );

  prob_num = p00_prob_num ( );

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

    printf ( "\n" );
    printf ( "  Problem %d\n", prob );
    printf ( "  Problem TITLE = \"%s\".\n", title );
    printf ( "\n" );
    printf ( "     N   Max ||Error||\n" );
    printf ( "\n" );

    for ( ndata = 2; ndata <= max_data; ndata = ndata + 4 )
    {
      xdata = ( double * ) malloc ( ndata * sizeof ( double ) );
      ydata = ( double * ) malloc ( ndata * sizeof ( double ) );
/*
  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 ) );
      }

      printf ( "  %4d  %14g\n", ndata, error_max );

      free ( xdata );
      free ( ydata );
    }
  }
  return;
}
/******************************************************************************/

void test06 ( )

/******************************************************************************/
/*
  Purpose:

    TEST06 uses Overhauser spline interpolation on all problems.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 February 2012

  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;
  char title[100];
  double *xdata;
  double xval;
  double *ydata;
  double yval;

  printf ( "\n" );
  printf ( "TEST06\n" );
  printf ( "  Overhauser spline interpolation.\n" );

  prob_num = p00_prob_num ( );

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

    p00_lim ( prob, &a, &b );

    ndata = 11;
    xdata = ( double * ) malloc ( ndata * sizeof ( double ) );
    ydata = ( double * ) malloc ( ndata * sizeof ( double ) );

    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] );
    }
    printf ( "\n" );
    printf ( "  Problem %d\n", prob );
    printf ( "  Problem TITLE = \"%s\".\n", title );
    printf ( "\n" );
    printf ( "  X   Y\n" );
    printf ( "\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 = ' ';
        }
        printf ( "  %c  %14g  %14g\n", mark, xval, yval );
      }
    }
    free ( xdata );
    free ( ydata );
  }

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

void test07 ( )

/******************************************************************************/
/*
  Purpose:

    TEST07 uses cubic spline interpolation on all problems.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 February 2012

  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;
  char title[100];
  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;

  printf ( "\n" );
  printf ( "TEST07\n" );
  printf ( "  Cubic spline interpolation.\n" );

  prob_num = p00_prob_num ( );

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

    ndata = 11;
    xdata = ( double * ) malloc ( ndata * sizeof ( double ) );
    ydata = ( double * ) malloc ( ndata * sizeof ( double ) );
    ypp = ( double * ) malloc ( ndata * sizeof ( double ) );

    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 );

    printf ( "\n" );
    printf ( "  Problem %d\n", prob );
    printf ( "  Problem TITLE = \"%s\".\n", title );
    printf ( "\n" );
    printf ( "    X   Y\n" );
    printf ( "\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 = ' ';
        }
        printf ( "  %c  %14g  %14g\n", mark, xval, yval );
      }

    }
    free ( xdata );
    free ( ydata );
    free ( ypp );
  }
  return;
}
/******************************************************************************/

void test08 ( )

/******************************************************************************/
/*
  Purpose:

    TEST08 uses B spline approximation on all problems.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 February 2012

  Author:

    John Burkardt
*/
{
  double a;
  double b;
  int i;
  int j;
  int jhi;
  int jmax;
  char mark;
  int ndata;
  int prob;
  int prob_num;
  char title[100];
  double *xdata;
  double xval;
  double *ydata;
  double yval;

  printf ( "\n" );
  printf ( "TEST08\n" );
  printf ( "  B spline approximation.\n" );

  prob_num = p00_prob_num ( );

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

    ndata = 11;
    xdata = ( double * ) malloc ( ndata * sizeof ( double ) );
    ydata = ( double * ) malloc ( ndata * sizeof ( double ) );

    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] );
    }

    printf ( "\n" );
    printf ( "  Problem %d\n", prob );
    printf ( "  Problem TITLE = \"%s\".\n", title );
    printf ( "\n" );
    printf ( "       X        Y\n" );
    printf ( "\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 = ' ';
        }
        printf ( "  %c  %14g  %14g\n", mark, xval, yval );
      }

    }
    free ( xdata );
    free ( ydata );
  }
  return;
}