PROGRAM e02affe

!      E02AFF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : e02aef, e02aff, nag_wp, x01aaf
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: fit, pi, piby2n
       INTEGER                         :: i, ifail, j, n, nplus1, r
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: a(:), f(:), xcap(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          real, sin
!      .. Executable Statements ..
       WRITE (nout,*) 'E02AFF Example Program Results'

!      Skip heading in data file
       READ (nin,*)

       READ (nin,*) n
       nplus1 = n + 1
       ALLOCATE (a(nplus1),f(nplus1),xcap(nplus1))

       piby2n = 0.5E0_nag_wp*x01aaf(pi)/real(n,kind=nag_wp)

       READ (nin,*) (f(r),r=1,nplus1)

       DO r = 1, nplus1
          i = r - 1

!         The following method of evaluating  XCAP = cos(PI*I/N)
!         ensures that the computed value has a small relative error
!         and, moreover, is bounded in modulus by unity for all
!         I = 0, 1, ..., N.  (It is assumed that the sine routine
!         produces a result with a small relative error for values
!         of the argument between  -PI/4  and  PI/4).

          IF (4*i<=n) THEN
             xcap(i+1) = 1.0E0_nag_wp - 2.0E0_nag_wp*sin(piby2n*real(i,kind= &
                nag_wp))**2
          ELSE IF (4*i>3*n) THEN
             xcap(i+1) = 2.0E0_nag_wp*sin(piby2n*real(n-i,kind=nag_wp))**2 - &
                1.0E0_nag_wp
          ELSE
             xcap(i+1) = sin(piby2n*real(n-2*i,kind=nag_wp))
          END IF

       END DO

       ifail = 0
       CALL e02aff(nplus1,f,a,ifail)

       WRITE (nout,*)
       WRITE (nout,*) '         Chebyshev'
       WRITE (nout,*) '  J   coefficient A(J)'
       WRITE (nout,99998) (j,a(j),j=1,nplus1)
       WRITE (nout,*)
       WRITE (nout,*) '  R    Abscissa   Ordinate      Fit'

       DO r = 1, nplus1

          ifail = 0
          CALL e02aef(nplus1,a,xcap(r),fit,ifail)

          WRITE (nout,99999) r, xcap(r), f(r), fit
       END DO

99999  FORMAT (1X,I3,3F11.4)
99998  FORMAT (1X,I3,F14.7)
    END PROGRAM e02affe