Program e02affe

!     E02AFF Example Program Text

!     Mark 28.6 Release. NAG Copyright 2022.

!     .. 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 Procedures ..
      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