Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/bmd/bmdx68.for
There is 1 other file named bmdx68.for in the archive. Click here to see a list.
C BMDX68 - MULTIPLE TIME SERIES ANALYSIS MAY 17, 1965
C
COMMON UNIT
COMMON A,B,T,U,ISAVE,KSAVE,KSAVE2,N1,N2,N3,N4,NLINES,V,W,S,FREQ,
1NSERS,FMULT
DIMENSION A(10,10),B(10,10),P( 9, 9), KIN1(9),KOUT1(9),KSAVE(10),
XFMOUT(12),E(9,9),KSAVE2(10),FMT(120),U(10),V(5),S(5),W(5),G(9,9),
XT(10),ISAVE(9),CG(10,10),CP(10,10),JSAVE(10),JSAVE2(10)
EQUIVALENCE (PROB,IA1),(FIN,IA2),(RIN,IA3),(MULT,FMULT),
1(YES,IYES),(FNO,INO)
C.....PULLOUT EQUIVALENCE, DIMENSION, AND STATEMENT CARDS MARKED 'JERK',
C.....WHEN THE CARDS USING 'AIMAG' AND 'REAL' ARE RESTORED.
DIMENSION AX(2,10,10),BX(2,10,10)
EQUIVALENCE (A,AX),(B,BX)
DOUBLE PRECISION PROB,FIN,RIN,OUT,GOB,UNIT,FILL,CODE,AKRD
COMPLEX A,B,T,U
DATA PROB,FIN,RIN,OUT,YES,FNO/6HPROBLM,6HFINISH,6HINPUT ,6HOUTPUT,
X3HYES,3H NO/,PLUS,YMINUS/1H+,1H-/,FNOT,BLANKS/3HNOT,3H /
C
4004 FORMAT('1BMDX68 - MULTIPLE TIME SERIES ANALYSIS - REVISED ',
1'MAY 13, 1968'/
X41H HEALTH SCIENCES COMPUTING FACILITY, UCLA//)
C
C TOLERANCE FOR CHECK ON SINGULARITY
TOL1=1.0/(10.0**8)
C
MTAPE=5
CALL USAGEB('BMDX68')
5 ASSIGN 906 TO IWRIT
C READ CONTROL CARDS (GENERALLY PROBLM OR FINISH CARDS)
10 READ(5,1000)CODE,GOB,NSERS,NFREQ,INPUT,ICOH,ICOND,FMULT,MATFRQ,
XINVERS,IGAIN,IBANDS,NDEGF, CONST,FREQ,DELTA,UNIT,KVF,NTAPE,
XKREW,IALT
IF(CODE.EQ.FIN)GO TO 900
50 IF(CODE.NE.PROB)GO TO 905
C PRINT PROGRAM NAME AND VERSION DATE
60 WRITE(6,4004)
C CHECK PARAMETERS AND SET FLAGS
110 ASSIGN 906 TO IWRIT
ASSIGN 650 TO KSKIP
IF(INVERS.EQ.IYES)GO TO 115
ASSIGN 382 TO KSKIP
INVERS=INO
115 IF(KREW.EQ.IYES)GO TO 120
KREW=INO
120 INPRNT=1
IF(INPUT.EQ.IYES)GO TO 157
156 INPUT=INO
INPRNT=2
157 ASSIGN 500 TO KCOND
IF(ICOND.EQ.IYES)GO TO 162
158 ICOND=INO
ASSIGN 335 TO KCOND
162 ASSIGN 340 TO IMULT
IF(FMULT.EQ.YES)GO TO 163
1625 MULT=INO
ASSIGN 350 TO IMULT
163 ASSIGN 360 TO IMAT
IF(MATFRQ.EQ.IYES)GO TO 1635
1630 MATFRQ=INO
ASSIGN 380 TO IMAT
1635 ASSIGN 384 TO IPHASE
IF(IGAIN.EQ.IYES)GO TO 1640
164 IGAIN=INO
ASSIGN 383 TO IPHASE
1640 ASSIGN 425 TO ICON
IF(IBANDS.EQ.IYES)GO TO 1650
1645 IBANDS=INO
ASSIGN 395 TO ICON
1650 ASSIGN 265 TO KPAIR
IF(ICOH.EQ.IYES)GO TO 1654
1652 ICOH =INO
ASSIGN 275 TO KPAIR
C PRINT PROBLEM CARD
1654 WRITE(6,4005)GOB,NSERS,NFREQ,INPUT,ICOH,ICOND,MULT,MATFRQ,INVERS,
XIGAIN,IBANDS,NDEGF,CONST,FREQ,DELTA,UNIT,KVF
IF(NSERS.LT.2.OR.NSERS.GT.10)GO TO 910
IF(NFREQ.LT.1.OR.NFREQ.GT.999)GO TO 910
INSERS=NSERS
C READ INPUT AND OUTPUT SERIES CARDS
125 READ(5,1001)CODE,NIN ,(KIN1(I),I=1,9),NOUT ,(KOUT1(I),I=1,9)
130 IF(CODE.NE.RIN)GO TO 915
135 IF(NIN.EQ.0.OR.NIN.GT.9)GO TO 910
IF(NOUT.EQ.0.OR.NOUT.GT.9)GO TO 910
IF((NIN+NOUT).GT.NSERS)GO TO 910
C
C ORDER THE INPUT AND OUTPUT SERIES FROM LOW TO HIGH.
C
LARGE=-1
J=1
136 IRANK=11
DO 142 I=1,NIN
IF(KIN1(I)-IRANK)137,955,140
137 IF(LARGE-KIN1(I))138,140,140
138 IRANK=KIN1(I)
KSAVE(J)=IRANK
140 DO 142 K=1,NOUT
IF(KIN1(I)-KOUT1(K))142,960,142
142 CONTINUE
LARGE=IRANK
J=J+1
IF(J.LE.NIN)GO TO 136
DO 145 I=1,NIN
145 KIN1(I)=KSAVE(I)
LARGE=-1
J=1
146 IRANK=11
DO 150 I=1,NOUT
IF(KOUT1(I)-IRANK)147,965,150
147 IF(LARGE-KOUT1(I))148,150,150
148 IRANK=KOUT1(I)
KSAVE(J)=IRANK
150 CONTINUE
LARGE=IRANK
J=J+1
IF(J.LE.NOUT)GO TO 146
DO 155 I=1,NOUT
155 KOUT1(I)=KSAVE(I)
169 NLINES=32
IF(KREW.NE.IYES)GO TO 1695
1693 CALL TPWD(NTAPE,MTAPE)
GO TO 165
1695 IF(NTAPE.EQ.0)GO TO 1697
MTAPE=NTAPE
GO TO 165
1697 MTAPE=5
C READ VARIABLE FORMAT CARD(S)
165 CALL VFCHCK(KVF)
NLINES=NLINES-KVF
KVF=18*KVF
READ(5,1002)(FMT(I),I=1,KVF)
WRITE(6,4035)MTAPE,(FMT(I),I=1,KVF)
IF(MTAPE.EQ.5)GO TO 1655
FILL=BLANKS
IF(KREW.EQ.IYES)GO TO 1658
FILL=FNOT
1658 WRITE(6,4037)FILL
NLINES=NLINES-1
1655 ASSIGN 337 TO KALT
IF(-IALT)166,170,170
166 IF(IALT.EQ.6.OR.IALT.EQ.MTAPE)GO TO 950
ASSIGN 550 TO KALT
167 READ (5,1002)(FMOUT(I),I=1,12)
WRITE(6,4036)IALT,(FMOUT(I),I=1,12)
NLINES=NLINES-2
170 WRITE(6,4007)RIN, (KIN1(I),I=1,NIN)
WRITE(6,4007)OUT,(KOUT1(I),I=1,NOUT)
171 DO 172 I=1,NIN
KSAVE(I)=KIN1(I)
172 KSAVE2(I)=KSAVE(I)
J=NIN
DO 174 I=1,NOUT
J=J+1
KSAVE(J)=KOUT1(I)
174 KSAVE2(J)=KSAVE(J)
N1=NIN
N2=NOUT
IF(J.EQ.NSERS)GO TO 1705
K=J
DO 1745 I=1,NSERS
DO 1743 L=1,J
IF(KSAVE(L).EQ.I)GO TO 1745
1743 CONTINUE
K=K+1
KSAVE(K)=I
KSAVE2(K)=I
IF(K.EQ.NSERS)GO TO 1705
1745 CONTINUE
C KSAVE CONTAINS THE SERIES NUMBERS OF THE INPUT SERIES IN THE FIRST
C N1 LOCATIONS AND THE SERIES NUMBERS OF THE OUTPUT SERIES IN THE
C NEXT N2 LOCATIONS.
C KSAVE2 CONTAINS THE COLUMN NUMBERS OF THE INPUT SERIES IN THE
C MATRIX A(I,J) IN THE FIRST N1 LOCATIONS AND THE COLUMN NUMBERS OF
C THE OUTPUT SERIES IN THE NEXT N2 LOCATIONS.
C MAIN LOOP TO READ IN DATA AND DO REQUIRED COMPUTATIONS
1705 KN1 = N1
KN2 = N2
DO 1710 I=1,NSERS
JSAVE(I) = KSAVE(I)
1710 JSAVE2(I) = KSAVE2(I)
DO 800 IFREQ=1,NFREQ
NSERS=INSERS
N1 = KN1
N2 = KN2
DO 1730 I=1,NSERS
KSAVE(I) = JSAVE(I)
KSAVE2(I) = JSAVE2(I)
READ(MTAPE,FMT)(A(I,J),J=1,NSERS)
1730 CONTINUE
200 N3=N1+1
N4=N1+N2
DO 210 I=1,NSERS
DO 202 J=1,NSERS
202 B(I,J)=A(I,J)
C IF(REAL(B(I,I)))203,209,210
IF(BX(1,I,I) )203,209,210
203 DO 207 J=1,N4
IF(I.EQ.(KSAVE2(J)))GO TO 940
207 CONTINUE
GO TO 210
C FORCE DIAGONALS TO BE GREATER THAN ZERO TO PREVENT DIVIDE CHECK.
209 B(I,I)=(1.0E-25,0.0)
A(I,I)=B(I,I)
210 CONTINUE
214 GO TO (220,260),INPRNT
220 CALL INPRIN(PLUS,YMINUS)
260 GO TO KPAIR,(265,275)
C CALCULATE PAIR-WISE COHERENCES OF THE INPUT SPECTRAL MATRIX
265 N=NSERS
NLINES=NLINES-N-9
IF(NLINES)266,267,267
266 WRITE(6,4009)
NLINES=54-N
GO TO 268
267 WRITE(6,4011)
268 WRITE(6,4012)FREQ,UNIT,(KSAVE(J),J=1,NSERS)
WRITE(6,4024)
DO 274 I=1,NSERS
K=KSAVE(I)
DO 270 J= I,NSERS
L=KSAVE(J)
C TOP=((REAL(A(K,L)))**2)+((AIMAG(A(K,L)))**2)
TOP=((AX(1,K,L) )**2)+((AX(2,K,L) )**2)
C DENOM=(REAL(A(K,K)))*(REAL(A(L,L)))
DENOM=(AX(1,K,K) )*(AX(1,L,L) )
IF((ABS(DENOM)).GT.TOL1)GO TO 269
IF(TOP.GT.TOL1)GO TO 2685
CG(I,J)=1.0
GO TO 270
2685 CG(I,J)=99.0
NLINES=NLINES-2
WRITE(6,4041)
GO TO 270
269 CG(I,J)=TOP/DENOM
270 CG(J,I)=CG(I,J)
274 WRITE(6,4031)K,(CG(I,J),J=1,NSERS)
275 CALL PIVOT(RIN,OUT,PLUS,YMINUS)
IF(NSERS)790,790,2755
C PRINT CONDITIONAL SPECTRAL MATRIX, IF DESIRED
2755 GO TO KCOND,(335,500)
C WRITE ALTERNATE OUTPUT TAPE OF THE CONDITIONAL SPECTRAL MATRIX,
C IF DESIRED.
335 GO TO KALT,(337,550)
C PRINT MULTIPLE COHERENCES, IF DESIRED
337 GO TO IMULT,(340,350)
340 NLINES=NLINES-N2-7
IF(NLINES)341,342,342
341 WRITE(6,4009)
NLINES=56-N2
GO TO 343
342 WRITE(6,4011)
343 WRITE(6,4017)N2
DO 345 I=N3,N4
L=KSAVE(I)
J=KSAVE2(I)
C COHMUL=(REAL(A(J,J))-REAL(B(J,J)))/REAL(A(J,J))
COHMUL=(AX(1,J,J) -BX(1,J,J) )/AX(1,J,J)
345 WRITE(6,4018)L,COHMUL
C CORRECT FREQUENCY RESPONSE MATRIX TO RESTORE PROPER VALUE, I.E.
C THE CONJUGATE IN SOME INSTANCES. WE ARE OPERATING ON MATRIX B12, BUT
C TAKING THE CONJUGATE TO YIELD THE TRANSPOSE OF B21, THE DESIRED MATRIX
350 DO 355 I=1,N1
II=KSAVE2(I)
DO 352 J=N3,N4
JJ=KSAVE2(J)
IF(II.GT.JJ)GO TO 352
B(II,JJ)=CONJG(B(II,JJ))
352 CONTINUE
355 CONTINUE
GO TO IMAT,(360,380)
360 CALL FRQPRN(PLUS,YMINUS,E,NDEGF)
C PRINT INVERSE OF A11, IF REQUESTED
380 GO TO KSKIP,(382,650)
382 GO TO IPHASE,(384,383)
383 IF(IBANDS.EQ.INO)GO TO 790
C CALCULATE GAIN AND PHASE
384 IF(MATFRQ .EQ.IYES)GO TO 385
GO TO 605
385 DO 390 K=1,N2
IJ=K+N1
DO 390 L=1,N1
J=KSAVE2(IJ)
I=KSAVE2(L)
IF(J.LE.I)GO TO 3855
M=J
J=I
I=M
C3855 G(K,L)=SQRT((REAL(B(J,I))**2)+(AIMAG(B(J,I))**2))
3855 G(K,L)=SQRT((BX(1,J,I) **2)+(BX(2,J,I) **2))
CG(K,L)=CONST*E(K,L)
C IF(REAL(B(J,I)).NE.0.0) GO TO 389
IF(BX(1,J,I) .NE.0.0) GO TO 389
C IF(AIMAG(B(J,I)))386,387,388
IF(BX(2,J,I) )386,387,388
386 P(K,L)=0.75
GO TO 3891
387 P(K,L)=0.0
GO TO 3891
388 P(K,L)=0.25
GO TO 3891
C 389 P(K,L)=(ATAN2(AIMAG(B(J,I)),REAL(B(J,I))))*(.159154943)
389 P(K,L)=(ATAN2(BX(2,J,I) ,BX(1,J,I) ))*(.159154943)
3891 IF(G(K,L).GT.0.0)GO TO 3892
3890 CP(K,L)=0.5
GO TO 390
3892 IF(G(K,L).LE.CG(K,L))GO TO 3890
CP(K,L)=(ATAN2(CG(K,L),SQRT((G(K,L)**2)-(CG(K,L)**2))))*.159154943
390 CONTINUE
391 GO TO ICON,(395,425)
395 CALL NOCONF(G,P)
GO TO 790
425 CALL GANPHZ(G,CG,P,CP)
GO TO 790
500 CALL PRINT(PLUS,YMINUS)
GO TO 335
C WRITE ALTERNATE OUTPUT TAPE
550 IF(ICOND.EQ.IYES)GO TO 565
N5=N4-1
DO 560 I=N3,N5
II=KSAVE2(I)
N6=I+1
DO 555 J=N6,N4
JJ=KSAVE2(J)
555 B(JJ,II)=CONJG(B(II,JJ))
560 CONTINUE
565 DO 580 II=N3,N4
I=KSAVE2(II)
L=0
DO 570 JJ=N3,N4
L=L+1
J=KSAVE2(JJ)
570 U(L)=B(I,J)
580 WRITE(IALT,FMOUT)(U(J),J=1,L)
GO TO 337
C
C PRINT MATRIX OF STANDARD ERRORS OF THE FREQUENCY RESPONSE MATRIX
C
605 BOT=NDEGF-N1-N1
DO 620 II=1,N2
K=II+N1
I=KSAVE2(K)
DO 610 JJ=1,N1
J=KSAVE2(JJ)
C THE NEGATIVE (REAL(B(J,J))) USED IN THE CALCULATION BELOW IS DUE TO
C THE PIVOTING ALGORITHM WHICH YIELDS A NEGATIVE INVERSE FOR THE INPUT
C SERIES MATRIX.
C 610 E(II,JJ)=SQRT(( 2.0*(REAL(B(I,I)))*(-REAL(B(J,J))))/BOT)
610 E(II,JJ)=SQRT(( 2.0*(BX(1,I,I) )*(-BX(1,J,J) ))/BOT)
620 CONTINUE
GO TO 385
650 CALL INVPRN(PLUS,YMINUS)
GO TO 382
790 FREQ=FREQ+DELTA
800 CONTINUE
GO TO 10
C ERROR RECOVERY
900 WRITE(6,4000)
9005 IF(IALT.LT.1)GO TO901
REWIND IALT
901 IF(MTAPE-5)902,904,902
902 REWIND MTAPE
904 STOP
905 GO TO IWRIT,(906,10)
906 WRITE(6,4001) PROB,CODE
907 WRITE(6,4002)
908 ASSIGN 10 TO IWRIT
GO TO 10
910 WRITE (6,4003)CODE
GO TO 907
915 WRITE(6,4001) RIN,CODE
GO TO 907
940 K=KSAVE(J)
NLINES=NLINES-5
IF(NLINES)941,942,942
941 WRITE(6,4009)
NLINES=58
GO TO 943
942 WRITE(6,4011)
943 WRITE(6,4025)K
N1=N1-1
AKRD=RIN
IF(J.GE.N3)GO TO 944
AKRD=OUT
N1=N1+1
N2=N2-1
944 WRITE (6,4015)AKRD
N=N4-1
DO 945 L=J,N
KSAVE2(L)=KSAVE2(L+1)
945 KSAVE (L)=KSAVE (L+1)
IF(N1.EQ.0)GO TO 946
IF(N2.GT.0)GO TO 200
946 WRITE(6,4016)AKRD
GO TO 790
950 WRITE(6,4032)IALT
NLINES=NLINES-3
GO TO 167
955 FILL=RIN
957 WRITE(6,4038)FILL
IF(IALT.NE.5.AND.KREW.NE.IYES)GO TO 959
958 WRITE(6,4002)
GO TO 908
959 WRITE(6,4039)
GO TO 9005
960 WRITE(6,4040)RIN,OUT
GO TO 958
965 FILL=OUT
GO TO 957
C
C INPUT FORMATS
1000 FORMAT(2A6,I2,I3,8A3,I3,F7.0,2F6.0,A6,2I2,A3,I2)
1001 FORMAT(A6,10I2,4X,10I2)
1002 FORMAT(18A4)
C
C OUTPUT FORMATS
4000 FORMAT(45H0FINISH CARD ENCOUNTERED, PROGRAM TERMINATED.)
4001 FORMAT(23H0CONTROL CARD ERROR. A ,A6,25H CARD WAS EXPECTED BUT A ,
XA6,15H CARD WAS READ.)
4002 FORMAT(63H PROGRAM WILL GO TO THE NEXT PROBLM CARD (IF ANY) OR TER
XMINATE.)
4003 FORMAT(1H0A6,34H CARD PARAMETER OUTSIDE OF LIMITS.)
4005 FORMAT(14H PROBLEM NAME.13(1X,1H.),A6/40H TOTAL NUMBER OF SERIES R
XEAD BY PROGRAM.4X,I2/ 28H TOTAL NUMBER
X OF FREQUENCIES 6(1X,1H.),3X,I3 /30H INITIAL MATRIX OF SPECTRA AND
X/32H CROSS SPECTRA TO BE PRINTED4(1X,1H.),3X,A3/ 28H
X PRINT PAIR-WISE COHERENCES.6(1X,1H.),3X,A3/34H PRINT CONDITIONAL
XSPECTRAL MATRIX,3(1X,1H.),3X,A3/28H MULTIPLE COHERENCES DESIRED 6(
X1X,1H.),3X,A3/ 40H FREQUENCY RESPONSE MATRIX DESIRED . . .3X,A3/
X40H INVERSE OF INPUT MATRIX A11 DESIRED . .3X,A3/
X24H GAIN AND PHASE DESIRED.8(1X,1H.),3X,A3/ 26H CONFIDENCE BANDS
XDESIRED.,7(1X,1H.),3X,A3/20H DEGREES OF FREEDOM.,10(1X,1H.),3X,I3/
X26H NUMBER OF STANDARD ERRORS/ 26H FOR CONF
XIDENCE BANDS.4(1X,1H.),F12.5/
X 28H INITIAL FREQUENCY OF SERIES,3(1X,1H.),F12.5/ 34H FREQUENCY D
XIFFERENCE (DELTAFREQ).,F12.5/ 14H UNIT OF TIME.13(1X,1H.),A6/
X 40H NUMBER OF VARIABLE FORMAT CARDS . . . .3X,I3)
4007 FORMAT(1H0,33X,A6,11H ARE SERIES I3,8(1H,,I3))
4009 FORMAT(1H1)
4011 FORMAT(1H0//)
4012 FORMAT(1H ,48X,20HPAIR-WISE COHERENCES/36X,18HAT A FREQUENCY OF F1
X2.5,9H CYCLES/ A6,1H.//51X6HSERIES/11X,I2,9(10X,I2))
4015 FORMAT(22X,28HPROGRAM WILL ELIMINATE THIS A6,22H SERIES AND TRY AG
XAIN.)
4016 FORMAT(1H0,24X,3HNO A6,61H SERIES LEFT. PROGRAM WILL GO TO THE NEX
XT FREQUENCY (IF ANY).)
4017 FORMAT(1H ,37X,27HMULTIPLE COHERENCES FOR THEI3,14H OUTPUT SERIES/
X23X,8HMULTIPLE/22X,9HCOHERENCE//)
4018 FORMAT(1H 9X,6HSERIESI3,3X,F8.5)
4024 FORMAT(7H SERIES)
4025 FORMAT(1H ,19X,27HDIAGONAL ELEMENT FOR SERIESI3,50H IS NEGATIVE IN
XDICATING LEAKAGE AT THIS FREQUENCY.)
4028 FORMAT(3H I2,2X,5(F10.4,2H (F10.4,2H) ))
4030 FORMAT(53X,12HINPUT SERIES/7H OUTPUT,10X,I2,4(24X,I2))
4031 FORMAT(1H ,2X,I2,2X,10(F11.5,1X))
4032 FORMAT(1H0,I2,115H AS THE ALTERNATE TAPE NUMBER FOR THE CONDITIONA
XL MATRIX IS INCORRECT. PROGRAM CONTINUES WITHOUT WRITING THIS TAPE
X./)
4035 FORMAT(1H0,28X,32HMATRICES READ FROM LOGICAL TAPE I2,28H UNDER THE
X FOLLOWING FORMAT,/(10X,18A4))
4036 FORMAT(1H0,11X,67HCONDITIONAL SPECTRAL MATRICES TO BE WRITTEN ON L
XOGICAL TAPE NUMBER I2,28H UNDER THE FOLLOWING FORMAT,/(10X,18A4))
4037 FORMAT(1H ,28X,28HTHE ALTERNATE INPUT TAPE IS A3,15H TO BE REWOUND
X.)
4038 FORMAT(1H0,39X,4HTWO ,A6,29H SERIES HAVE THE SAME NUMBER.)
4039 FORMAT(1H0,15X,88H**WARNING** SINCE THE ALTERNATE INPUT TAPE WAS
XNOT TO BE REWOUND, THE NEXT PROBLEM CARD/15X,68HMAY CAUSE INAPPROP
XRIATE DATA TO BE READ. PROGRAM WILL BE TERMINATED.)
4040 FORMAT(1H0,25X,11HONE OF THE A6,16H AND ONE OF THE ,A6, 29H SERIES
X HAVE THE SAME NUMBER.)
4041 FORMAT(1H0,3X,112HTHE 99.0 FOR THE COHERENCE IN THE FOLLOWING LINE
X INDICATES A NON-COMPUTABLE COHERENCE DUE TO A ZERO DENOMINATOR.)
END
SUBROUTINE FRQPRN(PLUS,YMINUS,E,NDEGF)
C SUBROUTINE FRQPRN FOR BMDX68 MAY 12, 1965
C
DIMENSION A(10,10),B(10,10),T(10),U(10),ISAVE(9),KSAVE(10),KSAVE2(
X10),V(5),W(5),S(5),E(9,9)
COMMON UNIT
COMMON A,B,T,U,ISAVE,KSAVE,KSAVE2,N1,N2,N3,N4,NLINES,V,W,S,FREQ
C.....PULLOUT EQUIVALENCE, DIMENSION, AND STATEMENT CARDS MARKED 'JERK',
C.....WHEN THE CARDS USING 'AIMAG' AND 'REAL' ARE RESTORED.
DIMENSION AX(2,10,10),BX(2,10,10)
EQUIVALENCE (A,AX),(B,BX)
DOUBLE PRECISION UNIT
COMPLEX A,B,T,U
C
C PRINT FREQUENCY RESPONSE MATRIX
360 L=1
IF(N1-5)361,361,362
361 N=N1
GO TO 363
362 N=5
363 NLINES=NLINES-N2-9
IF(NLINES)3635,364,364
3635 WRITE(6,4009)
NLINES=54-N2
GO TO 365
364 WRITE(6,4011)
365 WRITE(6,4020)
WRITE(6,4019)FREQ,UNIT,(KSAVE(J),J=L,N)
WRITE(6,4024)
DO 369 IJ=N3,N4
K=0
DO 368 JI=L,N
I=KSAVE2(IJ)
J=KSAVE2(JI)
K=K+1
IF(I.LE.J)GO TO 3655
M=I
I=J
J=M
C3655 V(K)=REAL(B(I,J))
3655 V(K)=BX(1,I,J)
C W(K)=AIMAG(B(I,J))
W(K)=BX(2,I,J)
S(K)=PLUS
IF(W(K))366,367,368
366 W(K)=-W(K)
S(K)=YMINUS
GO TO 368
367 W(K)=0.0
368 CONTINUE
369 WRITE(6,4010)KSAVE(IJ),(V(J),S(J),W(J),J=1,K)
IF(N-N1)370,600,600
370 L=6
N=N1
GO TO 363
600 NLINES=NLINES-N2-9
IF(NLINES)601,602,602
601 WRITE(6,4009)
NLINES=54-N2
GO TO 603
602 WRITE(6,4011)
603 WRITE(6,4033)FREQ,UNIT
WRITE (6,4026)(KSAVE(J),J=1,N1)
WRITE (6,4024)
605 BOT=NDEGF-N1-N1
DO 620 II=1,N2
K=II+N1
I=KSAVE2(K)
DO 610 JJ=1,N1
J=KSAVE2(JJ)
C THE NEGATIVE (REAL(B(J,J))) USED IN THE CALCULATION BELOW IS DUE TO
C THE PIVOTING ALGORITHM WHICH YIELDS A NEGATIVE INVERSE FOR THE INPUT
C SERIES MATRIX.
C 610 E(II,JJ)=SQRT(( 2.0*(REAL(B(I,I)))*(-REAL(B(J,J))))/BOT)
610 E(II,JJ)=SQRT(( 2.0*(BX(1,I,I) )*(-BX(1,J,J) ))/BOT)
615 WRITE(6,4022)KSAVE(K),(E(II,JJ),JJ=1,N1)
620 CONTINUE
380 RETURN
4009 FORMAT(1H1)
4010 FORMAT(1H ,2X,I2,2X,5(F11.5,1X,A1,F11.5,1HI))
4011 FORMAT(1H0//)
4019 FORMAT(36X,18HAT A FREQUENCY OF F12.5,9H CYCLES/ A6,1H.//
X 53X,12HINPUT SERIES/7H OUTPUT11X,I2,4(23X,I2))
4020 FORMAT(1H ,46X,25HFREQUENCY RESPONSE MATRIX)
4022 FORMAT(3H ,I2,2X,10F12.5)
4024 FORMAT(7H SERIES)
4026 FORMAT(53X,12HINPUT SERIES/7H OUTPUT,5X,I2,9(10X,I2))
4033 FORMAT(1H 29X,60HSTANDARD ERROR MATRIX OF THE FREQUENCY RESPONSE C
XOEFFICIENTS/36X,18HAT A FREQUENCY OF F12.5,9H CYCLES/ A6,1H.//)
END
SUBROUTINE GANPHZ(G,CG,P,CP)
C SUBROUTINE GANPHZ FOR BMDX68 MAY 12, 1965
C
DIMENSION A(10,10),B(10,10),T(10),U(10),ISAVE(9),KSAVE(10),KSAVE2(
X10),V(5),W(5),S(5),G(9,9),CG(10,10),P(9,9),CP(10,10)
COMMON UNIT
COMMON A,B,T,U,ISAVE,KSAVE,KSAVE2,N1,N2,N3,N4,NLINES,V,W,S,FREQ
DOUBLE PRECISION UNIT
COMPLEX A,B,T,U
C
C PRINT GAIN AND PHASE WITH CONFIDENCE BANDS
425 L2=5
L1=1
IF(N1-5)427,428,428
427 L2=N1
428 NLINES=NLINES-N2-11
IF(NLINES)430,435,435
430 WRITE(6,4009)
NLINES=52-N2
GO TO 440
435 WRITE(6,4011)
440 WRITE(6,4020)
WRITE(6,4029)
WRITE(6,4021)FREQ,UNIT
WRITE(6,4030)(KSAVE(J),J=L1,L2)
WRITE(6,4024)
K=N1
DO 450 I=1,N2
K=K+1
450 WRITE(6,4028)KSAVE(K),(G(I,J),CG(I,J),J=L1,L2)
IF(L2-N1)455,460,460
455 L2=N1
L1=6
GO TO 428
460 L2=5
L1=1
IF(N1-5)462,464,464
462 L2=N1
464 NLINES=NLINES-N2-12
IF(NLINES)465,470,470
465 WRITE(6,4009)
NLINES=51-N2
GO TO 475
470 WRITE(6,4011)
475 WRITE(6,4020)
WRITE(6,4029)
WRITE(6,4023)FREQ,UNIT
WRITE(6,4030)(KSAVE(J),J=L1,L2)
WRITE(6,4024)
K=N1
DO 480 I=1,N2
K=K+1
480 WRITE(6,4028)KSAVE(K),(P(I,J),CP(I,J),J=L1,L2)
IF(L2-N1)485,790,790
485 L2=N1
L1=6
GO TO 464
790 RETURN
4009 FORMAT(1H1)
4011 FORMAT(1H0//)
4020 FORMAT(1H ,46X,25HFREQUENCY RESPONSE MATRIX)
4021 FORMAT(33X,23HGAIN AT A FREQUENCY OF F12.5,9H CYCLES/ A6,1H.//)
4023 FORMAT(33X,24HPHASE AT A FREQUENCY OF F12.5,9H CYCLES/ A6,1H./
X37X,45H* * PHASE GIVEN IN FRACTION OF A CIRCLE * *//)
4024 FORMAT(7H SERIES)
4028 FORMAT(3H I2,2X,5(F10.4,2H (F10.4,2H) ))
4029 FORMAT(1H 41X37HWITH CONFIDENCE LIMITS IN PARENTHESES)
4030 FORMAT(53X,12HINPUT SERIES/7H OUTPUT,10X,I2,4(24X,I2))
END
SUBROUTINE INPRIN(PLUS,YMINUS)
C SUBROUTINE INPRIN FOR BMDX68 MAY 12, 1965
C
DIMENSION A(10,10),B(10,10),T(10),U(10),ISAVE(9),KSAVE(10),KSAVE2(
X10),V(5),W(5),S(5)
COMMON UNIT
COMMON A,B,T,U,ISAVE,KSAVE,KSAVE2,N1,N2,N3,N4,NLINES,V,W,S,FREQ,
1NSERS
C.....PULLOUT EQUIVALENCE, DIMENSION, AND STATEMENT CARDS MARKED 'JERK',
C.....WHEN THE CARDS USING 'AIMAG' AND 'REAL' ARE RESTORED.
DIMENSION AX(2,10,10),BX(2,10,10)
EQUIVALENCE (A,AX),(B,BX)
DOUBLE PRECISION UNIT
COMPLEX A,B,T,U
C
220 IF(NSERS-5)222,222,250
222 N=NSERS
225 L=1
2250 NLINES=NLINES-NSERS-9
IF(NLINES)2205,221,221
2205 WRITE(6,4009)
NLINES=54-NSERS
GO TO 226
221 WRITE(6,4011)
C PRINT OUT INPUT MATRIX, IF DESIRED
226 WRITE(6,4008)FREQ,UNIT,(KSAVE(J),J=L,N)
WRITE(6,4024)
DO 230 II=1,NSERS
I=KSAVE2(II)
K=0
DO 229 JJ=L,N
J=KSAVE2(JJ)
K=K+1
C V(K)=REAL(A(I,J))
V(K)=AX(1,I,J)
C W(K)=AIMAG(A(I,J))
W(K)=AX(2,I,J)
S(K)=PLUS
IF(W(K))227,228,229
227 W(K)=-W(K)
S(K)=YMINUS
GO TO 229
228 W(K)=0.0
229 CONTINUE
230 WRITE(6,4010)KSAVE(II),(V(J),S(J),W(J),J=1,K)
IF(N-NSERS)245,260,260
245 L=6
N=NSERS
GO TO 2250
250 N=5
GO TO 225
260 RETURN
4008 FORMAT(1H ,48X,21HINPUT SPECTRAL MATRIX/36X,18HAT A FREQUENCY OF
XF12.5,9H CYCLES/ A6,1H.//51X,6HSERIES/20X,4(I2,22X),I2)
4009 FORMAT(1H1)
4010 FORMAT(1H ,2X,I2,2X,5(F11.5,1X,A1,F11.5,1HI))
4011 FORMAT(1H0//)
4024 FORMAT(7H SERIES)
END
SUBROUTINE INVPRN(PLUS,YMINUS)
C SUBROUTINE INVPRN FOR BMDX68 MAY 12, 1965
C
DIMENSION A(10,10),B(10,10),T(10),U(10),ISAVE(9),KSAVE(10),KSAVE2(
X10),V(5),W(5),S(5)
COMMON UNIT
COMMON A,B,T,U,ISAVE,KSAVE,KSAVE2,N1,N2,N3,N4,NLINES,V,W,S,FREQ
C.....PULLOUT EQUIVALENCE, DIMENSION, AND STATEMENT CARDS MARKED 'JERK',
C.....WHEN THE CARDS USING 'AIMAG' AND 'REAL' ARE RESTORED.
DIMENSION AX(2,10,10),BX(2,10,10)
EQUIVALENCE (A,AX),(B,BX)
DOUBLE PRECISION UNIT
COMPLEX A,B,T,U
C
C PRINT B11, THE INVERSE OF A11.
C THIS MATRIX AS CALCULATED IS THE NEGATIVE OF THE INVERSE, SO THE SIGNS
C ARE REVERSED IN THIS PRINTING.
C
650 L=1
IF(N1-5)651,651,652
651 N=N1
GO TO 653
652 N=5
653 NLINES=NLINES-N1-9
IF(NLINES)654,655,655
654 WRITE(6,4009)
NLINES=54-N1
GO TO 656
655 WRITE(6,4011)
656 WRITE(6,4034)FREQ,UNIT,(KSAVE(J),J=L,N)
WRITE (6,4024)
DO 665 II=1,N1
I=KSAVE2(II)
K=0
DO 663 JJ=L,N
J=KSAVE2(JJ)
K=K+1
IF(J.LT.I)GO TO 659
C V(K)=-REAL(B(I,J))
V(K)=-BX(1,I,J)
C W(K)=-AIMAG(B(I,J))
W(K)=-BX(2,I,J)
657 S(K)=PLUS
IF(W(K))658,660,663
658 W(K)=-W(K)
S(K)=YMINUS
GO TO 663
C 659 V(K)=-REAL(B(J,I))
659 V(K)=-BX(1,J,I)
C W(K)=+AIMAG(B(J,I))
W(K)= BX(2,J,I)
GO TO 657
660 W(K)=0.0
663 CONTINUE
665 WRITE (6,4010)KSAVE(II),(V(J),S(J),W(J),J=1,K)
IF(N-N1)670,382,382
670 L=6
N=N1
GO TO 653
382 RETURN
4009 FORMAT(1H1)
4010 FORMAT(1H ,2X,I2,2X,5(F11.5,1X,A1,F11.5,1HI))
4011 FORMAT(1H0//)
4024 FORMAT(7H SERIES)
4034 FORMAT(1H ,39X,39HINVERSE OF INPUT SERIES SPECTRAL MATRIX/36X,18H
XAT A FREQUENCY OF F12.5,9H CYCLES/ A6,1H.//51X,6HSERIES/20X,4(I2,
X22X),I2)
END
SUBROUTINE NOCONF(G,P)
C SUBROUTINE NOCONF FOR BMDX68 MAY 12, 1965
C
DIMENSION A(10,10),B(10,10),T(10),U(10),ISAVE(9),KSAVE(10),KSAVE2(
X10),V(5),W(5),S(5),G(9,9),P(9,9)
COMMON UNIT
COMMON A,B,T,U,ISAVE,KSAVE,KSAVE2,N1,N2,N3,N4,NLINES,V,W,S,FREQ
DOUBLE PRECISION UNIT
COMPLEX A,B,T,U
C
C PRINT GAIN AND PHASE WITHOUT CONFIDENCE BANDS.
395 NLINES=NLINES-N2-9
IF(NLINES)396,397,397
396 WRITE(6,4009)
NLINES=54-N2
GO TO 398
397 WRITE(6,4011)
398 WRITE(6,4020)
WRITE(6,4021)FREQ,UNIT
WRITE(6,4026)(KSAVE(J),J=1,N1)
WRITE(6,4024)
DO 400 I=1,N2
K=I+N1
400 WRITE(6,4022)KSAVE(K),(G(I,J),J=1,N1)
NLINES=NLINES-N2-10
IF(NLINES)401,402,402
401 WRITE(6,4009)
NLINES=53-N2
GO TO 403
402 WRITE(6,4011)
403 WRITE(6,4020)
WRITE(6,4023)FREQ,UNIT
WRITE(6,4026)(KSAVE(J),J=1,N1)
WRITE(6,4024)
DO 405 I=1,N2
K=I+N1
405 WRITE(6,4022)KSAVE(K),(P(I,J),J=1,N1)
RETURN
4009 FORMAT(1H1)
4011 FORMAT(1H0//)
4020 FORMAT(1H ,46X,25HFREQUENCY RESPONSE MATRIX)
4021 FORMAT(33X,23HGAIN AT A FREQUENCY OF F12.5,9H CYCLES/ A6,1H.//)
4022 FORMAT(3H ,I2,2X,10F12.5)
4023 FORMAT(33X,24HPHASE AT A FREQUENCY OF F12.5,9H CYCLES/ A6,1H./
X37X,45H* * PHASE GIVEN IN FRACTION OF A CIRCLE * *//)
4024 FORMAT(7H SERIES)
4026 FORMAT(53X,12HINPUT SERIES/7H OUTPUT,5X,I2,9(10X,I2))
END
SUBROUTINE PIVOT(RIN,OUT,PLUS,YMINUS)
C SUBROUTINE PIVOT FOR BMDX68 MAY 11, 1965
C
COMMON UNIT
COMMON A,B,T,U,ISAVE,KSAVE,KSAVE2,N1,N2,N3,N4,NLINES,V,W,S,FREQ,
1NSERS,FMULT
DIMENSION A(10,10),B(10,10),T(10),U(10),ISAVE(9),KSAVE(10),KSAVE2(
X10),V(5),W(5),S(5),JSAVE(9),JSAVE2(9)
C.....PULLOUT EQUIVALENCE, DIMENSION, AND STATEMENT CARDS MARKED 'JERK',
C.....WHEN THE CARDS USING 'AIMAG' AND 'REAL' ARE RESTORED.
DIMENSION AX(2,10,10),BX(2,10,10),XAB(2)
EQUIVALENCE (A,AX),(B,BX),(XA3,XAB)
DOUBLE PRECISION UNIT,RIN,OUT,AKRD
COMPLEX A,B,T,U,XA3
DATA YES/3HYES/
C
C TOLERANCE CHECK FOR SINGULARITY
C
TOL = 1.0 / (10.0**6)
C FIND MAXIMUM DIAGONAL ELEMENT OF A11 AND
C PIVOT TO OBTAIN REQUIRED MATRICES
C CODING OBTAINED FROM SUBROUTINE STEP OF BMD02R DATED MARCH 19,1964
C THIS METHOD YIELDS A NEGATIVE INVERSE IN THE A11 LOCATION OF THE
C PARTITIONED MATRIX AND CORRECT VALUES ELSEWHERE.
10 N3=N1+1
N4=N1+N2
DO 325 K=1,N1
ISAVE(K)=0
BIG=-10.0
DO 2765 IK=1,N1
II=KSAVE2(IK)
C IF(BIG.GE.(REAL(B(II,II))/REAL(A(II,II))))GO TO 2765
IF(BIG.GE.(BX(1,II,II) /AX(1,II,II) ))GO TO 2765
DO 2763 J=1,K
IF(ISAVE(J).EQ.II)GO TO 2765
2763 CONTINUE
C BIG=REAL(B(II,II))/REAL(A(II,II))
BIG=BX(1,II,II) /AX(1,II,II)
ISAVE(K)=II
2765 CONTINUE
KAY=ISAVE(K)
XA3=B(KAY,KAY)
KAY1=KAY-1
KAY2=KAY+1
C CHECK FOR SINGULARITY
IF(BIG .LT. TOL) GO TO 920
IF(-KAY1)280,290,290
280 DO 285 I=1,KAY1
U(I)=B(I,KAY)
T(I)=CONJG(U(I))
B(I,KAY)=(0.0,0.0)
285 CONTINUE
290 U(KAY)=(-1.0,0.0)
T(KAY)=U(KAY)
B(KAY,KAY)=(0.0,0.0)
IF(KAY2-NSERS)295,295,305
295 DO 300 I=KAY2,NSERS
T(I)=B(KAY,I)
U(I)=CONJG(T(I))
B(KAY,I)=(0.0,0.0)
300 CONTINUE
305 DO 310 I=1,NSERS
B(I,I) = B(I,I) - (U(I)*T(I))/XA3
DO 306 KK=N3,N4
IF(KSAVE2(KK) .EQ. I) GO TO 307
306 CONTINUE
GO TO 308
C307 IF(REAL(B(I,I)) .LT. 0.0) GO TO 312
307 IF(BX(1,I,I) .LT. 0.0) GO TO 312
308 KK = I + 1
DO 310 J=KK,NSERS
B(I,J) = B(I,J) - (U(I)*T(J)) / XA3
310 CONTINUE
GO TO 325
312 DO 315 I=1,N1
KK = KSAVE2(I)
IF(KK .EQ. ISAVE(K)) GO TO 316
315 CONTINUE
316 WRITE(6,4027) KSAVE(I)
WRITE(6,4015) RIN
N1 = N1 - 1
N4 = N4 - 1
DO 318 II=I,N4
KSAVE(II) = KSAVE(II+1)
318 KSAVE2(II) = KSAVE2(II+1)
NLINES = 0
DO 320 I=1,NSERS
DO 320 J=1,NSERS
320 B(I,J) = A(I,J)
IF(N1 .GT. 0) GO TO 10
GO TO 939
325 CONTINUE
326 IF(FMULT .NE. YES) GO TO 950
NLINES = NLINES - 5
IF(NLINES) 5,101,101
5 WRITE(6,4009)
NLINES = 58
GO TO 15
101 WRITE(6,4011)
C15 TOL = -1.0 / (REAL(A(KAY,KAY))*REAL(B(KAY,KAY)))
15 TOL = -1.0 / (AX(1,KAY,KAY) * BX(1,KAY,KAY))
WRITE(6,4008) TOL
GO TO 950
920 N = K - 1
L = 0
DO 922 I=1,N1
KK = KSAVE2(I)
DO 921 J=1,N
IF(ISAVE(J) .EQ. KK) GO TO 922
921 CONTINUE
L = L + 1
JSAVE2(L) = KK
JSAVE(L) = KSAVE(I)
922 CONTINUE
WRITE(6,4013)
WRITE(6,4014) (JSAVE(I), I=1,L)
K = 0
DO 925 I=1,N1
DO 924 J=1,L
IF(JSAVE(J) .EQ. KSAVE(I)) GO TO 925
924 CONTINUE
K = K + 1
KSAVE(K) = KSAVE(I)
KSAVE2(K) = KSAVE2(I)
925 CONTINUE
930 WRITE(6,4015) RIN
N3 = K + 1
N4 = K + N2
DO 935 I=N3,N4
J = N1 - K + I
KSAVE(I) = KSAVE(J)
935 KSAVE2(I) = KSAVE2(J)
N1 = K
NLINES = 55 - N
KAY = ISAVE(N)
IF(N1 .GT. 0) GO TO 326
939 WRITE(6,4016) RIN
940 NSERS = -NSERS
950 RETURN
4008 FORMAT(1H ,37X,33HTOLERANCE FOR THE INPUT SERIES = F8.5)
4009 FORMAT(1H1)
4011 FORMAT(1H0//)
4013 FORMAT(1H1,48X,21H** SINGULAR MATRIX **/2X,6HSERIESI3,37H IS A LIN
XEAR COMBINATION OF THE FIRSTI3,26H SERIES AT A FREQUENCY OF F12.5,
X8H CYCLES/A6,12H AS FOLLOWS.//)
4014 FORMAT(1H ,35X,6HSERIESI3,13H CONSTANT IS ,F12.5,A1,F12.5,1HI)
4015 FORMAT(22X,28HPROGRAM WILL ELIMINATE THIS A6,22H SERIES AND TRY AG
XAIN.)
4016 FORMAT(1H0,24X,3HNO A6,61H SERIES LEFT. PROGRAM WILL GO TO THE NEX
XT FREQUENCY (IF ANY).)
4027 FORMAT(22X,6HSERIESI3,51H CAUSES THE MATRIX TO BE NON POSITIVE DEF
XINITE. THE)
END
SUBROUTINE PRINT(PLUS,YMINUS)
C SUBROUTINE PRINT FOR BMDX68 MAY 12, 1965
C
DIMENSION A(10,10),B(10,10),T(10),U(10),ISAVE(9),KSAVE(10),KSAVE2(
X10),V(5),W(5),S(5)
COMMON UNIT
COMMON A,B,T,U,ISAVE,KSAVE,KSAVE2,N1,N2,N3,N4,NLINES,V,W,S,FREQ
C.....PULLOUT EQUIVALENCE, DIMENSION, AND STATEMENT CARDS MARKED 'JERK',
C.....PULLOUT EQUIVALENCE, DIMENSION, AND STATEMENT CARDS MARKED 'JERK',
C.....WHEN THE CARDS USING 'AIMAG' AND 'REAL' ARE RESTORED.
DIMENSION AX(2,10,10),BX(2,10,10)
EQUIVALENCE (A,AX),(B,BX)
DOUBLE PRECISION UNIT
COMPLEX A,B,T,U
C
C PRINT CONDITIONAL SPECTRAL MATRIX
C
500 IF(N2-5)501,501,540
501 N=N4
502 L=N3
504 NLINES=NLINES-N2-8
IF(NLINES)505,510,510
505 WRITE(6,4009)
NLINES=55-N2
GO TO 515
510 WRITE(6,4011)
515 WRITE(6,4006)FREQ,UNIT,(KSAVE(K),K=L,N)
WRITE(6,4024)
DO 525 I=N3,N4
II=KSAVE2(I)
K=0
DO 520 J=L,N
K=K+1
JJ=KSAVE2(J)
IF(I.GE.J)GO TO 517
B(JJ,II)=CONJG(B(II,JJ))
C 517 V(K)=REAL(B(II,JJ))
517 V(K)=BX(1,II,JJ)
C W(K)=AIMAG(B(II,JJ))
W(K)=BX(2,II,JJ)
S(K)=PLUS
IF(W(K))518,519,520
518 W(K)=-W(K)
S(K)=YMINUS
GO TO 520
519 W(K)=0.0
520 CONTINUE
525 WRITE(6,4010)KSAVE(I),(V(J),S(J),W(J),J=1,K)
IF(N-N4)530,335,335
530 L=N1+6
N=N4
GO TO 504
540 N=N1+5
GO TO 502
335 RETURN
4006 FORMAT(1H 22X,46HCONDITIONAL SPECTRAL MATRIX AT A FREQUENCY OF F12
X.5,9H CYCLES/ A6,1H.//51X,6HSERIES/20X,4(I2,22X),I2)
4009 FORMAT(1H1)
4010 FORMAT(1H ,2X,I2,2X,5(F11.5,1X,A1,F11.5,1HI))
4011 FORMAT(1H0//)
4024 FORMAT(7H SERIES)
END
SUBROUTINE TPWD(NT1,NT2)
C SUBROUTINE TPWD JULY 15, 1964
C
IF(NT1)40,10,12
10 NT1=5
12 IF(NT1-NT2)14,19,14
14 IF(NT2-5)15,19,17
15 REWIND NT2
GO TO 19
17 REWIND NT2
19 IF(NT1-5)18,24,18
18 IF(NT1-6)22,40,22
22 REWIND NT1
24 NT2=NT1
28 RETURN
40 WRITE (6,49)
49 FORMAT(25H ERROR ON TAPE ASSIGNMENT)
STOP
END
SUBROUTINE VFCHCK(NVF)
C SUBROUTINE VFCHCK JULY 15, 1964
C
IF(NVF)10,10,20
20 IF (NVF-10) 50,50,10
10 WRITE (6,4000)
NVF=1
50 RETURN
C
C
C
4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
XIED, ASSUMED TO BE 1.)
END