program main

!*****************************************************************************80
!
!! MAIN is the main program for TEST_APPROX_PRB.
!
!  Discussion:
!
!    TEST_APPROX_PRB calls the TEST_APPROX tests.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 June 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp (  )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST_APPROX_PRB'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Demonstrate the TEST_APPROX'
  write ( *, '(a)' ) '  test interpolation and approximation functions.'

  call test001
  call test002
  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
  call test07
  call test08
  call test09

  call test10
  call test11
  call test12

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST_APPROX_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine test001

!*****************************************************************************80
!
!! TEST001 shows how P00_TITLE can be called.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    19 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  integer prob
  integer prob_num
  character ( len = 80 ) title
  integer type

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST001'
  write ( *, '(a)' ) '  Demonstrate some of the bookkeeping routines.'
  write ( *, '(a)' ) '  GET_PROB_NUM returns the number of problems.'
  write ( *, '(a)' ) '  P00_TITLE returns the problem title.'
  write ( *, '(a)' ) '  P00_TYPE returns the problem type.'
  write ( *, '(a)' ) '  P00_LIMIT returns the problem limits.'

  call get_prob_num ( prob_num )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Number of problems = ', prob_num

  do prob = 1, prob_num

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Problem ', prob
    call p00_title ( prob, title )
    write ( *, '(a)' ) '  Problem TITLE = "' // trim ( title ) // '".'

    call p00_type ( prob, type )
    write ( *, '(a,i8)' ) '  Problem TYPE = ', type

    if ( type == 0 ) then

      call p00_lim ( prob, a, b )
      write ( *, '(a,g14.6)' ) '  Problem lower limit A = ', a
      write ( *, '(a,g14.6)' ) '  Problem upper limit B = ', b

    end if

  end do

  return
end
subroutine test002

!*****************************************************************************80
!
!! TEST002 shows how P00_STORY can be called.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    19 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer prob
  integer prob_num

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST002'
  write ( *, '(a)' ) '  P00_STORY prints the problem "story".'

  call get_prob_num ( prob_num )

  do prob = 1, prob_num

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Problem ', prob

    call p00_story ( prob )

  end do

  return
end
subroutine test01

!*****************************************************************************80
!
!! TEST01 uses equally spaced polynomial interpolation on functional problems.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    19 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: max_tab = 8

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) diftab(max_tab)
  integer i
  integer imax
  integer ntab
  integer prob
  integer prob_num
  real ( kind = 8 ) p00_fun
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ) x
  real ( kind = 8 ) xtab(max_tab)
  real ( kind = 8 ) yapprox
  real ( kind = 8 ) ytab(max_tab)
  real ( kind = 8 ) ytrue

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  Equally spaced polynomial interpolation.'

  call get_prob_num ( prob_num )

  do prob = 1, prob_num

    call p00_title ( prob, title )
    call p00_type ( prob, type )

    if ( type /= 0 ) then

      write ( *, '(a)' ) ' '
      write ( *, '(a,i8)' ) '  Skipped problem ', prob
      write ( *, '(a)' ) '  Wrong type.'

    else

      call p00_lim ( prob, a, b )

      write ( *, '(a)' ) ' '
      write ( *, '(a,i8)' ) '  Problem ', prob
      write ( *, '(2x,a)' ) trim ( title )
      write ( *, '(a)' ) ' '

      do ntab = 2, max_tab

        write ( *, '(a)' ) ' '
        write ( *, '(a,i8)' ) '  Interpolating polynomial order = ', ntab
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) &
          '       X            Y(X)        Y_Interp        Error'
        write ( *, '(a)' ) ' '
!
!  Evaluate the function at N equally spaced points.
!
        do i = 1, ntab

          if ( ntab == 1 ) then
            xtab(i) = 0.5D+00 * ( a + b )
          else
            xtab(i) = ( real ( ntab - i,     kind = 8 ) * a &
                      + real (        i - 1, kind = 8 ) * b &
                    ) / real ( ntab     - 1, kind = 8 )
          end if

          ytab(i) = p00_fun ( prob, xtab(i) )

        end do
!
!  Construct the interpolating polynomial via finite differences.
!
        call data_to_dif ( ntab, xtab, ytab, diftab )
!
!  Now examine the approximation error.
!
        imax = 10
        do i = 0, imax

          if ( imax == 0 ) then
            x = 0.5D+00 * ( a + b )
          else
            x = ( real ( imax - i, kind = 8 ) * a    &
                + real (        i, kind = 8 ) * b  ) &
                / real ( imax,     kind = 8 )
          end if
          
          ytrue = p00_fun ( prob, x )
          call dif_val ( ntab, xtab, diftab, x, yapprox )
 
          write ( *, '(2x,3g14.6,g10.2)' ) x, ytrue, yapprox, ytrue - yapprox

        end do

      end do

    end if

  end do

  return
end
subroutine test02

