C	WESTERN MICHIGAN UNIVERSITY
C	PFRAC.F4 (FILE NAME ON LIBRARY DECTAPE)
C	PFRAC, 2.10.1 (CALLING NAME, SUBLST NO.)
C	PARTIAL FRACTION EXPANSION
C	THIS PROGRAM WAS ADAPTED BY B.G. HOUCHARD FROM "BASIC
C	 CIRCUIT THEORY" BY L.P. HUELSMAN.  CHANGES ARE:
C	 (1) PROG. IS IN CONVERSATIONAL MODE  (2)  DIMENSION
C	 EXPANDED  (3) SUBROUTINE ROOT WAS SUBSTITUTED BY
C	 POLRT.  (4) OUTPUT FORMATTING WAS CHANGED TO FIT
C	 INTO TTY: OUTPUT.
C	LIBRARY DECTAPE PROG. USED:  USAGE.MAC
C	APLIB PROGS. USED:  GETFOR
C	INTERNAL SUBR. USED:  PDIV, PFEXSD, DIFPL, CVALPL, POLRT
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C
CC    MAIN PROGRAM FOR FINDING PARTIAL FRACTION EXPANSION
C                M - DEGREE OF NUMERATOR POLYNOMIAL
C                N - DEGREE OF DENOMINATOR POLYNOMIAL
C                C - ARRAY OF NUMERATOR POLYNOMIAL COEFFICIENTS
C                B - ARRAY OF DENOMINATOR POLYNOMIAL COEFFICIENTS
C                P - COMPLEX ARRAY OF POLE LOCATIONS
C               R1 - COMPLEX ARRAY OF RESIDUES FOR 1ST-ORDER POLES
C               R2 - COMPLEX ARRAY OF CONSTANTS FOR 2ND-ORDER POLES
      DIMENSION C(36), B(36), G(36), A(36),NFT(16),ID(16)
      COMPLEX P(36), R1(36), R2(36)
	DOUBLE PRECISION NUM,DEN
	COMMON/IOBLK/IDLG,ICC,IOUT
	DATA NUM,DEN/'NUMERATOR','DENOMINATO'/
	DATA NUM1,DEN1/' ','R'/
	IDLG=-1
	ICC=-4
	IOUT=30
C	CALL USAGE('PFRAC')
C---------------NFT, ISTD ARE RETURNED.  OTHER ARGS. ARE INPUT.  
C---------------16=NO. OF FMT. WORDS FOR OBJ. TIME FORMAT
C---------------(1 LINE).  4 MEANS FORMAT NOT RESTRICTED.
1020	CALL GETFOR(IDLG,ICC,NFT,ISTD,16,4)
	IF (ISTD.NE.1) GO TO 102
	NFT(1)='(10E)'
	DO 1021 I=2,16
1021	NFT(I)=' '
102	WRITE(IDLG,510)
510	FORMAT(' ENTER HEADER'/)
	READ(ICC,511) ID
511	FORMAT(16A5)
	WRITE(IDLG,100)
