Program e02rafe

!     E02RAF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: c02agf, e02raf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: l = 4, m = 4, nout = 6
      Integer, Parameter               :: ia = l + 1
      Integer, Parameter               :: ib = m + 1
      Integer, Parameter               :: ic = ia + ib - 1
      Integer, Parameter               :: jw = ib*(2*ib+3)
      Logical, Parameter               :: scale = .True.
!     .. Local Scalars ..
      Integer                          :: i, ifail
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: a(ia), b(ib), c(ic), dd(ia+ib),      &
                                          w(jw), work(2*(l+m+1)), z(2,l+m)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'E02RAF Example Program Results'

!     Power series coefficients in C

      c(1) = 1.0E0_nag_wp

      Do i = 1, ic - 1
        c(i+1) = c(i)/real(i,kind=nag_wp)
      End Do

      ifail = 0
      Call e02raf(ia,ib,c,ic,a,b,w,jw,ifail)

      Write (nout,*)
      Write (nout,*) 'The given series coefficients are'
      Write (nout,99999) c(1:ic)
      Write (nout,*)
      Write (nout,*) 'Numerator coefficients'
      Write (nout,99999) a(1:ia)
      Write (nout,*)
      Write (nout,*) 'Denominator coefficients'
      Write (nout,99999) b(1:ib)

!     Calculate zeros of the approximant using C02AGF
!     First need to reverse order of coefficients

      dd(ia:1:-1) = a(1:ia)

      ifail = 0
      Call c02agf(dd,l,scale,z,work,ifail)

      Write (nout,*)
      Write (nout,*) 'Zeros of approximant are at'
      Write (nout,*)
      Write (nout,*) '    Real part    Imag part'
      Write (nout,99998)(z(1,i),z(2,i),i=1,l)

!     Calculate poles of the approximant using C02AGF
!     Reverse order of coefficients

      dd(ib:1:-1) = b(1:ib)

      ifail = 0
      Call c02agf(dd,m,scale,z,work,ifail)

      Write (nout,*)
      Write (nout,*) 'Poles of approximant are at'
      Write (nout,*)
      Write (nout,*) '    Real part    Imag part'
      Write (nout,99998)(z(1,i),z(2,i),i=1,m)

99999 Format (1X,5E13.4)
99998 Format (1X,2E13.4)
    End Program e02rafe