* E02RAF Example Program Text * Mark 16 Revised. NAG Copyright 1993. * .. Parameters .. INTEGER L, M, IA, IB, IC, JW PARAMETER (L=4,M=4,IA=L+1,IB=M+1,IC=IA+IB-1,JW=IB*(2*IB+3)) INTEGER NOUT PARAMETER (NOUT=6) LOGICAL SCALE PARAMETER (SCALE=.TRUE.) * .. Local Scalars .. INTEGER I, IFAIL * .. Local Arrays .. DOUBLE PRECISION A(IA), B(IB), C(IC), DD(IA+IB), W(JW), + WORK(2*(L+M+1)), Z(2,L+M) * .. External Subroutines .. EXTERNAL C02AGF, E02RAF * .. Intrinsic Functions .. INTRINSIC DBLE * .. Executable Statements .. WRITE (NOUT,*) 'E02RAF Example Program Results' * Power series coefficients in C C(1) = 1.0D0 DO 20 I = 1, IC - 1 C(I+1) = C(I)/DBLE(I) 20 CONTINUE IFAIL = 1 * CALL E02RAF(IA,IB,C,IC,A,B,W,JW,IFAIL) * WRITE (NOUT,*) IF (IFAIL.NE.0) THEN WRITE (NOUT,99997) ' ** E02RAF returned with IFAIL = ', IFAIL GO TO 80 END IF * WRITE (NOUT,*) 'The given series coefficients are' WRITE (NOUT,99999) (C(I),I=1,IC) WRITE (NOUT,*) WRITE (NOUT,*) 'Numerator coefficients' WRITE (NOUT,99999) (A(I),I=1,IA) WRITE (NOUT,*) WRITE (NOUT,*) 'Denominator coefficients' WRITE (NOUT,99999) (B(I),I=1,IB) * Calculate zeros of the approximant using C02AGF * First need to reverse order of coefficients DO 40 I = 1, IA DD(IA-I+1) = A(I) 40 CONTINUE 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 DO 60 I = 1, IB DD(IB-I+1) = B(I) 60 CONTINUE 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) 80 CONTINUE * 99999 FORMAT (1X,5E13.4) 99998 FORMAT (1X,2E13.4) 99997 FORMAT (1X,A,I5) END