!   D01TDF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2017.

    Module d01tdfe_mod

!     D01TDF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: basic_types
!     .. Parameters ..
      Integer, Parameter, Public       :: chebyshev1 = 2, chebyshev2 = 3,      &
                                          hermite = 6, jacobi = 4,             &
                                          laguerre = 5, legendre = 1
    Contains
      Subroutine basic_types(rulekind,alpha,beta,n,a,b,c,muzero)
!       This procedure supplies the coefficients of the three term
!       recurrence relationship for various classical orthogonal
!       polynomials.

!       .. Use Statements ..
        Use nag_library, Only: s14aaf, x01aaf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: alpha, beta
        Real (Kind=nag_wp), Intent (Out) :: muzero
        Integer, Intent (In)           :: n, rulekind
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: a(1:n), b(1:n), c(1:n)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: abl, mypi
        Integer                        :: i, ifail
!       .. Intrinsic Procedures ..
        Intrinsic                      :: real, sqrt
!       .. Executable Statements ..
        mypi = x01aaf(mypi)

        Select Case (rulekind)

        Case (legendre)
          muzero = 2.0_nag_wp
!         P(x) in [-1, 1), w(x) = 1.0
          Do i = 1, n
            a(i) = (real(2*i-1,kind=nag_wp))/real(i,kind=nag_wp)
            b(i) = 0.0_nag_wp
            c(i) = real(i-1,kind=nag_wp)/real(i,kind=nag_wp)
          End Do

        Case (chebyshev1)

          muzero = mypi
!         c(i)= (i-1)/i
!         T (x) in [-1,1], w(x) = (1-x**2)**(-0.5)
          Do i = 1, n
            a(i) = 2.0_nag_wp
            b(i) = 0.0_nag_wp
            c(i) = 1.0_nag_wp
          End Do
          If (n>0) Then
            a(1) = 1.0_nag_wp
          End If

        Case (chebyshev2)

          muzero = mypi/2.0_nag_wp
