subroutine chebyshev_coefficients ( a, b, n, f, c )

!*****************************************************************************80
!
!! chebyshev_coefficients() determines Chebyshev interpolation coefficients.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    14 September 2011
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Roger Broucke,
!    Algorithm 446:
!    Ten Subroutines for the Manipulation of Chebyshev Series,
!    Communications of the ACM,
!    Volume 16, Number 4, April 1973, pages 254-256.
!
!  Parameters:
!
!    Input, real ( kind = rk ) A, B, the domain of definition.
!
!    Input, integer N, the order of the interpolant.
!
!    Input, real ( kind = rk ), external :: F ( X ), an external function.
!
!    Output, real ( kind = rk ) C(N), the Chebyshev coefficients.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) a
  real ( kind = rk ) angle
  real ( kind = rk ) b
  real ( kind = rk ) c(n)
  real ( kind = rk ), external :: f
  real ( kind = rk ) fx(n)
  integer i
  integer j
  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  real ( kind = rk ) x

  do i = 1, n
    angle = real ( 2 * i - 1, kind = rk ) * pi / real ( 2 * n, kind = rk )
    x = cos ( angle )
    x = 0.5D+00 * ( a + b ) + x * 0.5D+00 * ( b - a )
    fx(i) = f ( x );
  end do

  do i = 1, n
    c(i) = 0.0D+00
    do j = 1, n
      angle = real ( ( i - 1 ) * ( 2 * j - 1 ), kind = rk ) * pi &
        / real ( 2 * n, kind = rk )
      c(i) = c(i) + fx(j) * cos ( angle )
    end do
  end do

  c(1:n) = 2.0D+00 * c(1:n) / real ( n, kind = rk )

  return
end
subroutine chebyshev_interpolant ( a, b, n, c, m, x, cf )

!*****************************************************************************80
!
!! chebyshev_interpolant() evaluates a Chebyshev interpolant.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    14 September 2011
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Roger Broucke,
!    Algorithm 446:
!    Ten Subroutines for the Manipulation of Chebyshev Series,
!    Communications of the ACM,
!    Volume 16, Number 4, April 1973, pages 254-256.
!
!  Parameters:
!
!    Input, real ( kind = rk ) A, B, the domain of definition.
!
!    Input, integer N, the order of the polynomial.
!
!    Input, real ( kind = rk ) C(N), the Chebyshev coefficients.
!
!    Input, integer M, the number of points.
!
!    Input, real ( kind = rk ) X(M), the point at which the polynomial is
!    to be evaluated.
!
!    Output, real ( kind = rk ) CF(M), the value of the Chebyshev
!    polynomial at X.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m
  integer n

  real ( kind = rk ) a
  real ( kind = rk ) b
  real ( kind = rk ) c(n)
  real ( kind = rk ) cf(m)
  real ( kind = rk ) di
  real ( kind = rk ) dip1
  real ( kind = rk ) dip2
  integer i
  integer j
  real ( kind = rk ) x(m)
  real ( kind = rk ) y

  do j = 1, m

    dip1 = 0.0D+00
    di = 0.0D+00
    y = ( 2.0D+00 * x(j) - a  - b ) / ( b - a )

    do i = n, 2, -1
      dip2 = dip1
      dip1 = di
      di = 2.0D+00 * y * dip1 - dip2 + c(i)
    end do

    cf(j) = y * di - dip1 + 0.5D+00 * c(1)

  end do

  return
end
subroutine chebyshev_zeros ( n, x )

!*****************************************************************************80
!
!! chebyshev_zeros() returns zeroes of the Chebyshev polynomial T(N,X).
!
!  Discussion:
!
!    We produce the Chebyshev zeros in ascending order.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    14 September 2011
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the order of the polynomial.
!
!    Output, real ( kind = rk ) X(N), the zeroes of T(N)(X).
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) angle
  integer i
  real ( kind = rk ), parameter :: pi = 3.141592653589793D+00
  real ( kind = rk ) x(n)

  do i = 1, n
    angle = real ( 2 * ( n - i ) + 1, kind = rk ) * pi / real ( 2 * n, kind = rk )
    x(i) = cos ( angle )
  end do

  return
end