Program e02affe
! E02AFF Example Program Text
! Mark 30.1 Release. NAG Copyright 2024.
! .. 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