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