!*****************************************************************************80
!
!! TEST02 uses polynomial interpolation on data vector problems.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    19 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: max_tab = 12

  real ( kind = 8 ) diftab(max_tab)
  integer i
  integer j
  integer jhi
  character mark
  integer ntab
  integer ndata
  integer prob
  integer prob_num
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ) x
  real ( kind = 8 ) xdata(max_tab)
  real ( kind = 8 ) yapprox
  real ( kind = 8 ) ydata(max_tab)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  Polynomial interpolation to a vector of data.'

  call get_prob_num ( prob_num )

  do prob = 1, prob_num

    call p00_title ( prob, title )
    call p00_type ( prob, type )

    if ( type /= 1 ) then

      write ( *, '(a,i8)' ) '  Skipped problem ', prob
      write ( *, '(a)' ) '  Wrong type.'

    else

      write ( *, '(a)' ) ' '
      write ( *, '(a,i8)' ) '  Problem ', prob
      write ( *, '(a)' ) trim ( title )

      call p00_ndata ( prob, ndata )

      if ( max_tab < ndata ) then

        write ( *, '(a)' ) ' '
        write ( *, '(a,i8)' ) '  Skipped problem ', prob
        write ( *, '(a)' ) '  Too big.'

      else

        call p00_dat ( prob, xdata, ydata, ndata )

        ntab = ndata

        write ( *, '(a)' ) ' '
        write ( *, '(a,i8)' ) '  Interpolating polynomial order = ', ntab
        write ( *, '(a)' ) ' '
!
!  Construct the interpolating polynomial via finite differences.
!
        call data_to_dif ( ntab, xdata, ydata, diftab )
!
!  Print out the approximation, including midpoints of the intervals.
!
        do i = 1, ntab

          if ( i < ntab ) then
            jhi = 2
          else
            jhi = 1
          end if

          do j = 1, jhi

            if ( i < ntab ) then
              x = ( real ( jhi - j + 1, kind = 8 ) * xdata(i)     &
                  + real (       j - 1, kind = 8 ) * xdata(i+1) ) &
                  / real ( jhi,         kind = 8 ) 
            else
              x = xdata(ntab)
            end if

            if ( j == 1 ) then
              mark = '*'
            else
              mark = ' '
            end if

            call dif_val ( ntab, xdata, diftab, x, yapprox )

            write ( *, '(2x,a,2g14.6)' ) mark, x, yapprox

          end do

        end do

      end if

    end if

  end do

  return
end
subroutine test03

!*****************************************************************************80
!
!! TEST03 uses Bernstein polynomial approximation on functional problems.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    19 May 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: max_tab = 5

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  integer i
  integer imax
  integer n
  integer prob
  integer prob_num
  real ( kind = 8 ) p00_fun
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ) xdata(0:max_tab)
  real ( kind = 8 ) xval
  real ( kind = 8 ) ydata(0:max_tab)
  real ( kind = 8 ) ytrue
  real ( kind = 8 ) yval

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)') '  Bernstein polynomial approximation.'

  call get_prob_num ( prob_num )

  do prob = 1, prob_num

    call p00_title ( prob, title )
    call p00_type ( prob, type )

    if ( type /= 0 ) then

      write ( *, '(a,i8)' ) '  Skipped problem ', prob
      write ( *, '(a)' ) '  Wrong type.'

    else

      call p00_lim ( prob, a, b )

      write ( *, '(a)' ) ' '
      write ( *, '(a,i8)' ) '  Problem ', prob
      write ( *, '(a)' ) trim ( title )
      write ( *, '(a)' ) ' '

      do n = 2, max_tab

        write ( *, '(a)' ) ' '
        write ( *, '(a,i8)' ) '  Approximating polynomial degree = ', n
        write ( *, '(a)' ) ' '
!
!  Evaluate the function at N + 1 equally spaced points.
!
        do i = 0, n

          if ( n == 0 ) then
            xdata(i) = 0.5D+00 * ( a + b )
          else
            xdata(i) = ( real ( n - i, kind = 8 ) * a    &
                       + real (     i, kind = 8 ) * b  ) &
                       / real ( n,     kind = 8 )
          end if

          ydata(i) = p00_fun ( prob, xdata(i) )

        end do
!
!  Now examine the approximation error.
!
        imax = 10
        do i = 0, imax

          if ( imax == 0 ) then
            xval = 0.5D+00 * ( a + b )
          else
            xval = ( real ( imax - i, kind = 8 ) * a    &
                   + real (        i, kind = 8 ) * b  ) &
                   / real ( imax,     kind = 8 )
          end if
          
          ytrue = p00_fun ( prob, xval )

          call bp_approx ( n, a, b, ydata, xval, yval )
 
          write ( *, '(2x,3g14.6,g10.2)' ) xval, ytrue, yval, ytrue - yval

        end do

      end do

    end if

  end do

  return
end
subroutine test04

