```!   D01TDF Example Program Text
!   Mark 27.0 Release. NAG Copyright 2019.

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
```