Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/pfrac.f4
There are no other files named pfrac.f4 in the archive.
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