!*****************************************************************************80
!
!! TEST04 uses linear spline interpolation on all problems.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 June 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  integer i
  integer imax
  character mark
  integer ndata
  integer prob
  integer prob_num
  real ( kind = 8 ) p00_fun
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ), allocatable, dimension ( : ) :: xdata
  real ( kind = 8 ) xval
  real ( kind = 8 ), allocatable, dimension ( : ) :: ydata
  real ( kind = 8 ) ypval
  real ( kind = 8 ) yval

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  Linear spline interpolation.'

  call get_prob_num ( prob_num )

  do prob = 1, prob_num

    call p00_title ( prob, title )
    call p00_type ( prob, type )

    call p00_lim ( prob, a, b )

    if ( type == 0 ) then

      ndata = 11
      allocate ( xdata(1:ndata) )
      allocate ( ydata(1:ndata) )

      do i = 1, ndata
        xdata(i) = ( real ( ndata - i,     kind = 8 ) * a   &
                   + real (         i - 1, kind = 8 ) * b ) &
                   / real ( ndata     - 1, kind = 8 )
        ydata(i) = p00_fun ( prob, xdata(i) )
      end do

    else

      call p00_ndata ( prob, ndata )
      allocate ( xdata(1:ndata) )
      allocate ( ydata(1:ndata) )

      call p00_dat ( prob, xdata, ydata, ndata )

    end if

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Problem ', prob
    write ( *, '(2x,a)' ) trim ( title )
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '       X          Y          Y'' '
    write ( *, '(a)' ) ' '
!
!  Evaluate the interpolation function.
!
    imax = 2 * ndata - 1

    do i = 1, imax

      xval = ( real ( imax - i,     kind = 8 ) * a   &
             + real (        i - 1, kind = 8 ) * b ) &
             / real ( imax     - 1, kind = 8 )
          
      call spline_linear_val ( ndata, xdata, ydata, xval, yval, ypval )
 
      if ( mod ( i, 2 ) == 1 ) then
        mark = '*'
      else
        mark = ' '
      end if

      write ( *, '(2x,a,3g14.6)' ) mark, xval, yval, ypval

    end do

    deallocate ( xdata )
    deallocate ( ydata )

  end do

  return
end
subroutine test05

!*****************************************************************************80
!
!! TEST05 uses Overhauser spline interpolation on all problems.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 June 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: num_dim = 1

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  integer i
  integer j
  integer jhi
  integer jmax
  character mark
  integer ndata
  integer prob
  integer prob_num
  real ( kind = 8 ) p00_fun
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ), allocatable, dimension ( : ) :: xdata
  real ( kind = 8 ) xval
  real ( kind = 8 ), allocatable, dimension ( : ) :: ydata
  real ( kind = 8 ) yval

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  Overhauser spline interpolation.'

  call get_prob_num ( prob_num )

  do prob = 1, prob_num

    call p00_title ( prob, title )

    call p00_type ( prob, type )

    call p00_lim ( prob, a, b )

    if ( type == 0 ) then

      ndata = 11
      allocate ( xdata(1:ndata) )
      allocate ( ydata(1:ndata) )

      do i = 1, ndata
        xdata(i) = ( real ( ndata - i,     kind = 8 ) * a   &
                   + real (         i - 1, kind = 8 ) * b ) &
                   / real ( ndata     - 1, kind = 8 )
        ydata(i) = p00_fun ( prob, xdata(i) )
      end do

    else

      call p00_ndata ( prob, ndata )

      allocate ( xdata(1:ndata) )
      allocate ( ydata(1:ndata) )

      call p00_dat ( prob, xdata, ydata, ndata )

    end if

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Problem ', prob
    write ( *, '(a)' ) trim ( title )
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  X   Y'
    write ( *, '(a)' ) ' '
!
!  Evaluate the interpolation function.
!
    do i = 1, ndata - 1
 
      jmax = 3

      if ( i == ndata - 1 ) then
        jhi = jmax
      else
        jhi = jmax - 1
      end if

      do j = 1, jhi
        
        xval = ( real ( jmax - j,     kind = 8 ) * xdata(i)     &
               + real (        j - 1, kind = 8 ) * xdata(i+1) ) &
               / real ( jmax     - 1, kind = 8 )
          
        call spline_overhauser_val ( num_dim, ndata, xdata, ydata, xval, &
          yval )
 
        if ( j == 1 .or. j == 3 ) then
          mark = '*'
        else
          mark = ' '
        end if

        write ( *, '(2x,a,2g14.6)' ) mark, xval, yval

      end do

    end do

    deallocate ( xdata )
    deallocate ( ydata )

  end do

  return
end
subroutine test06

