subroutine chebyshev_mantle (T, kmax, radius)

c  This subroutine calculates chebyshev polynomial in the 
c  mantle (at a given radius).

c  July 31, 1999.
c  Last Modified:  July 31, 1999.

cccccccccc
c  Input:
c	kmax = maximum order of Chebyshev polynomials
c	radius = radius within mantle in kilometres
c  Output:
c	T = chebyshev polynomial at the given radius
cccccccccc


      integer kmax
c  maximum order of chebyshev polynomial

      real T(0:kmax)
c  Chebyshev polynomial at a given position

      real radius
c  radius = radius at which T is calculated.  This
c  	    needs to be in km.

      real rMOHO, rCMB
      parameter (rMOHO = 6346.6, rCMB = 3480.0)
c  rMOHO = radius at MOHO in PREM
c  rCMB = radius of CMB in PREM.

      real x
c  x = position between -1 and 1.

      integer k
c  do loop variable

c==========PROGRAM===============

      do k = 0, kmax
	 T(k) = 0.0  ! initial cleaning up
      end do  ! end k loop

c  first, check radius to be within the mantle.
      if ((radius.gt.rMOHO).or.(radius.lt.rCMB)) then
	 print*, '!!Error: input radius = ', radius, ' km (sub_chebyshev.f)'
	 print*, '     Need radius within the mantle. (sub_chebyshev.f)'
	 print*, 'Quitting sub_chebyshev.f ...'
	 go to 15
      end if  ! if ((radius.gt.rMOHO).or.(radius.lt.rCMB))

c  convert radius to x (i.e., between -1 and 1 value)
      x = (2.0*radius - rMOHO - rCMB)/(rMOHO - rCMB)
c      print*, 'x = ', x, ' (sub_chebyshev.f)'

      call chebyshev(T, kmax, x)

 15   return

      END

c******************************************************************

        subroutine chebyshev (T, kmax, x)

c  This subroutine calculates chebyshev polynomial at a 
c  given position.

c  July 31, 1999.
c  Last Modified:  July 31, 1999.

cccccccccc
c  Input:
c	kmax = maximum order of Chebyshev polynomials
c	x = value between -1 and 1
c  Output:
c	T = chebyshev polynomial at the given radius
cccccccccc

cccccccccc
c  Note:  Chebyshev Polynomial Recursion Relationship
c    This subroutine uses the recursion relationship of Chebyshev
c    polynomials (not satisfying the closure)
c	T_n(x) = 2*x*T_{n-1}(x) - T_{n-2}(x)
c  Reference:
c    Gradshteyn, I.S., & Ryzhik, I.M., 1980.
c    Table of Integrals, Series, and Products. (corrected and
c    enlarged edition)
c    p1032, Academic Press, Inc., London, England.
cccccccccc

cccccccccc
c  Note:  Chebyshev Polynomial Closure Relationship
c    Chebyshev polynomials have the relation
c	int_{-1}^{1} [T_n(x)]**2 dx = 1 - 1/(4n*n-1)
c				    = (4n*n-2)/(4n*n-1)
c  Reference:
c    Gradshteyn, I.S., & Ryzhik, I.M., 1980.
c    Table of Integrals, Series, and Products. (corrected and
c    enlarged edition)
c    p833, Academic Press, Inc., London, England.
cccccccccc

      integer kmax
c  maximum order of chebyshev polynomial to be made

      real T(0:kmax)
c  Chebyshev polynomial at a given position

      real x
c  x = position between -1 and 1.

      integer N_Tn
      parameter (N_Tn = 20)
c  This is the maximum order allowed in Tn.

      real Tn(0:N_Tn)
c  This is the chebyshev polynomial without orthonormalisation.

      integer k
c  do loop variable

c===========PROGRAM===============

      do k = 0, kmax
	 T(k) = 0.0  ! initial cleaning up
      end do  ! end k loop

      if (kmax.gt.N_Tn) then  ! need to change dimension of Tn
	 print*, '!!Errror: Need to change dimensions of Tn '//
     +			' in subroutine chebyshev (sub_chebyshev.f)'
	 print*, '     kmax = ', kmax, ' but N_Tn = ', N_Tn
	 print*, 'Quitting subroutine chebyshev ...'
	 go to 15
      end if  ! if (kmax.gt.N_Tn)

      do k = 0, N_Tn
	 Tn(k) = 0.0  ! cleaning up
      end do  ! end k loop

c  determining un-orthonormalised Chebyshev polynomials
      Tn(0) = 1.0
      Tn(1) = x
      do k = 2, kmax
	 Tn(k) = 2.0*x*Tn(k-1) - Tn(k-2)
      end do  ! end k loop

c  orthonormalising...
      do k = 0, kmax
	 T(k) = Sqrt( (4.0*Float(k*k)-1.0)/(4.0*Float(k*k)-2.0) )*Tn(k)
      end do  ! end k loop

 15   return

      END