TD05AD

Evaluation of a transfer function G(jW) for a specified frequency

Purpose

  Given a complex valued rational function of frequency (transfer
  function) G(jW) this routine will calculate its complex value or
  its magnitude and phase for a specified frequency value.

      SUBROUTINE TD05AD( UNITF, OUTPUT, NP1, MP1, W, A, B, VALR, VALI,
     $                   INFO )
C     .. Scalar Arguments ..
      CHARACTER         OUTPUT, UNITF
      INTEGER           INFO, MP1, NP1
      DOUBLE PRECISION  VALI, VALR, W
C     .. Array Arguments ..
      DOUBLE PRECISION  A(*), B(*)

Mode Parameters

  UNITF   CHARACTER*1
          Indicates the choice of frequency unit as follows:
          = 'R':  Input frequency W in radians/second;
          = 'H':  Input frequency W in hertz.

  OUTPUT  CHARACTER*1
          Indicates the choice of co-ordinates for output as folows:
          = 'C':  Cartesian co-ordinates (output real and imaginary
                  parts of G(jW));
          = 'P':  Polar co-ordinates (output magnitude and phase
                  of G(jW)).

  NP1     (input) INTEGER
          The order of the denominator + 1, i.e. N + 1.  NP1 >= 1.

  MP1     (input) INTEGER
          The order of the numerator + 1, i.e. M + 1.  MP1 >= 1.

  W       (input) DOUBLE PRECISION
          The frequency value W for which the transfer function is
          to be evaluated.

  A       (input) DOUBLE PRECISION array, dimension (NP1)
          This array must contain the vector of denominator
          coefficients in ascending order of powers. That is, A(i)
          must contain the coefficient of (jW)**(i-1) for i = 1,
          2,...,NP1.

  B       (input) DOUBLE PRECISION array, dimension (MP1)
          This array must contain the vector of numerator
          coefficients in ascending order of powers. That is, B(i)
          must contain the coefficient of (jW)**(i-1) for i = 1,
          2,...,MP1.

  VALR    (output) DOUBLE PRECISION
          If OUTPUT = 'C', VALR contains the real part of G(jW).
          If OUTPUT = 'P', VALR contains the magnitude of G(jW)
                           in dBs.

  VALI    (output) DOUBLE PRECISION
          If OUTPUT = 'C', VALI contains the imaginary part of
                           G(jW).
          If OUTPUT = 'P', VALI contains the phase of G(jW) in
                           degrees.

  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value;
          = 1:  if the frequency value W is a pole of G(jW), or all
                the coefficients of the A polynomial are zero.

  By substituting the values of A, B and W in the following
  formula:

         B(1)+B(2)*(jW)+B(3)*(jW)**2+...+B(MP1)*(jW)**(MP1-1)
  G(jW) = ---------------------------------------------------.
         A(1)+A(2)*(jW)+A(3)*(jW)**2+...+A(NP1)*(jW)**(NP1-1)

  None.

  The algorithm requires 0(N+M) operations.

  None

Program Text

*     TD05AD EXAMPLE PROGRAM TEXT.
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NP1MAX, MP1MAX
      PARAMETER        ( NP1MAX = 20, MP1MAX = 20 )
*     .. Local Scalars ..
      DOUBLE PRECISION VALI, VALR, W
      INTEGER          I, INFO, MP1, NP1
      CHARACTER*1      UNITF, OUTPUT
*     .. Local Arrays ..
      DOUBLE PRECISION A(NP1MAX), B(MP1MAX)
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         TD05AD
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) NP1, MP1, W, UNITF, OUTPUT
      IF ( NP1.LE.0 .OR. NP1.GT.NP1MAX ) THEN
         WRITE ( NOUT, FMT = 99995 ) NP1
      ELSE
         READ ( NIN, FMT = * ) ( A(I), I = 1,NP1 )
         IF ( MP1.LE.0 .OR. MP1.GT.MP1MAX ) THEN
            WRITE ( NOUT, FMT = 99994 ) MP1
         ELSE
            READ ( NIN, FMT = * ) ( B(I), I = 1,MP1 )
*           Find the real and imaginary parts of G(jW), where
*           W = 1.0 radian.
            CALL TD05AD( UNITF, OUTPUT, NP1, MP1, W, A, B, VALR, VALI,
     $                   INFO )
*
            IF ( INFO.NE.0 ) THEN
               WRITE ( NOUT, FMT = 99998 ) INFO
            ELSE
               IF ( LSAME( OUTPUT, 'C' ) ) THEN
                  WRITE ( NOUT, FMT = 99997 ) VALR, VALI
               ELSE
                  WRITE ( NOUT, FMT = 99996 ) VALR, VALI
               END IF
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' TD05AD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TD05AD = ',I2)
99997 FORMAT (' Complex value of G(jW) = ',F8.4,1X,F8.4,'*j')
99996 FORMAT (' Magnitude of G(jW) = ',F8.4,' dBs, Phase of G(jW) = ',
     $       F8.4,' degrees ')
99995 FORMAT (/' NP1 is out of range.',/' NP1 = ',I5)
99994 FORMAT (/' MP1 is out of range.',/' MP1 = ',I5)
      END

 TD05AD EXAMPLE PROGRAM DATA
   6     4     1.0     R     C
   1.0   1.0   0.0   0.0   2.0   1.0
   6.0   2.0   3.0   1.0

 TD05AD EXAMPLE PROGRAM RESULTS

 Complex value of G(jW) =   0.8462  -0.2308*j

Return to index