!*****************************************************************************80
!
!! TEST06 uses cubic spline interpolation on all problems.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 June 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  integer i
  integer ibcbeg
  integer ibcend
  integer j
  integer jhi
  integer jmax
  character mark
  integer ndata
  integer prob
  integer prob_num
  real ( kind = 8 ) p00_fun
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ), allocatable, dimension ( : ) :: xdata
  real ( kind = 8 ) xval
  real ( kind = 8 ) ybcbeg
  real ( kind = 8 ) ybcend
  real ( kind = 8 ), allocatable, dimension ( : ) :: ydata
  real ( kind = 8 ), allocatable, dimension ( : ) :: ypp
  real ( kind = 8 ) yppval
  real ( kind = 8 ) ypval
  real ( kind = 8 ) yval

  ibcbeg = 0
  ibcend = 0
  ybcbeg = 0.0D+00
  ybcend = 0.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  Cubic spline interpolation.'

  call get_prob_num ( prob_num )

  do prob = 1, prob_num

    call p00_title ( prob, title )
    call p00_type ( prob, type )

    call p00_lim ( prob, a, b )

    if ( type == 0 ) then

      ndata = 11
      allocate ( xdata(1:ndata) )
      allocate ( ydata(1:ndata) )
      allocate ( ypp(1:ndata) )

      do i = 1, ndata
        xdata(i) = ( real ( ndata - i,     kind = 8 ) * a   &
                   + real (         i - 1, kind = 8 ) * b ) &
                   / real ( ndata     - 1, kind = 8 )
        ydata(i) = p00_fun ( prob, xdata(i) )
      end do

    else

      call p00_ndata ( prob, ndata )
      allocate ( xdata(1:ndata) )
      allocate ( ydata(1:ndata) )
      allocate ( ypp(1:ndata) )

      call p00_dat ( prob, xdata, ydata, ndata )

    end if
!
!  Set up the interpolation function.
!
      call spline_cubic_set ( ndata, xdata, ydata, ibcbeg, ybcbeg, &
        ibcend, ybcend, ypp )

      write ( *, '(a)' ) ' '
      write ( *, '(a,i8)' ) '  Problem ', prob
      write ( *, '(2x,a)' ) trim ( title )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '    X   Y'
      write ( *, '(a)' ) ' '
!
!  Evaluate the interpolation function.
!
      do i = 1, ndata - 1
 
        jmax = 3

        if ( i == ndata - 1 ) then
          jhi = jmax
        else
          jhi = jmax - 1
        end if

        do j = 1, jhi
        
          xval = ( real ( jmax - j,     kind = 8 ) * xdata(i)     &
                 + real (        j - 1, kind = 8 ) * xdata(i+1) ) &
                 / real ( jmax     - 1, kind = 8 )
          
          call spline_cubic_val ( ndata, xdata, ydata, ypp, xval, yval, ypval, &
            yppval )
 
          if ( j == 1 .or. j == 3 ) then
            mark = '*'
          else
            mark = ' '
          end if

          write ( *, '(2x,a,2g14.6)' ) mark, xval, yval

        end do

      end do

    deallocate ( xdata )
    deallocate ( ydata )
    deallocate ( ypp )

  end do

  return
end
subroutine test07