!         u(x) in [-1, 11, W(x) = (1-x**2)** 0.5;
          Do i = 1, n
            a(i) = 2.0_nag_wp
            b(i) = 0.0_nag_wp
            c(i) = 1.0_nag_wp
          End Do

        Case (jacobi)
          ifail = 0
          muzero = 2**(alpha+beta+1)*s14aaf(alpha+1,ifail)*                    &
            s14aaf(beta+1,ifail)/s14aaf(alpha+beta+2,ifail)
!         P(alpha,beta)(x) in [-1, 11, w(x) = (1-x)**alpha*(l+x)**beta
!         alpha> -1 and beta > -1
          If (n>0) Then
            a(1) = 0.5_nag_wp*(alpha+beta+2)
            b(1) = 0.5_nag_wp*(alpha-beta)
            c(1) = 0.0_nag_wp
          End If
          Do i = 2, n
            abl = 2.0_nag_wp*real(i,kind=nag_wp)*(real(i,kind=nag_wp)+alpha+   &
              beta)
            a(i) = (2.0_nag_wp*real(i,kind=nag_wp)+alpha+beta-1.0_nag_wp)*     &
              (2.0_nag_wp*real(i,kind=nag_wp)+alpha+beta)/abl
            abl = (2.0_nag_wp*real(i,kind=nag_wp)+alpha+beta-2.0_nag_wp)*abl
            b(i) = (2.0_nag_wp*real(i,kind=nag_wp)+alpha+beta-1.0_nag_wp)*     &
              (alpha**2-beta**2)/abl
            c(i) = 2.0_nag_wp*(real(i,kind=nag_wp)-1.0_nag_wp+alpha)*          &
              (real(i,kind=nag_wp)-1.0_nag_wp+beta)*                           &
              (2.0_nag_wp*real(i,kind=nag_wp)+alpha+beta)/abl
          End Do

        Case (laguerre)
          ifail = 0
          muzero = s14aaf(alpha+1.0_nag_wp,ifail)
!         L(alpha)(x) in [0, infinity), w(x) = exp(-x)*x**alpha,
!         alpha > -1
          Do i = 1, n
            a(i) = -1.0_nag_wp/real(i,kind=nag_wp)
            b(i) = (2.0_nag_wp*real(i,kind=nag_wp)-1.0_nag_wp+alpha)/          &
              real(i,kind=nag_wp)
            c(i) = (real(i,kind=nag_wp)-1.0_nag_wp+alpha)/real(i,kind=nag_wp)
          End Do

        Case (hermite)
          muzero = sqrt(mypi)
!         H(x) in (-infinity,+infinity), w(x) = exp(-x**2)
          Do i = 1, n
            a(i) = 2.0_nag_wp
            b(i) = 0.0_nag_wp
            c(i) = 2.0_nag_wp*real(i-1,kind=nag_wp)
          End Do

        End Select
      End Subroutine basic_types
    End Module d01tdfe_mod

    Program d01tdfe

!     .. Use Statements ..
      Use d01tdfe_mod, Only: basic_types, chebyshev1, chebyshev2, hermite,     &
                             jacobi, laguerre, legendre
      Use nag_library, Only: d01tdf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: n = 4, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: alpha, beta, muzero
      Integer                          :: ifail, j, rule
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: a(1:n), b(1:n), c(1:n), t(1:n),      &
                                          w(1:n)
!     .. Executable Statements ..

      Write (nout,*) 'D01TDF Example Program Results '
      Write (6,*)
      rule = legendre
      Do rule = 1, 6
        ifail = -1

        Select Case (rule)

        Case (legendre)
          Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
          Write (nout,*) 'Using the Gauss-Legendre Rule'
          Call d01tdf(n,a,b,c,muzero,w,t,ifail)
          If (ifail==0) Then
            Write (nout,99998)
            Write (nout,99997)(j,t(j),w(j),j=1,n)
            Write (6,*)
          End If

        Case (chebyshev1)
          Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
          Write (nout,*) 'Using the Chebyshev Rule 1'
          Call d01tdf(n,a,b,c,muzero,w,t,ifail)
          If (ifail==0) Then
            Write (nout,99998)
            Write (nout,99997)(j,t(j),w(j),j=1,n)
            Write (6,*)
          End If

        Case (chebyshev2)
          Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
          Write (nout,*) 'Using the Chebyshev Rule 2'
          Call d01tdf(n,a,b,c,muzero,w,t,ifail)
          If (ifail==0) Then
            Write (nout,99998)
            Write (nout,99997)(j,t(j),w(j),j=1,n)
            Write (6,*)
          End If

        Case (jacobi)
          alpha = 0.5_nag_wp
          beta = 0.5_nag_wp
          Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
          Write (nout,99999) alpha, beta
          Call d01tdf(n,a,b,c,muzero,w,t,ifail)
          If (ifail==0) Then
            Write (nout,99998)
            Write (nout,99997)(j,t(j),w(j),j=1,n)
            Write (6,*)
          End If

        Case (laguerre)
          Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
          Write (nout,*) 'Using the Laguerre Rule'
          Call d01tdf(n,a,b,c,muzero,w,t,ifail)
          If (ifail==0) Then
            Write (nout,99998)
            Write (nout,99997)(j,t(j),w(j),j=1,n)
            Write (6,*)
          End If

        Case (hermite)
          Call basic_types(rule,alpha,beta,n,a,b,c,muzero)
          Write (nout,*) 'Using the Hermite Rule'
          Call d01tdf(n,a,b,c,muzero,w,t,ifail)
          If (ifail==0) Then
            Write (nout,99998)
            Write (nout,99997)(j,t(j),w(j),j=1,n)
            Write (6,*)
          End If

        End Select

      End Do
99999 Format (' Using the Jacobi Rule: alpha = ',F10.5,'  beta = ',F10.5)

99998 Format (/,'   j    ','         Abscissa        ','          Weight')
99997 Format (I8,D25.15,D25.15)
    End Program d01tdfe