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