!*****************************************************************************80
!
!! TEST07 plots Overhauser spline interpolants.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 June 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: num_dim = 1
  integer, parameter :: max_point = 500

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  integer i
  integer ierror
  integer iunit
  integer j
  integer jhi
  integer jmax
  integer npoint
  integer ndata
  integer nx
  integer ny
  integer prob
  character ( len = 80 ) ps_file_name
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ), allocatable, dimension ( : ) :: xdata
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xmaxg
  real ( kind = 8 ) xming
  real ( kind = 8 ) xplot(max_point)
  real ( kind = 8 ) xval
  real ( kind = 8 ), allocatable, dimension ( : ) :: ydata
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin
  real ( kind = 8 ) ymaxg
  real ( kind = 8 ) yming
  real ( kind = 8 ) yplot(max_point)
  real ( kind = 8 ) yrange
  real ( kind = 8 ) yval

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  Plot an Overhauser spline interpolant.'

  prob = 9

  call p00_title ( prob, title )

  call p00_type ( prob, type )

  call p00_lim ( prob, a, b )

  call p00_ndata ( prob, ndata )

  allocate ( xdata(1:ndata) )
  allocate ( ydata(1:ndata) )

  call p00_dat ( prob, xdata, ydata, ndata )

  call r8vec2_print ( ndata, xdata, ydata, '  Data points:' )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) 'Problem ', prob
  write ( *, '(a)' ) trim ( title )

  iunit = 1
  ps_file_name = 'test07.ps'

  call ps_file_open ( ps_file_name, iunit, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST07'
    write ( *, '(a,i8)' ) '  File creation error ', ierror
    stop
  end if
!
!  Evaluate the interpolation function.
!
  npoint = 0

  do i = 1, ndata - 1
 
    jmax = 7

    if ( i == ndata - 1 ) then
      jhi = jmax
    else
      jhi = jmax - 1
    end if

    do j = 1, jhi
        
      xval = ( real ( jmax - j,     kind = 8 ) * xdata(i)     &
             + real (        j - 1, kind = 8 ) * xdata(i+1) ) &
             / real ( jmax     - 1, kind = 8 )
          
      call spline_overhauser_val ( num_dim, ndata, xdata, ydata, xval, yval )
 
      npoint = npoint + 1
      xplot(npoint) = xval
      yplot(npoint) = yval

    end do

  end do

  xmin = min ( minval ( xplot(1:npoint) ), minval ( xdata(1:ndata) ) )
  xmax = max ( maxval ( xplot(1:npoint) ), maxval ( xdata(1:ndata) ) )

  ymin = min ( minval ( yplot(1:npoint) ), minval ( ydata(1:ndata) ) )
  ymax = max ( maxval ( yplot(1:npoint) ), maxval ( ydata(1:ndata) ) )

  yrange = ymax - ymin
  ymin = ymin - 0.1D+00 * yrange
  ymax = ymax + 0.1D+00 * yrange

  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  XMIN, XMAX = ', xmin, xmax
  write ( *, '(a,2g14.6)' ) '  YMIN, YMAX = ', ymin, ymax

  call r8vec2_print ( npoint, xplot, yplot, '  Plot points:' )

  call ps_file_head ( ps_file_name )
  write ( *, *) 'DEBUG: About to call PS_PAGE_HEAD'
  call ps_page_head ( xmin, ymin, xmax, ymax )
  write ( *, *) 'DEBUG: About to call PS_LINE_OPEN'
  call ps_line_open ( npoint, xplot, yplot )

  do i = 1, ndata
    call ps_mark_circle ( xdata(i), ydata(i) )
  end do

  xming = xmin
  xmaxg = xmax
  nx = 21
  yming = ymin
  ymaxg = ymax
  ny = 21

  call ps_grid_cartesian ( xming, xmaxg, nx, yming, ymaxg, ny )

  call ps_page_tail ( )

  call ps_file_tail ( )

  call ps_file_close ( iunit )

  deallocate ( xdata )
  deallocate ( ydata )

  return
end
subroutine test08

!*****************************************************************************80
!
!! TEST08 plots cubic spline interpolants.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 June 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: max_point = 500

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  integer i
  integer ibcbeg
  integer ibcend
  integer ierror
  integer iunit
  integer j
  integer jhi
  integer jmax
  integer npoint
  integer ndata
  integer nx
  integer ny
  integer prob
  character ( len = 80 ) ps_file_name
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ), allocatable, dimension ( : ) :: xdata
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xmaxg
  real ( kind = 8 ) xming
  real ( kind = 8 ) xplot(max_point)
  real ( kind = 8 ) xval
  real ( kind = 8 ) ybcbeg
  real ( kind = 8 ) ybcend
  real ( kind = 8 ), allocatable, dimension ( : ) :: ydata
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin
  real ( kind = 8 ) ymaxg
  real ( kind = 8 ) yming
  real ( kind = 8 ) yplot(max_point)
  real ( kind = 8 ), allocatable, dimension ( : ) :: ypp
  real ( kind = 8 ) yppval
  real ( kind = 8 ) ypval
  real ( kind = 8 ) yrange
  real ( kind = 8 ) yval

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  Plot a cubic spline interpolant.'

  prob = 9

  call p00_title ( prob, title )

  call p00_type ( prob, type )

  call p00_lim ( prob, a, b )

  call p00_ndata ( prob, ndata )

  allocate ( xdata(1:ndata) )
  allocate ( ydata(1:ndata) )
  allocate ( ypp(1:ndata) )

  call p00_dat ( prob, xdata, ydata, ndata )
