      function chbval(x,hc,m)
!
! +--------------------------------------------------------------------+
!  Evaluates the Chebyshev polynomial of order m (with coeffs hc(0:m)) !
!  at point x (-1 <= x <= 1: mapping must be made before calling).     !
! +--------------------------------------------------------------------+
!
      implicit none

      real*8 :: chbval

      integer, intent(in) :: m
      real*8, intent(in) :: x
      real*8, DIMENSION(:), intent(in) :: hC(0:m)

      real*8 :: d, dd, xx, sv
      integer :: j

      d = 0.
      dd = 0.
      xx = x + x
      do  j=m,1,-1
         sv = d
         d = xx*d - dd + hC(j)
         dd = sv
      enddo

      chbval = x*d - dd + 0.5d0*hC(0)

      return
      end function chbval

! ==================================================================== !

      subroutine chbder(a,b,hc,hcd,n)
!
! +--------------------------------------------------------------------+
! ! Construct the array hcd(0:n) for the Chebyshev approximation       !
! ! of the function F(x), a <= x <= b, which is the derivative of      !
! ! the function f(x) having Chebyshev approximation with              !
! ! coefficients hc(0:n) over the same interval.                       !
! +--------------------------------------------------------------------+
!
      implicit none

      integer, intent(in) :: n
      real*8, dimension(:), intent(in) :: hc(0:n)
      real*8, dimension(:), intent(out) :: hcd(0:n)
      real*8, intent(in) ::  a, b
     
      real*8 :: zab
      integer :: j

      hcd(n) = 0.
      hcd(n-1) = 2*n*hc(n)
      if(n.ge.2)  then
         do  j=n-1,1,-1
            hcd(j-1) = hcd(j+1) + 2*j*hc(j)
         enddo
      endif

      zab = 2/(b-a)
      do  j=0,n-1
         hcd(j) = zab*hcd(j)
      enddo

      return
      end subroutine chbder

! ==================================================================== !

      subroutine chbfit(a,b,hc,n,func)
!
! +--------------------------------------------------------------------+
! ! Construct the array hc(0:n) for the Chebyshev approximation        !
! ! of the external function func(x), a <= x <= b, which is passed.    !
! +--------------------------------------------------------------------+
!
! dmc 21 Mar 2003 -- based loosely on Numerical Recipes "chebft",
! Numerical Recipes in Fortran 77, p. 186; "n" is not quite the same
! here.
!
! real*8 with fortran-90 automatic arrays & use r8sincos to speed trig 
! function evaluations.
!
      implicit NONE

      real*8, intent(in) :: a,b         ! range over which func(x) is defined
      integer :: n                      ! max degree n of Chebyshev polynomial
      real*8, intent(out) :: hc(0:n)    ! Chebyshev Coefficients (out)

      real*8 :: func                    ! func(x) (passed)
      external func

! local stuff

      real*8, parameter :: pi = 3.1415926535897931D+00
      real*8, parameter :: HALF = 0.5D0
      real*8, parameter :: TWO = 2.0D0
      integer j,k,np1
      real*8 bma,bpa,fac,y,f(0:n),csum

      real*8 costab(n,0:n),sintab(n),cosarg(0:n)

! executable code

      bma=(b-a)*HALF
      bpa=(b+a)*HALF

      np1 = n+1

      do k=0,n
         cosarg(k) = PI*(k+HALF)/np1

         !  evaluate algebraic sin-cos ladder; save cos part only

         call r8sincos(cosarg(k),n,sintab,costab(1:n,k))

         !  now costab(j,k) = cos(j*cosarg(k))

         y=costab(1,k)

         !  evaluate func at the necessary points (cf num. recipies).

         f(k) = func(y*bma+bpa)
      enddo

      fac = TWO/np1

      do j=0,n
         if(j.eq.0) then
            csum = sum(f)               ! sum all f array elements
         else
            csum = 0
            do k=0,n
               csum = csum + f(k)*costab(j,k)
            enddo
         endif

         hc(j)= fac*csum

      enddo

      return
      end
