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

