!
!  Set up the interpolation function.
!
  ibcbeg = 0
  ibcend = 0
  ybcbeg = 0.0D+00
  ybcend = 0.0D+00

  call spline_cubic_set ( ndata, xdata, ydata, ibcbeg, ybcbeg, ibcend, &
    ybcend, ypp )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) 'Problem ', prob
  write ( *, '(a)' ) trim ( title )

  iunit = 1
  ps_file_name = 'test08.ps'

  call ps_file_open ( ps_file_name, iunit, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST08'
    write ( *, '(a,i8)' ) '  File creation error ', ierror
    stop
  end if
!
!  Evaluate the interpolation function.
!
  npoint = 0
  xmin = a
  xmax = b
  ymin = 0.0D+00
  ymax = 0.0D+00

  do i = 1, ndata - 1
 
    jmax = 7

    if ( i == ndata - 1 ) then
      jhi = jmax
    else
      jhi = jmax - 1
    end if

    do j = 1, jhi
        
      xval = ( real ( jmax - j,     kind = 8 ) * xdata(i)     &
             + real (        j - 1, kind = 8 ) * xdata(i+1) ) & 
             / real ( jmax     - 1, kind = 8 )
          
      call spline_cubic_val ( ndata, xdata, ydata, ypp, xval, yval, ypval, &
        yppval )
 
      npoint = npoint + 1
      xplot(npoint) = xval
      yplot(npoint) = yval

      if ( npoint == 1 ) then
        ymin = yplot(npoint)
        ymax = yplot(npoint)
      else if ( j == 1 .or. j == jmax ) then
        ymin = min ( ymin, yplot(npoint) )
        ymax = max ( ymax, yplot(npoint) )
      end if

    end do

  end do

  call r8vec2_print ( npoint, xplot, yplot, '  Plot points:' )

  yrange = ymax - ymin
  ymin = ymin - 0.1D+00 * yrange
  ymax = ymax + 0.1D+00 * yrange

  call ps_file_head ( ps_file_name )

  call ps_page_head ( xmin, ymin, xmax, ymax )

  call ps_line_open ( npoint, xplot, yplot )

  do i = 1, ndata
    call ps_mark_circle ( xdata(i), ydata(i) )
  end do

  xming = xmin
  xmaxg = xmax
  nx = 21
  yming = ymin
  ymaxg = ymax
  ny = 21

  call ps_grid_cartesian ( xming, xmaxg, nx, yming, ymaxg, ny )

  call ps_page_tail ( )

  call ps_file_tail ( )

  call ps_file_close (iunit )

  deallocate ( xdata )
  deallocate ( ydata )
  deallocate ( ypp )

  return
end
subroutine test09

!*****************************************************************************80
!
!! TEST09 uses B spline approximation on all problems.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 June 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  integer i
  integer j
  integer jhi
  integer jmax
  character mark
  integer ndata
  integer prob
  integer prob_num
  real ( kind = 8 ) p00_fun
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ), allocatable, dimension ( : ) :: xdata
  real ( kind = 8 ) xval
  real ( kind = 8 ), allocatable, dimension ( : ) :: ydata
  real ( kind = 8 ) yval

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  B spline approximation.'

  call get_prob_num ( prob_num )

  do prob = 1, prob_num

    call p00_title ( prob, title )
    call p00_type ( prob, type )

    call p00_lim ( prob, a, b )

    if ( type == 0 ) then

      ndata = 11
      allocate ( xdata(1:ndata) )
      allocate ( ydata(1:ndata) )

      do i = 1, ndata
        xdata(i) = ( real ( ndata - i,     kind = 8 ) * a   &
                   + real (         i - 1, kind = 8 ) * b ) &
                   / real ( ndata     - 1, kind = 8 )
        ydata(i) = p00_fun ( prob, xdata(i) )
      end do

    else

      call p00_ndata ( prob, ndata )
      allocate ( xdata(1:ndata) )
      allocate ( ydata(1:ndata) )

      call p00_dat ( prob, xdata, ydata, ndata )

    end if

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Problem ', prob
    write ( *, '(2x,a)' ) trim ( title )
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '       X        Y'
    write ( *, '(a)' ) ' '
!
!  Evaluate the interpolation function.
!
    do i = 1, ndata - 1
 
      jmax = 3

      if ( i == ndata - 1 ) then
        jhi = jmax
      else
        jhi = jmax - 1
      end if

      do j = 1, jhi
        
        xval = ( real ( jmax - j,     kind = 8 ) * xdata(i)     &
               + real (        j - 1, kind = 8 ) * xdata(i+1) ) &
               / real ( jmax     - 1, kind = 8 )
          
        call spline_b_val ( ndata, xdata, ydata, xval, yval )
 
        if ( j == 1 .or. j == 3 ) then
          mark = '*'
        else
          mark = ' '
        end if

        write ( *, '(2x,a,2g14.6)' ) mark, xval, yval

      end do

    end do

    deallocate ( xdata )
    deallocate ( ydata )

  end do

  return
end
subroutine test10