100	FORMAT(' ENTER DEGREE OF NUMERATOR POLYNOMIAL FOLLOWED BY THAT
     1 OF THE DENOMINATOR'/' SEPARATED BY COMMA, MAXIMUM IS 35'/)
C
C    READ THE DATA FOR THE NUMERATOR AND DENOMINATOR POLYNOMIALS
      READ(ICC,2) M,N
2	FORMAT(2I)
	IF (((M.GT.0).AND.(M.LE.35)).AND.((N.GT.0).AND.(N.LE.35)))
     1 GO TO 50
	WRITE(IDLG,101)
101	FORMAT(' DEGREE OF POLYNOMIAL OUTSIDE ALLOWABLE RANGE'/
     1 ' MAXIMUM IS 35, TRY AGAIN'/)
	GO TO 102
50	WRITE(IDLG,51) NUM,NUM1
51	FORMAT(' ENTER COEFFICIENTS OF THE ',A10,A1,' POLYNOMIAL
     1 STARTING WITH CONSTANT')
	IF (ISTD.NE.1) WRITE(IDLG,512)
	IF (ISTD.EQ.1) WRITE(IDLG,513)
512	FORMAT('+TERM'/)
513	FORMAT('+TERM, 10 PER LINE SEPARATED BY COMMAS'/)
      READ(ICC,NFT) (C(I),I=1,M+1)
	WRITE(IDLG,51) DEN,DEN1
	IF (ISTD.NE.1) WRITE(IDLG,512)
	IF (ISTD.EQ.1) WRITE(IDLG,513)
      READ(ICC,NFT) (B(I),I=1,N+1)
C
C    CALL SUBROUTINE PDIV TO REMOVE THE INFINITE FREQUENCY TERMS
      CALL PDIV (M, C, N, B, KI, G, KR, A)
      IF (KI.LT.0) GO TO 13
	WRITE(IOUT,123)
123	FORMAT('1')
	DO 120 I=1,16
	IF (ID(I).NE.' ') GO TO 121
120	CONTINUE
	GO TO 122
121	WRITE(IOUT,124) ID
124	FORMAT('-',16A5)
122      WRITE(IOUT,12) (G(I),I=1,KI+1)
12	FORMAT('-INFINITE TERMS =',5E13.5)
C
C    CALL SUBROUTINE POLRT TO FIND THE POLE LOCATIONS
13	CALL POLRT(B,N,P)
C    CALL SUBROUTINE PFEXSD TO FIND THE RESIDUES
      CALL PFEXSD (KR, A, N, B, P, R1, R2)
	WRITE(IOUT,15)
15	FORMAT('-',12X,'POLE LOCATION'/)
	DO 16 I=1,N
16	WRITE(IOUT,17) I, P(I)
17	FORMAT(1X,I3,2X,2(E11.4,'+J',E11.4,2X))
	WRITE(IOUT,18)
18	FORMAT('-',23X,'R E S I D U E     O F'/11X,'FIRST DEGREE TERM',
     1 8X,'SECOND DEGREE TERM'//)
	DO 19 I=1,N
19	WRITE(IOUT, 17) I,R1(I),R2(I)
	WRITE(IDLG,200)
200	FORMAT('-END OF CALCULATION, TYPE 1 TO CONTINUE OR 0 TO
     1 TERMINATE'/)
	READ(ICC,2) LAST
	IF (LAST.GT.0) GO TO 1020
	CALL EXIT
	END
C---------------M, C, N, B ARE INPUT.  OTHER ARGS. ARE RETURNED.
      SUBROUTINE PDIV (M, C, N, B, KI, G, KR, A)
C    SUBROUTINE TO REMOVE SINGULARITIES AT INFINITY FROM A RATIO OF
C    POLYNOMIALS C(S)/B(S). THIS RATIO IS EXPRESSED AS G(S) + A(S)/B(S)
C    WHERE ALL POLYNOMIALS HAVE THE FORM A(1)+A(2)*S+A(3)*S**2+...
C    THE ORIGINAL POLYNOMIALS ARE ALL PRESERVED
C           M - DEGREE OF NUMERATOR POLYNOMIAL C(S)
C           C - ARRAY OF NUMERATOR POLYNOMIAL COEFFICIENTS
C           N - DEGREE OF DENOMINATOR POLYNOMIAL B(S)
C           B - ARRAY OF DENOMINATOR POLYNOMIAL COEFFICIENTS
C          KI - DEGREE OF QUOTIENT POLYNOMIAL G(S) (= M - N)
C           G - OUTPUT ARRAY OF QUOTIENT POLYNOMIAL COEFFICIENTS
C          KR - DEGREE OF REMAINDER POLYNOMIAL A(S) (= N - 1)
C           A - OUTPUT ARRAY OF REMAINDER POLYNOMIAL COEFFICIENTS
      DIMENSION C(36), B(36), G(36), A(36), D(36)
C
C    STORE THE COEFFICIENTS OF C(S) IN THE D ARRAY
      DO 2 I=1,10
2     D(I)=C(I)
      KR=M
      KI=M-N
      IF (KI.LT.0) GO TO 17
      KR=N-1
      MT=M
8     K=MT-N+1
      IF (K.LT.1) GO TO 19
C
C    COMPUTE THE KTH COEFFICIENT OF THE POLYNOMIAL G(S)
      G(K)=D(MT+1)/B(N+1)
      L=K
C
C    COMPUTE THE COEFFICIENTS OF THE REDUCED POLYNOMIAL D
      DO 14  I=1,N
      D(K)=D(K)-G(L)*B(I)
14    K=K+1
      MT=MT-1
      GO TO 8
C
C    IF THE NUMERATOR IS LESS IN DEGREE THAN THE DENOMINATOR
C    SET THE COEFFICIENTS OF THE G ARRAY TO ZERO
17    DO 18 I=1,10
18    G(I)=0.0
      KI=0
19    KRP=KR+1
C
C    STORE THE COEFFICIENTS OF THE REDUCED POLYNOMIAL D IN
C    THE A ARRAY AS THE REMAINDER POLYNOMIAL
      DO 21 I=1,KRP
21    A(I)=D(I)
      RETURN
      END
C---------------R1, R2 ARE RETURNED.  OTHER ARGS. ARE INPUT.
C--------------- P IS RETURNED FROM POLRT.
      SUBROUTINE PFEXSD (M, A, N, B, P, R1, R2)
C    SUBROUTINE FOR FINDING THE COEFFICIENTS IN THE PARTIAL
C    FRACTION EXPANSION FOR SIMPLE AND DOUBLE POLES OF A
C    RATIONAL FUNCTION F(S)=A(S)/B(S) WHERE THE POLYNOMIALS
C    HAVE THE FORM A(1) + A(2)*S + A(3)*S**2+ ...
C           M - DEGREE OF NUMERATOR POLYNOMIAL A(S)
C           A - ARRAY OF COEFFICIENTS OF POLYNOMIAL A(S)
C           N - DEGREE OF DENOMINATOR POLYNOMIAL B(S)
C           B - ARRAY OF COEFFICIENTS OF POLYNOMIAL B(S)
C           P - ARRAY OF POLES FOR WHICH RESIDUES ARE TO BE FOUND
C          R1 - RESIDUE OF A SIMPLE POLE, COEFFICIENT OF FIRST
C                        DEGREE TERM FOR A DOUBLE POLE
C          R2 - RESIDUE FOR A DOUBLE POLE (SET TO ZERO FOR A
C                        SIMPLE POLE)
C    THE COEFFICIENTS OF THE POLYNOMIALS A(S) AND B(S) ARE PRESERVED
      COMPLEX P(36), PA(36), R1(36), R2(36), CV, AV, CDV, CDDV, SA, SB
      DIMENSION A(36), B(36), C(36), CD(36), CDD(36)
	COMMON/IOBLK/IDLG,ICC,IOUT
C
C    CONSTRUCT THE POLYNOMIAL C(S) EQUAL TO B(S) AND FIND
C    ITS DERIVATIVE
      NP=N+1
      DO 3 I=1,NP
3     C(I)=B(I)
      CALL DIFPL (N, C, NM, CD)
C
C    PERFORM THE PARTIAL-FRACTION EXPANSION PROCEDURE N TIMES
      DO 49 K=1,N
C
C    TEST TO MAKE CERTAIN THAT S = P(K) IS A ZERO OF C(S)
17    CALL CVALPL (N, C, P(K), CV)
      IF (CABS(CV).GT.0.001) GO TO 42
      CALL CVALPL (M, A, P(K), AV)
C
C    TEST THE DERIVATIVE OF C(S) TO DETERMINE THE ORDER
C    OF THE ZERO AT S = P(K)
      CALL CVALPL (NM, CD, P(K), CDV)
      IF (CABS(CDV).LT.0.001) GO TO 25
C
C    IF THE ZERO OF C(S) IS SIMPLE, COMPUTE THE RESIDUE
C    R1, SET R2 TO ZERO, AND RETURN
      R1(K)=AV/CDV
      R2(K)=(0.0,0.0)
      GO TO 49
C
C    TEST TO MAKE CERTAIN THAT THE ZERO OF C(S) IS NOT
C    OF HIGHER DEGREE THAN SECOND
25    CALL DIFPL (NM, CD, NC, CDD)
      CALL CVALPL (NC, CDD, P(K), CDDV)
      IF (CABS(CDDV).LT.0.001) GO TO 45
C
C    IF THE ZERO OF C(S) IS DOUBLE, COMPUTE THE COEFFICIENTS R1 AND R2
C    FOR THE TERMS IN THE PARTIAL FRACTION EXPANSION
      R2(K)=2.*AV/CDDV
      PRK=REAL(P(K))
      PIK=AIMAG(P(K))
      SA=(0.,0.)
      IF (M.EQ.0) GO TO 34
	CALL POLRT(A,M,PA)
      DO 33 I=1,M
33    SA=SA+(1.,0.)/(P(K)-PA(I))
34    SB=(0.,0.)
      DO 39 I=1,N
      IF (ABS(PRK-REAL(P(I))).GT.0.001) GO TO 38
      IF (ABS(PIK-AIMAG(P(I))).LT.0.001) GO TO 39
38    SB=SB+(1.,0.)/(P(K)-P(I))
39    CONTINUE
      R1(K)=R2(K)*(SA-SB)
      GO TO 49
C
C    PRINT ANY ERROR MESSAGES THAT ARE REQUIRED
42    WRITE(IDLG,43)
43    FORMAT (/40H0FUNCTION F(S) DOES NOT HAVE A POLE AT P/)
      GO TO 47
45    WRITE(IDLG,46)
46    FORMAT (/42H0POLE OF F(S) IS OF GREATER THAN 2ND ORDER/)
47    WRITE(IDLG,48) P(K)
48    FORMAT (4H P =,E11.3,3H +J,E11.3/)
49    CONTINUE
      RETURN
      END
C---------------N, A ARE INPUT.  M, B ARE RETURNED.
      SUBROUTINE DIFPL (N, A, M, B)
C    SUBROUTINE FOR DIFFERENTIATING A POLYNOMIAL OF THE FORM
C    A(1) + A(2)*P + A(3)*P**2 + ... + A(N+1)*P**N
C    TO OBTAIN THE RESULTING DIFFERENTIATED POLYNOMIAL
C    B(1) + B(2)*P + B(3)*P**2 + ... + A(M+1)*P**M
C           N - DEGREE OF INPUT POLYNOMIAL
C           A - ARRAY OF COEFFICIENTS OF INPUT POLYNOMIAL
C           M - DEGREE OF DIFFERENTIATED POLYNOMIAL (=N-1)
C           B - OUTPUT ARRAY OF COEFFICIENTS OF DIFFERENTIATED POLYNOMIAL
      DIMENSION A(36), B(36)
      DO 3 I=1,N
3     B(I)=A(I+1)*I
      M=N-1
      RETURN
      END
C---------------VALUE IS RETURNED.  M, A, P ARE INPUT.
      SUBROUTINE CVALPL (M, A, P, VALUE)
C    SUBROUTINE FOR DETERMINING THE VALUE OF A POLYNOMIAL
C    A(1) + A(2)*P + A(3)*P**2 + ... + A(M+1)*P**M
C    FOR A COMPLEX VALUE OF ITS VARIABLE
C         M - DEGREE OF POLYNOMIAL
C         A - ARRAY OF POLYNOMIAL COEFFICIENTS (IN ASCENDING ORDER)
C         P - COMPLEX VALUE OF FREQUENCY VARIABLE
C    VALUE - OUTPUT COMPLEX VALUE OF THE POLYNOMIAL
      DIMENSION A(36)
      COMPLEX B(36), P, VALUE
      MC=M+1
C
C    STORE COEFFICIENTS AS COMPLEX VARIABLES
      DO 3 I=1,MC
3     B(I)=CMPLX(A(I), 0.)
C
C    STORE VALUE OF HIGHEST DEGREE COEFFICIENT
      VALUE=B(MC)
      IF (M.LE.0) RETURN
C
C    ACCUMULATE EFFECT OF OTHER COEFFICIENTS
      DO 8 I=1,M
      K=MC-I
8     VALUE=VALUE*P+B(K)
      RETURN
      END
C
C     ..................................................................
C
C        SUBROUTINE POLRT
C
C        PURPOSE
C           COMPUTES THE REAL AND COMPLEX ROOTS OF A REAL POLYNOMIAL
C
C        USAGE
C           CALL POLRT(XCOF,M,P,IER)
C
C        DESCRIPTION OF PARAMETERS
C           XCOF -VECTOR OF M+1 COEFFICIENTS OF THE POLYNOMIAL
C                 ORDERED FROM SMALLEST TO LARGEST POWER
C           COF  -WORKING VECTOR OF LENGTH M+1
C           M    -ORDER OF POLYNOMIAL
C           ROOTR-RESULTANT VECTOR OF LENGTH M CONTAINING REAL ROOTS
C                 OF THE POLYNOMIAL
C           ROOTI-RESULTANT VECTOR OF LENGTH M CONTAINING THE
C                 CORRESPONDING IMAGINARY ROOTS OF THE POLYNOMIAL
C           P    -VECTOR CONTAINING THE REAL AND IMAGINARY PARTS OF THE
C                 ROOTS
C
C        REMARKS
C           LIMITED TO 36TH ORDER POLYNOMIAL OR LESS.
C           FLOATING POINT OVERFLOW MAY OCCUR FOR HIGH ORDER
C           POLYNOMIALS BUT WILL NOT AFFECT THE ACCURACY OF THE RESULTS.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           NEWTON-RAPHSON ITERATIVE TECHNIQUE.  THE FINAL ITERATIONS
C           ON EACH ROOT ARE PERFORMED USING THE ORIGINAL POLYNOMIAL
C           RATHER THAN THE REDUCED POLYNOMIAL TO AVOID ACCUMULATED
C           ERRORS IN THE REDUCED POLYNOMIAL.
C
C     ..................................................................
C
C---------------XCOF, M ARE INPUT.  P IS RETURNED.
      SUBROUTINE POLRT(XCOF,M,P)
      DIMENSION XCOF(1),COF(36),ROOTR(36),ROOTI(36)
	COMPLEX P(36)
      DOUBLE PRECISION XO,YO,X,Y,XPR,YPR,UX,UY,V,YT,XT,U,XT2,YT2,SUMSQ,
     1 DX,DY,TEMP,ALPHA
	COMMON/IOBLK/IDLG,ICC,IOUT
C
C        ...............................................................
C
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C        STATEMENT WHICH FOLLOWS.
C
C     DOUBLE PRECISION XCOF,COF,ROOTR,ROOTI
C
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C        ROUTINE.
C        THE DOUBLE PRECISION VERSION MAY BE MODIFIED BY CHANGING THE
C        CONSTANT IN STATEMENT 78 TO 1.0D-12 AND IN STATEMENT 122 TO
C        1.0D-10.  THIS WILL PROVIDE HIGHER PRECISION RESULTS AT THE
C        COST OF EXECUTION TIME
C
C        ...............................................................
C
      IFIT=0
      N=M
   35 NX=N
      NXX=N+1
      N2=1
      KJ1 = N+1
      DO 40 L=1,KJ1
      MT=KJ1-L+1
   40 COF(MT)=XCOF(L)
C
C        SET INITIAL VALUES
C
   45 XO=.00500101
      YO=0.01000101
C
C        ZERO INITIAL VALUE COUNTER
C
      IN=0
   50 X=XO
C
C        INCREMENT INITIAL VALUES AND COUNTER
C
      XO=-10.0*YO
      YO=-10.0*X
C
C        SET X AND Y TO CURRENT VALUE
C
      X=XO
      Y=YO
      IN=IN+1
      GO TO 59
   55 IFIT=1
      XPR=X
      YPR=Y
C
C        EVALUATE POLYNOMIAL AND DERIVATIVES
C
   59 ICT=0
   60 UX=0.0
      UY=0.0
      V =0.0
      YT=0.0
      XT=1.0
      U=COF(N+1)
      IF(U) 65,130,65
   65 DO 70 I=1,N
      L =N-I+1
      TEMP=COF(L)
      XT2=X*XT-Y*YT
      YT2=X*YT+Y*XT
      U=U+TEMP*XT2
      V=V+TEMP*YT2
      FI=I
      UX=UX+FI*XT*TEMP
      UY=UY-FI*YT*TEMP
      XT=XT2
   70 YT=YT2
      SUMSQ=UX*UX+UY*UY
      IF(SUMSQ) 75,110,75
   75 DX=(V*UY-U*UX)/SUMSQ
      X=X+DX
      DY=-(U*UY+V*UX)/SUMSQ
      Y=Y+DY
   78 IF(DABS(DY)+DABS(DX)-1.0D-05) 100,80,80
C
C        STEP ITERATION COUNTER
C
   80 ICT=ICT+1
      IF(ICT-500) 60,85,85
   85 IF(IFIT)100,90,100
   90 IF(IN-5) 50,95,95
95	WRITE(IDLG,200)
200	FORMAT(' UNABLE TO DETERMINE ROOT WITH 500 ITERATIONS ON
     1 5 STARTING VALUES'/)
20	DO 16 I=1,M
16	P(I)=CMPLX(ROOTR(I),ROOTI(I))
	RETURN
  100 DO 105 L=1,NXX
      MT=KJ1-L+1
      TEMP=XCOF(MT)
      XCOF(MT)=COF(L)
  105 COF(L)=TEMP
      ITEMP=N
      N=NX
      NX=ITEMP
      IF(IFIT) 120,55,120
  110 IF(IFIT) 115,50,115
  115 X=XPR
      Y=YPR
  120 IFIT=0
  122 IF(DABS(Y)-1.0D-4*DABS(X)) 135,125,125
  125 ALPHA=X+X
      SUMSQ=X*X+Y*Y
      N=N-2
      GO TO 140
  130 X=0.0
      NX=NX-1
      NXX=NXX-1
  135 Y=0.0
      SUMSQ=0.0
      ALPHA=X
      N=N-1
  140 COF(2)=COF(2)+ALPHA*COF(1)
  145 DO 150 L=2,N
  150 COF(L+1)=COF(L+1)+ALPHA*COF(L)-SUMSQ*COF(L-1)
  155 ROOTI(N2)=Y
      ROOTR(N2)=X
      N2=N2+1
      IF(SUMSQ) 160,165,160
  160 Y=-Y
      SUMSQ=0.0
      GO TO 155
  165 IF(N) 20,20,45
      END