!*****************************************************************************80
!
!! TEST10 plots B spline approximants.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 June 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: max_point = 500

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  integer i
  integer ierror
  integer iunit
  integer j
  integer jhi
  integer jmax
  integer npoint
  integer ndata
  integer nx
  integer ny
  integer prob
  character ( len = 80 ) ps_file_name
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ), allocatable, dimension ( : ) :: xdata
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xmaxg
  real ( kind = 8 ) xming
  real ( kind = 8 ) xplot(max_point)
  real ( kind = 8 ) xval
  real ( kind = 8 ), allocatable, dimension ( : ) :: ydata
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin
  real ( kind = 8 ) ymaxg
  real ( kind = 8 ) yming
  real ( kind = 8 ) yplot(max_point)
  real ( kind = 8 ) yrange
  real ( kind = 8 ) yval

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  Plot a B spline approximant.'

  prob = 9

  call p00_title ( prob, title )

  call p00_type ( prob, type )

  call p00_lim ( prob, a, b )

  call p00_ndata ( prob, ndata )

  allocate ( xdata(1:ndata) )
  allocate ( ydata(1:ndata) )

  call p00_dat ( prob, xdata, ydata, ndata )

  call r8vec2_print ( ndata, xdata, ydata, '  Data points:' )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Problem ', prob
  write ( *, '(2x,a)' ) trim ( title )

  iunit = 1
  ps_file_name = 'test10.ps'

  call ps_file_open ( ps_file_name, iunit, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST10'
    write ( *, '(a,i8)' ) '  File creation error ', ierror
    stop
  end if
!
!  Evaluate the interpolation function.
!
  npoint = 0
  xmin = a
  xmax = b
  ymin = 0.0D+00
  ymax = 0.0D+00

  do i = 1, ndata - 1
 
    jmax = 7

    if ( i == ndata - 1 ) then
      jhi = jmax
    else
      jhi = jmax - 1
    end if

    do j = 1, jhi
        
      xval = ( real ( jmax - j,     kind = 8 ) * xdata(i)     &
             + real (        j - 1, kind = 8 ) * xdata(i+1) ) &
             / real ( jmax     - 1, kind = 8 )
          
      call spline_b_val ( ndata, xdata, ydata, xval, yval )
 
      npoint = npoint + 1
      xplot(npoint) = xval
      yplot(npoint) = yval

      if ( npoint == 1 ) then
        ymin = yplot(npoint)
        ymax = yplot(npoint)
      else if ( j == 1 .or. j == jmax ) then
        ymin = min ( ymin, yplot(npoint) )
        ymax = max ( ymax, yplot(npoint) )
      end if

    end do

  end do

  call r8vec2_print ( npoint, xplot, yplot, '  Plot points:' )

  yrange = ymax - ymin
  ymin = ymin - 0.1D+00 * yrange
  ymax = ymax + 0.1D+00 * yrange

  call ps_file_head ( ps_file_name )

  call ps_page_head ( xmin, ymin, xmax, ymax )

  call ps_line_open ( npoint, xplot, yplot )

  do i = 1, ndata
    call ps_mark_circle ( xdata(i), ydata(i) )
  end do

  xming = xmin
  xmaxg = xmax
  nx = 21
  yming = ymin
  ymaxg = ymax
  ny = 21

  call ps_grid_cartesian ( xming, xmaxg, nx, yming, ymaxg, ny )

  call ps_page_tail ( )

  call ps_file_tail ( )

  call ps_file_close ( iunit )

  deallocate ( xdata )
  deallocate ( ydata )

  return
end
subroutine test11

!*****************************************************************************80
!
!! TEST11 plots beta spline approximants.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 June 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: max_point = 500

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) beta1
  real ( kind = 8 ) beta2
  integer i
  integer ierror
  integer iunit
  integer j
  integer jhi
  integer jmax
  integer npoint
  integer ndata
  integer nx
  integer ny
  integer prob
  character ( len = 80 ) ps_file_name
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ), allocatable, dimension ( : ) :: xdata
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xmaxg
  real ( kind = 8 ) xming
  real ( kind = 8 ) xplot(max_point)
  real ( kind = 8 ) xval
  real ( kind = 8 ), allocatable, dimension ( : ) :: ydata
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin
  real ( kind = 8 ) ymaxg
  real ( kind = 8 ) yming
  real ( kind = 8 ) yplot(max_point)
  real ( kind = 8 ) yrange
  real ( kind = 8 ) yval

  beta1 = 100.0D+00
  beta2 = 0.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  Plot a beta spline approximant.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  BETA1 = ', beta1
  write ( *, '(a,g14.6)' ) '  BETA2 = ', beta2

  prob = 9

  call p00_title ( prob, title )

  call p00_type ( prob, type )

  call p00_lim ( prob, a, b )

  call p00_ndata ( prob, ndata )

  allocate ( xdata(1:ndata) )
  allocate ( ydata(1:ndata) )

  call p00_dat ( prob, xdata, ydata, ndata )

  call r8vec2_print ( ndata, xdata, ydata, '  Data points:' )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Problem ', prob
  write ( *, '(2x,a)' ) trim ( title )

  iunit = 1
  ps_file_name = 'test11.ps'

  call ps_file_open ( ps_file_name, iunit, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST11'
    write ( *, '(a,i8)' ) '  File creation error ', ierror
    stop
  end if
!
!  Evaluate the interpolation function.
!
  npoint = 0
  xmin = a
  xmax = b
  ymin = 0.0D+00
  ymax = 0.0D+00

  do i = 1, ndata - 1
 
    jmax = 7

    if ( i == ndata - 1 ) then
      jhi = jmax
    else
      jhi = jmax - 1
    end if

    do j = 1, jhi
        
      xval = ( real ( jmax - j,     kind = 8 ) * xdata(i)     &
             + real (        j - 1, kind = 8 ) * xdata(i+1) ) &
             / real ( jmax     - 1, kind = 8 )
          
      call spline_beta_val ( beta1, beta2, ndata, xdata, ydata, xval, yval )
 
      npoint = npoint + 1
      xplot(npoint) = xval
      yplot(npoint) = yval

      if ( npoint == 1 ) then
        ymin = yplot(npoint)
        ymax = yplot(npoint)
      else if ( j == 1 .or. j == jmax ) then
        ymin = min ( ymin, yplot(npoint) )
        ymax = max ( ymax, yplot(npoint) )
      end if

    end do

  end do

  call r8vec2_print ( npoint, xplot, yplot, '  Plot points:' )

  yrange = ymax - ymin
  ymin = ymin - 0.1D+00 * yrange
  ymax = ymax + 0.1D+00 * yrange

  call ps_file_head ( ps_file_name )

  call ps_page_head ( xmin, ymin, xmax, ymax )

  call ps_line_open ( npoint, xplot, yplot )

  do i = 1, ndata
    call ps_mark_circle ( xdata(i), ydata(i) )
  end do

  xming = xmin
  xmaxg = xmax
  nx = 21
  yming = ymin
  ymaxg = ymax
  ny = 21

  call ps_grid_cartesian ( xming, xmaxg, nx, yming, ymaxg, ny )

  call ps_page_tail ( )

  call ps_file_tail ( )

  call ps_file_close ( iunit )

  deallocate ( xdata )
  deallocate ( ydata )

  return
end
subroutine test12

!*****************************************************************************80
!
!! TEST12 plots cubic spline interpolants.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    10 June 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: plot_num = 101

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  integer i
  integer ibcbeg
  integer ibcend
  integer ierror
  integer iunit
  integer j
  integer jhi
  integer jmax
  integer npoint
  integer ndata
  integer nx
  integer ny
  integer prob
  character ( len = 80 ) ps_file_name
  character ( len = 80 ) title
  integer type
  real ( kind = 8 ), allocatable, dimension ( : ) :: xdata
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xmaxg
  real ( kind = 8 ) xming
  real ( kind = 8 ) xplot(plot_num)
  real ( kind = 8 ) xval
  real ( kind = 8 ) ybcbeg
  real ( kind = 8 ) ybcend
  real ( kind = 8 ), allocatable, dimension ( : ) :: ydata
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin
  real ( kind = 8 ) ymaxg
  real ( kind = 8 ) yming
  real ( kind = 8 ) yplot(plot_num)
  real ( kind = 8 ), allocatable, dimension ( : ) :: ypp
  real ( kind = 8 ) yppval
  real ( kind = 8 ) ypval
  real ( kind = 8 ) yrange
  real ( kind = 8 ) yval

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  Plot a cubic spline interpolant.'

  prob = 10

  call p00_title ( prob, title )

  call p00_type ( prob, type )

  call p00_lim ( prob, a, b )
!
!  TEMPORARILY MODIFY B.  THINGS ARE BLOWING UP.
!
  b = 275.0

  call p00_ndata ( prob, ndata )

  allocate ( xdata(1:ndata) )
  allocate ( ydata(1:ndata) )
  allocate ( ypp(1:ndata) )

  call p00_dat ( prob, xdata, ydata, ndata )
!
!  Set up the interpolation function.
!
  ibcbeg = 0
  ibcend = 0
  ybcbeg = 0.0D+00
  ybcend = 0.0D+00

  call spline_cubic_set ( ndata, xdata, ydata, ibcbeg, ybcbeg, ibcend, &
    ybcend, ypp )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Problem ', prob
  write ( *, '(2x,a)' ) trim ( title )

  iunit = 1
  ps_file_name = 'test12.ps'

  call ps_file_open ( ps_file_name, iunit, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST12'
    write ( *, '(a,i8)' ) '  File creation error ', ierror
    stop
  end if
!
!  Evaluate the interpolation function.
!
  npoint = 0
  xmin = a
  xmax = b
  ymin =  huge ( ymin )
  ymax = -huge ( ymax )

  do j = 1, plot_num
        
    xval = ( real ( plot_num - j,     kind = 8 ) * xmin    &
           + real (            j - 1, kind = 8 ) * xmax ) & 
           / real ( plot_num     - 1, kind = 8 )
          
    call spline_cubic_val ( ndata, xdata, ydata, ypp, xval, yval, ypval, &
      yppval )
 
    xplot(j) = xval
    yplot(j) = yval

    ymin = min ( ymin, yplot(j) )
    ymax = max ( ymax, yplot(j) )

  end do

  call r8vec2_print ( plot_num, xplot, yplot, '  Plot points:' )

  yrange = ymax - ymin
  ymin = ymin - 0.1D+00 * yrange
  ymax = ymax + 0.1D+00 * yrange

  call ps_file_head ( ps_file_name )

  call ps_page_head ( xmin, ymin, xmax, ymax )

  call ps_line_open ( plot_num, xplot, yplot )

  do i = 1, ndata
    call ps_mark_circle ( xdata(i), ydata(i) )
  end do

  xming = xmin
  xmaxg = xmax
  nx = 21
  yming = ymin
  ymaxg = ymax
  ny = 21

  call ps_grid_cartesian ( xming, xmaxg, nx, yming, ymaxg, ny )

  call ps_page_tail ( )

  call ps_file_tail ( )

  call ps_file_close (iunit )

  deallocate ( xdata )
  deallocate ( ydata )
  deallocate ( ypp )

  return
end