CBMD02T  AUTOCOVARIANCE AND POWER SPECTRAL ANALYSIS  AUGUST 13, 1966
      DIMENSION Z(2000),DUMB(20),FMT(36),NSELL(20),RX(30),PX(30),
     1 LISA(20),ADDXY(30),SUBXY(30),CXY(30),QXY(30),SCXY(30),
     2 SQXY(30),AMXY(30),PHASXY(30),RXY(30),RYX(30),YP( 4),
     3 CRO(60),COXYSQ(30),IPRNT1(20),IPRNT2(20)
      DIMENSION X(100),Y(100),SPX(30),SPY(30),IPRNT3(20)
      COMMON UNIT
      COMMON  RXY    , RYX    , RX     , SUBXY
      COMMON  Z      , DUMB   , FMT    , NSELL  , LISA   , SCXY
      COMMON  SQXY   , YP              , CRO    , COXYSQ , IORDAT
      COMMON  IPLDAT , IDETRN , IPREWT , THETA  , KSER   , NPOINT
      COMMON  LAG    , NSEL   , IPLOT  , DELTAT          , INFORM
      COMMON  KVR    , IPRNT1 , IPRNT2 , IPRNT3 , IPOW   , ICCS
      COMMON  ICOH   , PI     , FNPT   , CONST  , PMD    , FLAG
      COMMON  FLLAG  , LLAG   , KIT    , KDUM   , MISTAK , MX
      COMMON  KX     , Q      , KPOINT
      COMMON SYM
C
      EQUIVALENCE (RXY,CXY,AMXY),(RYX,QXY,PHASXY),(RX,ADDXY),(PX,SUBXY),
     1(YES,IYES),(END,IEND),(PROBLM,IPROB),(SEL,ISEL)
C
 9002 FORMAT(1H177H BMD02T-AUTOCOVARIANCE AND POWER SPECTRAL ANALYSIS -
     1REVISED JANUARY 15, 1969/
     22X40HHEALTH SCIENCES COMPUTING FACILITY, UCLA//)
C
      DOUBLE PRECISION IPROB
	DOUBLE PRECISION PROBLM,SEL,END,JOB,PROB,UNIT
      DATA PROBLM/6HPROBLM/,SEL/6HSELECT/,END/6HFINISH/
      DATA IYES/3HYES/
      INTEGER SYM(4)
C
      PI=3.141593
      INFORM=5
	CALL USAGEB('BMD02T')
C
 999  READ (5,5001)JOB,PROB,IORDAT,IPLDAT,IDETRN,IPREWT,THETA, KSER,NPOI
     1NT,LAG,NSEL,IPLOT,DELTAT,UNIT,IOLD,NTAPE,KVR
 5001 FORMAT(2A6,4A3,F5.0,I2,I4,I3,2I2,F5.0,A6,A3,12X,2I2)
      IF( END .NE. JOB) GO TO 10
      WRITE (6,887)
  887 FORMAT(1H0 / / / '0FINISH CARD ENCOUNTERED.' )
  888 IF (INFORM .EQ. 5) GO TO 890
  889 REWIND INFORM
  890 STOP
   10 IF (IPROB .NE. JOB) GO TO 740
 150  WRITE (6,9002)
      WRITE (6,9003)PROB
 9003 FORMAT(1H ///40X,26HP R O B L E M  N U M B E R,4X,A6/40X,26(1H*)//
     1/)
      WRITE (6,9005)IORDAT,IPLDAT,IDETRN,IPREWT,THETA,KSER, NPOINT,LAG,N
     1SEL,IOLD
 9005 FORMAT(20X,31HINPUT DATA TO BE PRINTED OUT---A3//20X,33HINPUT SERI
     1ES TO BE PLOTTED OUT---A3//20X,13HDETRENDING---A3//20X,15HPREWHITE
     2NING---A3//20X,81HVALUE OF CONSTANT C USED IN THE PREWHITENING TRA
     3NSFORMATION Z(T)=X(T+1)-CX(T) ---F10.5//20X19HNUMBER OF SERIES = I
     42//20X,23HNUMBER OF DATA POINTS =I5//20X,24HNUMBER OF LAGS CHOSEN 
     5= I3//20X,28HNUMBER OF SELECTION CARDS = I2//20X,20HUSE PREVIOUS D
     6ATA---A3//)
      IF(KSER*(KSER-21))904,750,750
  904 WRITE (6,9013)DELTAT,UNIT
 9013 FORMAT(20X,25HCONSTANT TIME INTERVAL = F10.5,1X,A6//)
      IF(NPOINT*(NPOINT-101))909,760,760
  909 WRITE (6,9021)KVR
 9021 FORMAT(20X,34HNUMBER OF VARIABLE FORMAT CARDS = I3//)
      KPOINT=NPOINT
      CTHETA=ABS(THETA)
      IF(LAG-30)8,8,730
    8 IF((NPOINT*KSER)-2000)9,9,720
 9    IF(CTHETA-1.0) 99,99,710
 99   IF(IOLD.EQ.IYES)GO TO 30
 11   ASSIGN 75 TO NSKIP
 12   IF(KVR)40,705,20
 40   IF(KVR+1)42,45,42
 42   L=1
      ASSIGN 74 TO NSKIP
      GO TO 46
 45   L=2
 46   KVR=1
      GO TO (25,20),L
 20   CALL VFCHCK(KVR)
      KVR=KVR*18
      READ (5,1002)(FMT(I),I=1,KVR)
 1002 FORMAT (18A4)
      WRITE (6,1003) (FMT(I),I= 1, KVR)
 1003 FORMAT( '0VARIABLE FORMAT IS' / (20X,18A4))
 25   ASSIGN 80 TO KSKIP
      ASSIGN 30 TO LSKIP
      IF ( NTAPE) 71,72,73
 71   KADD=-NPOINT
       DO 84 K=1,KSER
      KADD=KADD+NPOINT
      GO TO NSKIP,(74,75,715)
  715 READ (INFORM)(X(I),I=1,NPOINT)
      GO TO 755
C
 74   READ (5,1002)(FMT(J),J=1,18)
      WRITE (6,1003) (FMT(I),I= 1, KVR)
 75   READ (5,FMT)(X(I),I=1,NPOINT)
  755 DO 84 I=1,NPOINT
      ND=I+KADD
   84 Z(ND)=X(I)
      GO TO 825
C
   73 ASSIGN 83 TO LSKIP
      IF(70- NTAPE)735,735,738
 735  NTAPE=NTAPE-70
      ASSIGN 715 TO NSKIP
      CALL TPWD(NTAPE,INFORM)
      GO TO 71
C
 738  CALL TPWD(NTAPE,INFORM)
      ASSIGN 81 TO KSKIP
   72 DO 82 I=1,NPOINT
      KADD=0
      GO TO KSKIP,(80,81)
 80   READ (5,FMT)(DUMB(J),J=1,KSER)
      GO TO 815
C
 81   READ (INFORM)(DUMB(J),J=1,KSER)
 815  DO 82 J=1,KSER
      L=KADD+I
      KADD=KADD+NPOINT
   82 Z(L)=DUMB(J)
  825 GO TO LSKIP,(30,83)
 83   REWIND INFORM
   30 LLAG=LAG+1
      FLLAG=LLAG
      FLAG=LAG
      FNPT=NPOINT
      Q=PI/FLAG
      CONST=2.0*DELTAT/PI
      PMD=0.5/(FLAG*DELTAT)
      IF(-NSEL)95,770,770
   95 DO 9000 KDUM=1,NSEL
      KIT=0
      DO 98 I=1,20
      NSELL(I)=0
      IPRNT1(I)=0
      IPRNT2(I)=0
      IPRNT3(I)=0
      LISA(I)=0
 98   CONTINUE
      READ (5,5003)JOB,IPOW,ICCS,ICOH,NNX,IPRNT1(1),LISA(1),NNS,(NSELL(I
     1),IPRNT1(I+1),IPRNT2(I+1),IPRNT3(I+1),LISA(I+1),I=1,NNS)
 5003 FORMAT(A6,1X,3A4,I2,A1,2I2,1X,6(I2,3A1,I2)/6X,9(I2,3A1,I2)/6X,
     14(I2,3A1,I2))
      IF (JOB .NE. SEL) GO TO 745
 101  CALL SMEAN (NNX,X,XBAR,XALPHA)
 102  IF (MISTAK) 9000,100,9000
 100  ASSIGN 145 TO KPOWR
      IF(IPOW.NE.IYES)GO TO 105
 110  ASSIGN 135 TO KPOWR
      CALL POWER (NNX,X,XBAR,XALPHA,SPX)
  105 DO 8000 I=1,NNS
      KIT=I
      NNY=NSELL(I)
      IF(NNX-NNY)120,9000,120
  120 CALL SMEAN (NNY,Y,YBAR,YALPHA)
      IF(MISTAK-17)122,740,122
 122  IF(MISTAK) 9000,125,9000
 125  GO TO KPOWR,(135,145)
  135 CALL POWER (NNY,Y,YBAR,YALPHA,SPY)
      IF(ICCS.NE.IYES)GO TO 8000
 140  IF(IPREWT.EQ.IYES)GO TO 130
  145 CALL CROSS (NNX,NNY,X,Y,SPX,SPY,XBAR,YBAR,XALPHA,YALPHA)
      GO TO 8000
  130 WRITE (6,1051)
 1051 FORMAT(1H1,10X,77HSINCE PREWHITENING IS DONE IN THE INPUT DATA , C
     1ROSS-SPECTRUM IS NOT COMPUTED/ 10X,66HPROGRAM WILL PROCEED TO NEXT
     2 SERIES IN THE SELECTION CARD (IF ANY)/ 10X,39H OR TO THE NEXT SEL
     3ECTION CARD (IF ANY))
 8000 CONTINUE
 9000 CONTINUE
      GO TO 999
C
  705 IF(INFORM)20,20,73
C
 710  WRITE (6,9200)THETA
      GO TO 751
C
 720  WRITE (6,9300)
      GO TO 751
C
 730  WRITE (6,9400)LAG
      GO TO 751
C
  740 CONTINUE
      WRITE (6,9501)
 9501 FORMAT('0PROGRAM EXPECTS A ''PROBLM'' CARD.')
      WRITE (6,9500) JOB,PROB,IORDAT,IPLDAT,IDETRN,IPREWT,THETA,KSER,
     1 NPOINT,LAG,NSEL,IPLOT,DELTAT,UNIT,IOLD,NTAPE,KVR
      GO TO 751
C
  745 CONTINUE
      WRITE (6,9551)
      WRITE (6,9550) JOB,IPOW,ICCS,ICOH,NNX,IPRNT1(1),LISA(1),NNS
      GO TO 751
C
C
  750 WRITE (6,9600)
  751 WRITE (6,9800)
      GO TO 888
C
  760 WRITE (6,9700)
      GO TO 751
C
  770 WRITE (6,9900)
      NSEL =1
      GO TO 95
C
 9200 FORMAT(1H0,20X,47HCONSTANT USED IN PREWHITENING TRANSFORMATION IS,
     1F7.5,23H,A VALUE NOT PERMITTED.)
 9300 FORMAT(1H0,37X,22HDATA STORAGE EXCEEDED.)
 9400 FORMAT(1H0,37X,29HNUMBER OF LAGS IS TOO LARGE =,I5)
 9500 FORMAT( '0CONTROL CARD OUT OF ORDER OR MISPUNCHED. CARD IS PRINTED
     1 BELOW. ' / 1X,2A6,4A3,F5.0,I2,I4,I3,2I2,F5.0,A6,A3,12X,2I2)
 9550 FORMAT('0CONTROL CARD OUT OF ORDER OR MISPUNCHED. CARD IS PRINTED 
     1BELOW.' / 1X,A6,1X,3A4,I2,A1,2I2)
 9551 FORMAT('0PROGRAM EXPECTS A ''SELECT'' CARD.')
 9600 FORMAT(1H0,37X,43HNUMBER OF SERIES IS OUTSIDE PROGRAM LIMITS.)
 9700 FORMAT(1H0,37X,43HNUMBER OF POINTS IS OUTSIDE PROGRAM LIMITS.)
 9800 FORMAT(1H0,37X,23HPROGRAM CANNOT PROCEED.)
 9900 FORMAT(1H0,37X,44HNO SELECTION CARD SPECIFIED. ONE IS ASSUMED.)
C
      END
CCROSS        SUBROUTINE CROSS FOR BMD02T            APRIL 20, 1966
      SUBROUTINE CROSS(NNX,NNY,X,Y,SPX,SPY,XBAR,YBAR,XALPHA,YALPHA)
C
      DIMENSION Z(2000),DUMB(20),FMT(36),NSELL(20),RX(30),PX(30),
     1 LISA(20),ADDXY(30),SUBXY(30),CXY(30),QXY(30),SCXY(30),
     2 SQXY(30),AMXY(30),PHASXY(30),RXY(30),RYX(30),YP( 4),
     3 CRO(60),COXYSQ(30),IPRNT1(20),IPRNT2(20)
      DIMENSION X(100),Y(100),SPX(30),SPY(30),IPRNT3(20)
      DOUBLE PRECISION UNIT
      INTEGER SYM(4)
      COMMON UNIT
      COMMON  RXY    , RYX    , RX     , SUBXY
      COMMON  Z      , DUMB   , FMT    , NSELL  , LISA   , SCXY
      COMMON  SQXY   , YP              , CRO    , COXYSQ , IORDAT
      COMMON  IPLDAT , IDETRN , IPREWT , THETA  , KSER   , NPOINT
      COMMON  LAG    , NSEL   , IPLOT  , DELTAT          , INFORM
      COMMON  KVR    , IPRNT1 , IPRNT2 , IPRNT3 , IPOW   , ICCS
      COMMON  ICOH   , PI     , FNPT   , CONST  , PMD    , FLAG
      COMMON  FLLAG  , LLAG   , KIT    , KDUM   , MISTAK , MX
      COMMON  KX     , Q      , KPOINT
      COMMON SYM
      EMDA(Q000FL,Q001FL,Q002FL)=((Q000FL*Q001FL)**2)*(1.0-2.0*(Q002FL/Q
     1000FL)-2.0*(Q002FL/Q000FL)**2)
C
      EQUIVALENCE (RXY,CXY,AMXY),(RYX,QXY,PHASXY),(RX,ADDXY),(PX,SUBXY),
     1(YES,IYES),(B,IB),(C,IC),(SS,IS),(H,IH),(T,IT)
C
      DATA B/1HB/,BLANK/1H /,C/1HC/,DBLSTR/2H**/,H/1HH/,SS/1HS/,STAR/
     11HX/,T/1HT/,YES/3HYES/
      INTEGER X007CT
      DATA X007CT/1H*/
      SYM(1)=X007CT
C
      DO 501 I=1,LLAG
      SXY=0.0
      SYX=0.0
      M=I-1
      L=NPOINT-M
      FL=L
      DO 502 J=1,L
      K=M+J
      SXY=SXY+X(J)*Y(K)
  502 SYX=SYX+X(K)*Y(J)
      RXY(I)=SXY/FL
  501 RYX(I)=SYX/FL
 455  WRITE (6,1501)NNX,NNY,NNX,NNY,UNIT
 1501 FORMAT(1H1,21X,3HLAG,14X,15HCROSSCOVARIANCE,15X,15HCROSSCOVARIANCE
     1/37X,9HOF SERIES,I3,4H AND,I3,11X,9HOF SERIES,I3,4H AND,I3/20X,1H(
     2A6,1H),11X,15H(POSITIVE  TAU),15X,15H(NEGATIVE  TAU)//)
      TIME=-DELTAT
      DO 503 I=1,LLAG
      TIME=TIME+DELTAT
  503 WRITE (6,1502)TIME,RXY(I),RYX(I)
 1502 FORMAT(19X,F9.4,10X,F15.6,15X,F15.6)
      TIME1=TIME
  509 DO510 I=1,LLAG
      K=LLAG-I+1
  510 CRO(I)=RXY(K)
      KL=LLAG
      DO 511 I=2,LLAG
      KL=KL+1
  511 CRO(KL)=RYX(I)
      LLLAG=2*LLAG-1
      IF(IPRNT2(KIT+1).EQ.IB)GO TO 520
 515  IF(IPRNT2(KIT+1).NE.IC)GO TO 655
 520  FR=-10.0**10
      RO=10**10
      DO630 I=1,LLLAG
      FR=AMAX1(FR,CRO(I))
  630 RO=AMIN1(RO,CRO(I))
      WRITE (6,1503)NNX,NNY,TIME1,UNIT
 1503 FORMAT(1H1, 7X,44HGRAPH OF CROSS-COVARIANCE FUNCTION OF SERIES,I3,
     14H AND,I3,36H PLOTTED AGAINST TIME UP TO A LAG OF,F9.4,1X,A6//)
      ANY=-TIME1-DELTAT
      IF(LLAG-50) 640,645,645
 640  K1=0
      K2=1
      ASSIGN 642 TO KSKIP
      GO TO 641
C
 645  K1=-1
      K2=K1
      ASSIGN 643 TO KSKIP
      KZ=LLLAG+1
 641   DO 650 I=1,LLLAG
      ANY=ANY+DELTAT
      GO TO KSKIP,(642,643)
 642  YP(1)=CRO(I)
      GO TO 650
C
 643  KZ=KZ-1
      YP(1)=CRO(KZ)
 650  CALL PLOTR(ANY,-TIME1,TIME1,YP,SYM,RO,FR,1,K2)
      CALL PLOTR(ANY,-TIME1,TIME1,YP,SYM,RO,FR,K1,K2)
  655 DO 15 JP=1,LLAG
      ADDXY(JP)=RXY(JP)+RYX(JP)
   15 SUBXY(JP)=RXY(JP)-RYX(JP)
      ADDXY(1)=0.5*ADDXY(1)
      ADDXY(LLAG)=0.5*ADDXY(LLAG)
      SUBXY(1)=0.5*SUBXY(1)
      SUBXY(LLAG)=0.5*SUBXY(LLAG)
      S1=-Q
      DO 17 IH=1,LLAG
      CXY (IH)=0.0
      QXY(IH)=0.0
      S1=S1+Q
      S=-S1
      DO 17 JP=1,LLAG
      S=S+S1
      CXY(IH)=CXY(IH)+COS(S)*ADDXY(JP)
   17 QXY(IH)=QXY(IH)+SIN(S)*SUBXY(JP)
      DO 16 I=1,LLAG
      CXY (I)=CXY(I)*CONST *.5
   16 QXY( I)=QXY(I)*CONST*.5
      SCXY(1)=.54*CXY(1)+.46*CXY(2)
      SCXY(LLAG)=.54*CXY(LLAG)+.46*CXY(LLAG-1)
      SQXY(1)=.54*QXY(1)+.46*QXY(2)
      SQXY(LLAG)=.54*QXY(LLAG)+.46*QXY(LLAG-1)
      KK=LLAG-1
      DO 18 J=2,KK
      SCXY(J)=.54*CXY(J)+.23*(CXY(J+1)+CXY(J-1))
   18 SQXY(J)=.54*QXY(J)+.23*(QXY(J+1)+QXY(J-1))
      DO 19 J=1,LLAG
   19 AMXY(J)=SQRT(SCXY(J)**2+SQXY(J)**2)
      CALL HONG
      WRITE (6,1005)UNIT,NNX,NNY,NNX,NNY,NNX,NNY,NNX,NNY
 1005 FORMAT(12H1  FREQUENCY,6X,11HCO-SPECTRUM,6X,19HQUADRATURE SPECTRUM
     1,5X,27HAMPLITUDE OF CROSS-SPECTRUM,4X,23HPHASE OF CROSS-SPECTRUM/1
     2X,8H(CYCLES/A6,1H),6X,2HOF,19X,2HOF,26X,2HOF,27X,2HOF/16X,6HSERIES
     3,I3,4H AND,I3,5X,6HSERIES,I3,4H AND,I3,12X,6HSERIES,I3,4H AND,I3,1
     43X,6HSERIES,I3,4H AND,I3,2H *//)
      ANY=-PMD
      DO 22 I=1,LLAG
      ANY=ANY+PMD
   22 WRITE (6,1006)ANY,SCXY(I),SQXY(I),AMXY(I),PHASXY(I)
 1006 FORMAT(3X,F8.3,5X,F14.7,7X,F14.7,14X,F14.7,15X,F14.7)
      WRITE (6,543)
  543 FORMAT(1H0,10X,37H* PHASE GIVEN IN FRACTION OF A CIRCLE)
      IF(IPRNT2(KIT+1).EQ.IB)GO TO 204
 201  IF(IPRNT2(KIT+1).NE.IS)GO TO 208
 204  DO 205 I=1,LLAG
  205 SCXY(I)=ALOG(AMXY(I))
      WRITE (6,1010)NNX,NNY,UNIT
 1010 FORMAT(1H1, 8X,46HGRAPH OF AMPLITUDE OF CROSS-SPECTRUM OF SERIES,I
     13,4H ANDI3,42H (IN LOG SCALE) AGAINST FREQUENCY (CYCLES/A6,1H)//)
      CALL PLUG (SCXY)
      WRITE (6,1011)NNX,NNY,UNIT
 1011 FORMAT(1H1,3X,33HPHASE OF CROSS-SPECTRUM OF SERIES,I3,4H AND,I3,63
     1H -- PLOTTED IN FRACTIONS OF A CIRCLE AGAINST FREQUENCY (CYCLES/,A
     26,1H)//)
      CALL KONG
 208  IF(ICOH.NE.IYES)GO TO 777
 215  IF(IPOW.EQ.IYES)GO TO 216
 210  WRITE (6,1015)NNX,NNY,KDUM
 1015 FORMAT(1H1,10X,5H*****/15X,29HTHE POWER SPECTRUM OF SERIES I3,5H A
     1ND I3,54H ARE NOT COMPUTED ACCORDING TO THIS SELECTION CARD NO. I3
     2/15X,68HTHE COHERENCE SQUARE AND THE TRANSFER FUNCTIONS CANNOT BE 
     3CALCULATED/ 15X,5H*****)
      GO TO 777
C
  216 ASSIGN 224 TO KADD1
      ASSIGN 250 TO KADD2
      SMAX=-9999999.0
      SCXMAX=-9999999.0
      SQXMAX=-9999999.0
      DO 220 I=1,LLAG
      IF (SPX(I))212,214,212
  212 RX(I)=BLANK
      SQXY(I)=AMXY(I)/SPX(I)
      SQXMAX=AMAX1(SQXMAX,SQXY(I))
  213 IF(SPY(I))217,218,217
  217 CRO(I)=BLANK
      COXYSQ(I)=(AMXY(I)**2)/(SPX(I)*SPY(I))
      IF(1.1-COXYSQ(I))2175,2180,2180
 2175 CRO(I)=DBLSTR
      ASSIGN 240 TO KADD2
      GO TO 221
C
 2180 SMAX=AMAX1(SMAX,COXYSQ(I))
      GO TO 221
C
  214 RX(I)=STAR
      ASSIGN 223 TO KADD1
      COXYSQ(I)=0.0
      CRO(I)=STAR
      SQXY(I)=0.0
  219 IF(SPY(I))221,222,221
  221 PX(I)=BLANK
      SCXY(I)=AMXY(I)/SPY(I)
      SCXMAX=AMAX1(SCXMAX,SCXY(I))
      GO TO 220
C
  218 COXYSQ(I)=0.0
      CRO(I)=STAR
  222 PX(I)=STAR
      SCXY(I)=0.0
      ASSIGN 223 TO KADD1
  220 CONTINUE
      WRITE (6,1012)UNIT,NNX,NNY,NNX,NNY,NNY,NNX,NNX,NNY
 1012 FORMAT(13H1   FREQUENCY7X16HCOHERENCE SQUARE,11X,12HAMPLITUDE OF,1
     13X,12HAMPLITUDE OF,15X,8HPHASE OF/9H (CYCLES/A6,1H),11X,2HOF,15X,3
     2(18HTRANSFER  FUNCTION7X)/20X,6HSERIES,I3,4H AND,I3,11X2(4HFROMI3,
     33H TO,I3,11X),4HFROMI3,3H TO,I3,2H *//)
      ANY=-PMD
      DO 230 I=1,LLAG
      ANY=ANY+PMD
      WRITE (6,1016)ANY,COXYSQ(I),CRO(I),SQXY(I),RX(I),SCXY(I),PX(I),PHA
     1SXY(I)
 1016 FORMAT(3X,F9.4,9X,F14.7,3(A2,9X,F14.7))
      IF(CRO(I).EQ.DBLSTR)GO TO 228
  227 IF(CRO(I).NE.STAR)GO TO 229
  228 COXYSQ(I)=SMAX
  229 IF(RX(I).EQ.STAR)GO TO 2292
 2291 IF(PX(I).EQ.STAR)GO TO 2294
	GO TO 230
 2292 SQXY(I)=SQXMAX
C
      GO TO 2291
 2294 SCXY(I)=SCXMAX
  230 CONTINUE
      WRITE (6,543)
      GO TO KADD2,(240,250)
  240 WRITE (6,1017)
 1017 FORMAT(1H0,10X,84H** INDICATES THAT THIS VALUE IS TOO HIGH DUE TO 
     1SAMPLING ERROR. IT WILL BE SET EQUAL/14X,71HTO THE MAXIMUM VALUE O
     2F THE REMAINING COHERENCES FOR PLOTTING PURPOSES.)
  250 GO TO KADD1,(223,224)
  223 WRITE (6,1018)
 1018 FORMAT(1H0,10X,91HX INDICATES THIS VALUE IS NOT COMPUTABLE DUE TO 
     1A NEGATIVE OR ZERO POWER SPECTRAL ESTIMATE./13X,82HIT WILL BE SET 
     2EQUAL TO THE MAXIMUM OF THE REMAINING VALUES FOR PLOTTING PURPOSES
     3.)
 224  IF(IPRNT3(KIT+1).EQ.IB)GO TO 232
 231  IF(IPRNT3(KIT+1).NE.IH)GO TO 233
 232  WRITE (6,1013)NNX,NNY,UNIT
 1013 FORMAT(1H1,16X,40HGRAPH OF COHERENCE SQUARE BETWEEN SERIESI3,4H AN
     1DI3,27H AGAINST FREQUENCY (CYCLES/A6,1H)//)
      CALL PLUG (COXYSQ)
 233  IF(IPRNT3(KIT+1).EQ.IB)GO TO 235
 234  IF(IPRNT3(KIT+1).NE.IT)GO TO 777
 235  DO 225 I=1,LLAG
      SQXY(I)=ALOG(SQXY(I))
  225 SCXY(I)=ALOG(SCXY(I))
      WRITE (6,1014)NNX,NNY,UNIT
 1014 FORMAT(1H1, 5X,51HGRAPH OF AMPLITUDE OF TRANSFER FUNCTION FROM SER
     1IESI3,3H TOI3,42H (IN LOG SCALE) AGAINST FREQUENCY (CYCLES/A6,1H)/
     2/)
      CALL PLUG (SQXY)
      WRITE (6,1014)NNY,NNX,UNIT
      CALL PLUG (SCXY)
  777 RETURN
C
      END
CHONG    SUBROUTINE HONG FOR BMD02T                  APRIL 15, 1963
      SUBROUTINE HONG
      DIMENSION Z(2000),DUMB(20),FMT(36),NSELL(20),RX(30),PX(30),
     1 LISA(20),ADDXY(30),SUBXY(30),CXY(30),QXY(30),SCXY(30),
     2 SQXY(30),AMXY(30),PHASXY(30),RXY(30),RYX(30),YP( 4),
     3 CRO(60),COXYSQ(30),IPRNT1(20),IPRNT2(20),IPRNT3(20)
      DOUBLE PRECISION UNIT
      INTEGER SYM(4)
      COMMON UNIT
      COMMON  RXY    , RYX    , RX     , SUBXY
      COMMON  Z      , DUMB   , FMT    , NSELL  , LISA   , SCXY
      COMMON  SQXY   , YP              , CRO    , COXYSQ , IORDAT
      COMMON  IPLDAT , IDETRN , IPREWT , THETA  , KSER   , NPOINT
      COMMON  LAG    , NSEL   , IPLOT  , DELTAT          , INFORM
      COMMON  KVR    , IPRNT1 , IPRNT2 , IPRNT3 , IPOW   , ICCS
      COMMON  ICOH   , PI     , FNPT   , CONST  , PMD    , FLAG
      COMMON  FLLAG  , LLAG   , KIT    , KDUM   , MISTAK , MX
      COMMON  KX     , Q      , KPOINT
      COMMON SYM
      EQUIVALENCE (RXY,CXY,AMXY),(RYX,QXY,PHASXY),(RX,ADDXY),(PX,SUBXY)
C
C     PI2 IS EQUAL TO 2 TIMES PI
C
      PI2=6.2831853071
C
      DO 10 I=1,LLAG
      AB=ABS(SQXY(I)/SCXY(I))
      PHI=ATAN(AB)
      IF( SCXY(I)) 11,12,13
   11 IF(SQXY(I)) 17,30,18
   17 PHASXY(I)=PI+PHI
      GO TO 10
C
 30   PHASXY(I)=PI
      GO TO 10
C
   18 PHASXY(I)=PI-PHI
      GO TO 10
C
 12   IF(SQXY(I))35,15,40
 35   PHASXY(I)=1.5*PI
      GO TO 10
C
   15 PHASXY(I)=0.0
      GO TO 10
C
C     STATEMENT 40 SETS PHASXY(I) EQUAL TO PI DIVIDED BY 2
C
   40 PHASXY(I)=1.5707963268
      GO TO 10
C
   13 IF (SQXY(I)) 14,15,16
 14   PHASXY(I)=PI2-PHI
      GO TO 10
C
   16 PHASXY(I)=PHI
   10 CONTINUE
      DO 100 I=1,LLAG
  100 PHASXY(I)=PHASXY(I)/PI2
      RETURN
C
      END
CKONG    SUBROUTINE KONG FOR BMD02T               FEBRUARY 17, 1964
      SUBROUTINE KONG
      INTEGER Q000CT,Q001CT,Q002CT,Q003CT
      DIMENSION Z(2000),DUMB(20),FMT(36),NSELL(20),RX(30),PX(30),
     1 LISA(20),ADDXY(30),SUBXY(30),CXY(30),QXY(30),SCXY(30),
     2 SQXY(30),AMXY(30),PHASXY(30),RXY(30),RYX(30),YP( 4),
     3 CRO(60),COXYSQ(30),IPRNT1(20),IPRNT2(20),IPRNT3(20)
      DOUBLE PRECISION UNIT
      INTEGER SYM(4)
      COMMON UNIT
      COMMON  RXY    , RYX    , RX     , SUBXY
      COMMON  Z      , DUMB   , FMT    , NSELL  , LISA   , SCXY
      COMMON  SQXY   , YP              , CRO    , COXYSQ , IORDAT
      COMMON  IPLDAT , IDETRN , IPREWT , THETA  , KSER   , NPOINT
      COMMON  LAG    , NSEL   , IPLOT  , DELTAT          , INFORM
      COMMON  KVR    , IPRNT1 , IPRNT2 , IPRNT3 , IPOW   , ICCS
      COMMON  ICOH   , PI     , FNPT   , CONST  , PMD    , FLAG
      COMMON  FLLAG  , LLAG   , KIT    , KDUM   , MISTAK , MX
      COMMON  KX     , Q      , KPOINT
      COMMON SYM
      EQUIVALENCE (RXY,CXY,AMXY),(RYX,QXY,PHASXY),(RX,ADDXY),(PX,SUBXY)
C
C
      DATA Q000CT/1H1/
      SYM(3)=Q000CT
      SYM(4)=SYM(3)
C
      IF(LLAG-35) 10,10,15
   10 NTIMES=2
      GO TO 30
   15 NTIMES=1
   30 TIMES=NTIMES+1
      DX=PMD/TIMES
      XMAX=1.0/(2.0*DELTAT)
      YMAX=1.5
      YMIN=-.5
      YP(3)=0.0
      YP(4)=1.0
      BNY=-PMD
      DO 100 I=1,LLAG
      BNY=BNY+PMD
      ANY=BNY
      FACE=PHASXY(I)+1.0
      IF( FACE-1.5) 105,105,110
  110 FACE=FACE-2.0
  105 YP(1)=PHASXY(I)
      YP(2)=FACE
      DATA Q001CT/2HX0/
      SYM(1)=Q001CT
      DATA Q002CT/2H*0/
      SYM(2)=Q002CT
      CALL PLOTR (ANY,0.0,XMAX,YP,SYM,YMIN,YMAX,4,-1)
      DATA Q003CT/2H 0/
      SYM(1)=Q003CT
      SYM(2)=SYM(1)
  125 DO 200 J=1,NTIMES
      ANY=ANY+DX
  200 CALL PLOTR (ANY,0.0,XMAX,YP,SYM,YMIN,YMAX, 4,-1)
  100 CONTINUE
      CALL PLOTR (ANY,0.0,XMAX,YP,SYM,YMIN,YMAX,-1,-1)
      RETURN
C
      END
CPLUG         SUBROUTINE PLUG FOR BMD02T              JULY 22, 1964
      SUBROUTINE PLUG (W)
      INTEGER Q000CT
      DIMENSION Z(2000),DUMB(20),FMT(36),NSELL(20),RX(30),PX(30),
     1 LISA(20),ADDXY(30),SUBXY(30),CXY(30),QXY(30),SCXY(30),
     2 SQXY(30),AMXY(30),PHASXY(30),RXY(30),RYX(30),YP( 4),
     3 CRO(60),COXYSQ(30),IPRNT1(20),IPRNT2(20)
      DIMENSION W(30),IPRNT3(20)
      DOUBLE PRECISION UNIT
      INTEGER SYM(4)
      COMMON UNIT
      COMMON  RXY    , RYX    , RX     , SUBXY
      COMMON  Z      , DUMB   , FMT    , NSELL  , LISA   , SCXY
      COMMON  SQXY   , YP              , CRO    , COXYSQ , IORDAT
      COMMON  IPLDAT , IDETRN , IPREWT , THETA  , KSER   , NPOINT
      COMMON  LAG    , NSEL   , IPLOT  , DELTAT          , INFORM
      COMMON  KVR    , IPRNT1 , IPRNT2 , IPRNT3 , IPOW   , ICCS
      COMMON  ICOH   , PI     , FNPT   , CONST  , PMD    , FLAG
      COMMON  FLLAG  , LLAG   , KIT    , KDUM   , MISTAK , MX
      COMMON  KX     , Q      , KPOINT
      COMMON SYM
      EQUIVALENCE (RXY,CXY,AMXY),(RYX,QXY,PHASXY),(RX,ADDXY),(PX,SUBXY)
C
C
      DATA Q000CT/2H*0/
      SYM(1)=Q000CT
C
      YMIN=9999999.0
      YMAX=-9999999.0
      DO 450 I=1,LLAG
      YMIN=AMIN1(YMIN,W(I))
  450 YMAX=AMAX1(YMAX,W(I))
      IF( LLAG-50) 455,460,460
  455 KA=1
      KB=0
      GO TO 470
C
 460  KA=-1
      KB=-1
 470  XMAX=0.5/DELTAT
 615  ANY=-PMD
      DO 630 I=1,LLAG
      ANY=ANY+PMD
      YP(1)=W(I)
  630 CALL PLOTR (ANY,0.00,XMAX,YP,SYM,YMIN,YMAX,1,KA)
      CALL PLOTR (ANY,0.00,XMAX,YP,SYM,YMIN,YMAX,KB,KA)
      RETURN
C
      END
CPOWER        SUBROUTINE POWER FOR BMD02T            APRIL 20, 1966
      SUBROUTINE POWER (NNX,X,XBAR,XALPHA,SPX)
C
      DIMENSION Z(2000),DUMB(20),FMT(36),NSELL(20),RX(30),PX(30),
     1 LISA(20),ADDXY(30),SUBXY(30),CXY(30),QXY(30),SCXY(30),
     2 SQXY(30),AMXY(30),PHASXY(30),RXY(30),RYX(30),YP( 4),
     3 CRO(60),COXYSQ(30),IPRNT1(20),IPRNT2(20)
      INTEGER SYM(4)
      DOUBLE PRECISION UNIT
      DIMENSION X(100),SPX(30),IPRNT3(20)
      COMMON UNIT
      COMMON  RXY    , RYX    , RX     , SUBXY
      COMMON  Z      , DUMB   , FMT    , NSELL  , LISA   , SCXY
      COMMON  SQXY   , YP              , CRO    , COXYSQ , IORDAT
      COMMON  IPLDAT , IDETRN , IPREWT , THETA  , KSER   , NPOINT
      COMMON  LAG    , NSEL   , IPLOT  , DELTAT          , INFORM
      COMMON  KVR    , IPRNT1 , IPRNT2 , IPRNT3 , IPOW   , ICCS
      COMMON  ICOH   , PI     , FNPT   , CONST  , PMD    , FLAG
      COMMON  FLLAG  , LLAG   , KIT    , KDUM   , MISTAK , MX
      COMMON  KX     , Q      , KPOINT
      COMMON SYM
      AMDA(Q000FL,Q001FL,Q002FL)=((Q000FL*Q001FL)**2)*(1.0-2.0*(Q002FL/Q
     1000FL)-2.0*(Q002FL/Q000FL)**2)
C
C
      EQUIVALENCE (RXY,CXY,AMXY),(RYX,QXY,PHASXY),(RX,ADDXY),(PX,SUBXY),
     1(YES,IYES),(IA,A),(IB,B),(IP,P)
C
      DATA A,B,BLANK,P,STAR,YES/1HA,1HB,1H ,1HP,1H*,3HYES/
      INTEGER I005CT
      DATA I005CT/2H*0/
      SYM(1)=I005CT
C
      DO 301 I=1,LLAG
      SX=0.0
      M=I-1
      NX=NPOINT-M
      FNX=NX
      DO 305 J=1,NX
      K=M+J
  305 SX=SX+X(J)*X(K)
  301 RX(I)=SX/FNX
 303  WRITE (6,1301)UNIT,NNX
 1301 FORMAT(1H1,21X,3HLAG,30X,14HAUTOCOVARIANCE/20X,1H(,A6,1H),28X,9HOF
     1 SERIES,I3//)
      TIME=-DELTAT
      DO 320 I=1,LLAG
      TIME=TIME+DELTAT
 320  WRITE (6,1302)TIME,RX(I)
 1302 FORMAT(19X,F9.4,26X,F13.6)
      TIME1=TIME
      IF(IPRNT1(KIT+1).EQ.IB)GO TO 325
 321  IF(IPRNT1(KIT+1).NE.IA)GO TO 615
  325 RXMIN=10**10
      RXMAX=-10**10
      WRITE (6,1303)NNX,TIME1,UNIT
 1303 FORMAT(1H1,9X,42HGRAPH OF AUTOCOVARIANCE FUNCTION OF SERIES,I3,37H
     1 PLOTTED AGAINST TIME UP TO A LAG OF F9.4,1X,A6//)
      DO 335 I=1,LLAG
      RXMIN=AMIN1(RXMIN,RX(I))
  335 RXMAX=AMAX1(RXMAX,RX(I))
      ANY=-DELTAT
      IF(LLAG-50) 600,605,605
 600  K1=0
      K2=1
      GO TO 606
C
 605  K1=-1
      K2=K1
 606  DO 610 I=1,LLAG
      ANY=ANY+DELTAT
      YP(1)=RX(I)
  610 CALL PLOTR(ANY,0.0,TIME1,YP,SYM,RXMIN,RXMAX,1,K2)
      CALL PLOTR(ANY,0.0,TIME1,YP,SYM,RXMIN,RXMAX,K1,K2)
  615 SX=RX(1)
      RX(1)=0.5*RX(1)
      RX(LLAG)=0.5*RX(LLAG)
      S1=-Q
      DO 401 IH=1,LLAG
      PX(IH)=0.0
      S1=S1+Q
      S=-S1
      DO 402 JP=1,LLAG
      S=S+S1
  402 PX(IH)=PX(IH)+RX(JP)*COS(S)
  401 PX(IH)=PX(IH)*CONST
      SPX(1)=.54*PX(1)+.46*PX(2)
      SPX(LLAG)=.54*PX(LLAG)+.46*PX(LLAG-1)
      KK=LLAG-1
      DO 415 J=2,KK
  415 SPX(J)=.54*PX(J)+.23*(PX(J+1)+PX(J-1))
      OMEGA=-Q
      THET1=1.0+(THETA*THETA)
      THET2=2.0*THETA
      IF(IPREWT.NE.IYES)GO TO 420
  425 DO 100 I=1,LLAG
      OMEGA=OMEGA+Q
      AA=THET1-THET2*COS(OMEGA)
  100 SPX(I)=SPX(I)/AA
 420  WRITE (6,1401)UNIT,NNX
 1401 FORMAT(1H119X9HFREQUENCY22X24HPOWER SPECTRAL ESTIMATES/17X,8H(CYCL
     1ES/A6,1H),25X, 9HOF SERIES I3//)
      ANY=-PMD
      ASSIGN 442 TO KADD1
      ASSIGN 433 TO KADD2
      RXMAX=-(0.5*(SPX(1)+SPX(LLAG)))
      DO 435 I=1,LLAG
      ANY=ANY+PMD
      IF(SPX(I))421,422,422
  421 RXY(I)=STAR
      GO TO 423
  422 RXY(I)=BLANK
  423 WRITE (6,1402)ANY,SPX(I),RXY(I)
 1402 FORMAT(21X,F8.3,26X,F14.7,A1)
      RXMAX=RXMAX+SPX(I)
      IF(RXY(I).EQ.BLANK)GO TO 435
  432 GO TO KADD2,(433,434)
  433 SPMAX=-9999999.0
      DO 4335 J=1,LLAG
      SPMAX=AMAX1(SPMAX,SPX(J))
 4335 CONTINUE
      ASSIGN 440 TO KADD1
      ASSIGN 434 TO KADD2
  434 SPX(I)=SPMAX
  435 CONTINUE
      GO TO KADD1,(440,442)
  440 WRITE (6,1405)
 1405 FORMAT(1H0, 6X,88H* THIS ESTIMATE IS NEGATIVE INDICATING SOME LEAK
     1AGE. IT WILL BE SET EQUAL TO THE MAXIMUM/ 9X,85HVALUE OF THE REMAI
     2NING ESTIMATES FOR FUTURE PLOTS AND EQUAL TO ZERO FOR CALCULATIONS
     3.)
  442 SPMAX=(Q*RXMAX)/DELTAT
      ANY=SPMAX-SX
      WRITE (6,1406)SPMAX,SX,ANY
 1406 FORMAT(1H0,6X,44HTHE CHECK SUM OF POWER SPECTRAL ESTIMATES IS,F14.
     17,14H AND SHOULD BE,F14.7/7X,17HTHE DIFFERENCE IS,F14.7)
      IF(IPRNT1(KIT+1).EQ.IB)GO TO 448
 443  IF(IPRNT1(KIT+1).NE.IP)GO TO 476
 448  IF(IPLOT) 450,460,450
 450  WRITE (6,1404)NNX,UNIT
 1404 FORMAT (1H1,18X,47HGRAPH OF THE POWER SPECTRAL ESTIMATES OF SERIES
     1I3,27H AGAINST FREQUENCY (CYCLES/A6,1H)//)
      CALL PLUG (SPX)
      IF(IPLOT) 476,460,460
 460  WRITE (6,1403)NNX,UNIT
 1403 FORMAT(1H1,10X,34HPOWER SPECTRAL ESTIMATES OF SERIESI3,6X49HPLOTTE
     1D IN A LOG SCALE AGAINST FREQUENCY (CYCLES/A6,1H)//)
      DO 475 I=1,LLAG
  475 RX(I)=ALOG(SPX(I))
      CALL PLUG (RX)
  476 DO 471 I=1,LLAG
      IF(RXY(I).EQ.BLANK)GO TO 471
  472 SPX(I)=0.0
  471 CONTINUE
  470 RETURN
C
      END
CSMEAN   SUBROUTINE SMEAN FOR BMD02T                 APRIL 20, 1966
      SUBROUTINE SMEAN (NNX,X,XBAR,XALPHA)
C
      DIMENSION Z(2000),DUMB(20),FMT(36),NSELL(20),RX(30),PX(30),
     1 LISA(20),ADDXY(30),SUBXY(30),CXY(30),QXY(30),SCXY(30),
     2 SQXY(30),AMXY(30),PHASXY(30),RXY(30),RYX(30),YP( 4),
     3 CRO(60),COXYSQ(30),IPRNT1(20),IPRNT2(20)
      DIMENSION X(100),IPRNT3(20)
      INTEGER  SYM(4), Q001CT
      DOUBLE PRECISION  UNIT, ACUMX(2)
      COMMON UNIT
      COMMON  RXY    , RYX    , RX     , SUBXY
      COMMON  Z      , DUMB   , FMT    , NSELL  , LISA   , SCXY
      COMMON  SQXY   , YP              , CRO    , COXYSQ , IORDAT
      COMMON  IPLDAT , IDETRN , IPREWT , THETA  , KSER   , NPOINT
      COMMON  LAG    , NSEL   , IPLOT  , DELTAT          , INFORM
      COMMON  KVR    , IPRNT1 , IPRNT2 , IPRNT3 , IPOW   , ICCS
      COMMON  ICOH   , PI     , FNPT   , CONST  , PMD    , FLAG
      COMMON  FLLAG  , LLAG   , KIT    , KDUM   , MISTAK , MX
      COMMON  KX     , Q      , KPOINT
      COMMON SYM
      EQUIVALENCE (RXY,CXY,AMXY),(RYX,QXY,PHASXY),(RX,ADDXY),(PX,SUBXY),
     1(YES,IYES)
C
      DATA IYES/3HYES/
      DATA Q001CT/2H*0/
      SYM(1)=Q001CT
C
      NPOINT=KPOINT
      IF (INFORM) 100,105,105
  100  DO 110 I=1,NPOINT
      ND=I+(NNX-1)*NPOINT
 110  X(I) = Z(ND)
      GO TO 135
  105 ISX=NPOINT*(NNX-1)
      DO 115 I=1,NPOINT
      K=I+ISX
  115 X(I)=Z(K)
 135  IF(IORDAT.NE.IYES)GO TO 140
 130  WRITE (6,1102)NNX
 1102 FORMAT(1H15X24HORIGINAL DATA OF SERIES I3//)
      WRITE (6,1104)(X(I),I=1,NPOINT)
 1104 FORMAT(10F11.5)
 140  IF(IPLDAT.NE.IYES)GO TO 141
  200 HU=-10**10
      SM=10**10
      DO 210 I=1,NPOINT
      HU=AMAX1(HU,X(I))
  210 SM=AMIN1(SM,X(I))
      WRITE (6,2001)NNX
 2001 FORMAT(1H1,10X,22HGRAPH OF INPUT SERIES I3///)
      FI=0.0
      DO 215 I=1,NPOINT
      FI=FI+1.0
      YP(1)=X(I)
      CALL PLOTR (FI,1.0,FNPT,YP,SYM,SM,HU,1,-1)
  215 CONTINUE
      CALL PLOTR (FI,1.0,FNPT,YP,SYM,SM,HU,-1,-1)
  141 MISTAK=0
  150 LL=KIT+1
      IF(NNX-LISA(LL)) 145,155,145
  155 CALL TRANS(NNX,X)
      IF (MISTAK) 160,145,160
 145  IF(IPREWT.NE.IYES)GO TO 305
 340  NT = NPOINT - 1
      DO 350 I=1,NT
 350  X(I)=X(I+1)-THETA*X(I)
      NPOINT=NT
 305  FNPT=NPOINT
      ACUMX(1) = 0.
      ACUMX(2) = 0.
      DO 330 I1=1,NPOINT
  330 ACUMX(1) = ACUMX(1) + X(I1)
      XBAR = ACUMX(1) / FNPT
      IF(IDETRN.NE.IYES)GO TO 360
 300  ACUMX(1) = 0.
      ACUMX(2) = 0.
      TBAR = NPOINT - 1
      TBAR = TBAR / 2.
      DENOM = (FNPT-1.) * FNPT * (FNPT+1.) / 6.
      DO 310 I1=1,NPOINT
      ACUMX(2) = 2 * I1 - NPOINT + 1
      ACUMX(2) = (X(I1) - XBAR) * ACUMX(2)
  310 ACUMX(1) = ACUMX(1) + ACUMX(2)
      XALPHA = ACUMX(1) / DENOM
      XBETA = XBAR - XALPHA * TBAR
      DO 320 I1=1,NPOINT
      FI1 = I1 - 1
 320  X(I1) = X(I1) - FI1 * XALPHA - XBETA
      GO TO 160
 360  DO 370 I1=1,NPOINT
 370  X(I1) = X(I1) - XBAR
      FNPT=NPOINT
  160 RETURN
C
      END
CTPWD         SUBROUTINE TPWD FOR BMDXXX SERIES
      SUBROUTINE TPWD(NT1,NT2)
      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)
      STOP
 49   FORMAT(25H ERROR ON TAPE ASSIGNMENT)
      END
CTRANS        SUBROUTINE TRANS FOR BMD02T             JUNE 2, 1964
      SUBROUTINE TRANS(NNX,X)
C
      DIMENSION Z(2000),DUMB(20),FMT(36),NSELL(20),RX(30),PX(30),
     1 LISA(20),ADDXY(30),SUBXY(30),CXY(30),QXY(30),SCXY(30),
     2 SQXY(30),AMXY(30),PHASXY(30),RXY(30),RYX(30),YP( 4),
     3 CRO(60),COXYSQ(30),CON( 8),IBIN( 8),IPRNT1(20),IPRNT2(20)
      INTEGER SYM(4)
      DOUBLE PRECISION UNIT
      DOUBLE PRECISION JOB,SPC
       DATA SPC/6HSPECTG/
      DIMENSION X(100),IPRNT3(20)
      COMMON UNIT
      COMMON  RXY    , RYX    , RX     , SUBXY
      COMMON  Z      , DUMB   , FMT    , NSELL  , LISA   , SCXY
      COMMON  SQXY   , YP              , CRO    , COXYSQ , IORDAT
      COMMON  IPLDAT , IDETRN , IPREWT , THETA  , KSER   , NPOINT
      COMMON  LAG    , NSEL   , IPLOT  , DELTAT          , INFORM
      COMMON  KVR    , IPRNT1 , IPRNT2 , IPRNT3 , IPOW   , ICCS
      COMMON  ICOH   , PI     , FNPT   , CONST  , PMD    , FLAG
      COMMON  FLLAG  , LLAG   , KIT    , KDUM   , MISTAK , MX
      COMMON  KX     , Q      , KPOINT
      COMMON SYM
      ASN(Q000FL)=ATAN(Q000FL/SQRT(1.0-Q000FL**2))
C
      EQUIVALENCE (RXY,CXY,AMXY),(RYX,QXY,PHASXY),(RX,ADDXY),(PX,SUBXY),
     1(SPC,ISPC)
C
C
      READ (5,1002)JOB,NTRAN,(IBIN(I),CON(I),I=1,NTRAN)
 1002 FORMAT(A6,I1,8(I2,F6.0))
      IF(JOB.NE. SPC)GO TO 700
 602  DO 500 I=1,NTRAN
      IF(IBIN(I)-17) 605,600,605
  605 JESUS=IBIN(I)
      IF(JESUS-6)610,905,610
C
  600 JESUS=6
  610 CC=CON(I)
      IF(JESUS*(JESUS-11)) 620,905,905
  620 DO 150 K=1,NPOINT
      GO TO (10,20,30,40,50,60,70,80,90,100) ,JESUS
   10 IF(X(K))200,150,14
C
   14 X(K)=SQRT(X(K))
      GO TO 150
C
   20 IF(X(K))200,22,23
   22 X(K)=1.0
      GO TO 150
C
   23 X(K)=SQRT(X(K))+SQRT(X(K)+1.0)
      GO TO 150
   30 IF(X(K))200,200,31
   31 X(K)=.434294481E+00*ALOG(X(K))
      GO TO 150
C
   40 X(K)=EXP(X(K))
      GO TO 150
C
   50 IF(X(K))200,150,53
C
   53 IF (X(K)-1.0) 54,55,200
  54  ARG =SQRT(X(K))
      X(K)=ASN(ARG)
      GO TO 150
C
   55 X(K)=PI/2.0
      GO TO 150
C
   60 IF(X(K))200,200,61
   61 X(K)=ALOG(X(K))
      GO TO 150
C
   70 IF(X(K))71,200,71
   71 X(K)=1.0/X(K)
      GO TO 150
C
   80 X(K)=X(K)+CC
      GO TO 150
C
   90 X(K)=X(K)*CC
      GO TO 150
C
  100 IF(X(K))200,200,101
  101 X(K)=X(K)**CC
      GO TO 150
  200 WRITE (6,1001)K,NNX,IBIN(I)
 1001 FORMAT(11H0DATA POINT I5,10H OF SERIESI3,
     157H VIOLATES THE RESTRICTION FOR TRANSGENERATION OF THE TYPEI3,
     252H. THE PROGRAM CONTINUES LEAVING THE VALUE UNCHANGED.)
  150 CONTINUE
  500 CONTINUE
  300 RETURN
  700 WRITE (6,1003)JOB
      WRITE (6,1006) JOB,NTRAN, (IBIN(I),CON(I),I=1,NTRAN)
 1006 FORMAT('0IMPROPER SPECTG CARD IS PRINTED BELOW.' / 1X,A6,I1,8(I2,
     1 F6.0))
  901 WRITE (6,1004)
      MISTAK=11
      GO TO 300
 1003 FORMAT(58H0CONTROL CARD ERROR. PROGRAM EXPECTED A SPECTG CARD BUT 
     1A A6,16H CARD WAS FOUND.)
 1004 FORMAT(55H0PROGRAM WILL GO TO THE NEXT SELECTION OR PROBLEM CARD.)
  905 WRITE (6,1005)
      GO TO 901
 1005 FORMAT(40H0ILLEGAL TRANSGENERATION CODE SPECIFIED.)
      END
CVFCHCK    SUBROUTINE TO CHECK FOR PROPER NUMBER OF VARIABLE FORMAT CRDS
      SUBROUTINE VFCHCK(NVF)
      IF(NVF)10,10,20
 10   WRITE (6,4000)
      NVF=1
 50   RETURN
C
   20 IF (NVF.LE.2) GO TO 50
      GO TO 10
C
 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
     1IED, ASSUMED TO BE 1.)
      END
C             SUBROUTINE PLOTR (IBM 360)              AUGUST 13, 1966
C
C     'PLOTR' IS A UTILITY SUBPROGRAM FOR THE BMD... PROGRAMS WHICH
C     PLOTS EITHER SINGLE-LINE OR WHOLE-PAGE GRAPHS AND SETS UP
C     APPROPRIATE SCALING.  THE CALLING PARAMETERS ARE AS FOLLOWS -
C
C     X - THE VALUE OF THE INDEPENDENT VARIABLE
C     ZMIN - THE MINIMUM VALUE OF X FOR THIS PLOT
C     ZMAX - THE MAXIMUM VALUE OF X FOR THIS PLOT
C     Y - THE ARRAY CONTAINING THE VALUES OF UP TO 15 DEPENDENT VAR.'S
C     SYM - THE ARRAY CONTAINING THE SYMBOLS TO BE PLOTTED
C     WMIN - THE MINIMUM VALUE OF ALL Y'S FOR THIS PLOT
C     WMAX - THE MAXIMUM VALUE OF ALL Y'S FOR THIS PLOT
C     NC - THE NUMBER OF DEPENDENT VARIABLES
C               NC=-1 CLOSES A SINGLE-LINE PLOT
C               NC= 0 PRINTS AND CLOSES A WHOLE-PAGE PLOT
C     NP - THE CONTROL VARIABLE
C               NP=-1 PRINTS A SINGLE LINE
C               NP=0 OR NP=1 SETS UP A WHOLE-PAGE PLOT
C
C     THE PLOTTING ROUTINE MUST BE CALLED ONCE FOR EACH VALUE OF THE
C     INDEPENDENT VARIABLE THAT IS TO BE PLOTTED NO MATTER WHETHER IN
C     THE SINGLE-LINE OR WHOLE-PAGE MODE
C
      SUBROUTINE PLOTR(X,ZMIN,ZMAX,Y,SYM,WMIN,WMAX,NC,NP)
C
      DIMENSION Y(15),CLAB(12),GF(10),FMT(12),XY(51,101),SYM(15)
      INTEGER XY,BLANKS,SYM,SYMB
      DATA TC,TP,BLANKS/1H.,1H+,1H /
      DATA GF/            4H 1X,,4H 2X,,4H 3X,,4H 4X,,4H 5X,,4H 6X,,
     14H 7X,,4H 8X,,4H 9X,,4H 10X/
      DATA FMT/'(17X',' ','5(F1','2.3,','8X)/','7X, ',' ','4(F1','2.3,',
     1'8X),','F12.','3)  '/
C
 100   FORMAT(1H 6X5(F12.3,8X),F12.3/17X,5(F12.3,8X))
 101  FORMAT(1H F12.3,1X,103A1,F12.3)
  102 FORMAT (1H 13X,103A1)
 1000 FORMAT(1H  14X,101A1)
 1001 FORMAT(15X,20(5H+....),1H+)
C
      DATA NCC/2/
C    'NCC' ON THE INITIAL ENTRY TO PLOTR IS NON-ZERO BECAUSE OF THE DATA
C    STATEMENT ABOVE.
C
C    'NCC' IS 0 WHILE A PLOT IS BEING MADE.  IT IS 1 OR 2 AT OTHER TIMES
C
      IF(NCC) 50,48,50
C
C    THE VARIABLE 'KL' CONTROLS THE FUNCTIONING OF THE OPENING AND
C    CLOSING  SECTIONS OF PLOTR.  KL=0 INDICATES OPENING OF THE GRAPH,
C    KL=1 INDICATES CLOSING.
C
   50 KL=0
      CALL SCALE(WMIN,WMAX,100.0,JY,YMIN,YMAX,YIJ)
      YR=YMAX-YMIN
  230 J=JY
      IF(J*(J-10))204,201,201
C
C     THE FOLLOWING SECTION OPENS OR CLOSES A PLOT IN FIXED FORMAT
C     UNDER CONTROL OF KL
C
  201 IF(KL)220,220,231
C
  231 WRITE (6,1001)
      IF(KL)250,250,220
C
  220 CLAB(1)= YMIN
      DO 222 I=2,11
  222 CLAB(I)=CLAB(I-1)+YIJ
      WRITE (6,100)(CLAB(I),I=1,11,2),(CLAB(J),J=2,10,2)
      IF(KL)231,231,14
C
C     THE FOLLOWING SECTION OPENS OR CLOSES A PLOT IN A VARIABLE
C     FORMAT UNDER CONTROL OF KL AND JY FROM 'SCALE'
C
  204 IF(J-5)205,221,207
  207 J=J-5
  205 JYT=5-J
  221 CONTINUE
  226 FMT(2)=GF(JY)
      IF (KL) 225,225,227
C
  225 FMT(7)=GF(JY)
      TT=JY
      TT=TT*YIJ/10.0
      CLAB(1)= YMIN+TT
      DO 223 I=2,10
  223 CLAB(I)=CLAB(I-1) +YIJ
      WRITE (6,FMT)(CLAB(I),I=2,10,2),(CLAB(I),I=1,9 ,2)
      IF(KL)227,227,14
C
  227 IF(JY-5)208,209,208
  208 WRITE(6,1000)(TC,I=1,J    ),(TP,(TC,I=1,4),K=1,19),TP,(TC,I=1,JYT)
      IF (KL) 250,250,225
C
  209 WRITE (6,1001)
      IF (KL) 250,250,225
C
  250 CONTINUE
      NCC=0
      IC=0
      IF(NP)80,11,11
C
C    THIS SECTION PREPARES FOR A FULL PAGE PLOT BY FILLING IN XY WITH
C     BLANKS AND SETTING UP SCALING FOR THE INDEPENDENT VARIABLE 'X'
C
   11 DO 1 I=1,51
       DO 1 J=1,101
   1  XY(I,J)=BLANKS
      CALL SCALE (ZMIN,ZMAX,50.0,JX,XMIN,XMAX,XIJ)
      XR=XMAX-XMIN
      GO TO 48
C
C
C     ENTRY TO PLOTS CAN BE USED ONLY AFTER THE CALLING PARAMETERS
C     HAVE BEEN TRANSFERRED BY A CALL ON PLOTR.  THE CALL ON PLOTS
C     IS IDENTICAL WITH ENTRY TO PLOTR BUT IT ALLOWS THE PROGRAMMER TO
C     CALL THE PLOTTING ROUTINE WITHOUT HAVING TO INCLUDE THE PARAMETERS
C
C      ENTRY PLOTS
   48 IF(NC)52,13,49
   49 IF(NP)80,10,10
C    THE FOLLOWING SECTION SETS UP A FULL PAGE BUT DOES NO PRINTING.
C    THIS SECTION IS REACHED BY SPECIFYING NC POSITIVE AND NP POSITIVE.
C
   10 DO 9 N=1,NC
      SYMB=SYM(N)
      XDIFFR=XMAX-X
      IF(XDIFFR)105,106,106
  105 XDIFFR=0.0
  106 YDIFFR=YMAX-Y(N)
      IF(YDIFFR)107,108,108
  107 YDIFFR=0.0
  108 L=51.0-(50.0*XDIFFR)/XR+.5
      K=101.0-(100.0*YDIFFR)/YR+.5
      CALL FORM2(SYMB,XY(L,K))
 9     CONTINUE
      GO TO 15
C
C     THE FOLLOWING SECTION CONSTRUCTS AND PLOTS ONE LINE OF A MULTILINE
C     GRAPH.  LOCATION ALONG THE AXIS OF THE PAPER IS PRINTED AT EVERY
C     STEP.  THIS SECTION IS ACCESSIBLE BY SPECIFYING NC POSITIVE AND
C     NP NEGATIVE.
C
   80 DO 86 I=1,101
   86 XY(1,I)=BLANKS
       L=1
      DO 95 N=1,NC
      SYMB=SYM(N)
      YDIFFR=YMAX-Y(N)
      IF(YDIFFR)860,865,865
  860 YDIFFR=0.0
  865 K=101.0-(100.0*YDIFFR)/YR+.5
   95 CALL FORM2(SYMB,XY(L,K))
      IF(MOD(IC,5))97,96,97
   96 W=TP
      GO TO 98
   97 W=TC
   98 WRITE (6,101)X,W,(XY(1,N),N=1,101),W,X
      IC=IC+1
      GO TO 15
C
C    THIS SECTION PLOTS OUT THE PREVIOUSLY PREPARED WHOLE PAGE GRAPH.
C    IT PRINTS LOCATION ALONG THE PAPER'S AXIS EVERY FIFTH STEP.  THIS
C    SECTION IS ACCESSED BY SPECIFYING NC=0.
C
   13 M=6-JX
      LL=50+M
      T=JX
      IF(5-JX)131,131,135
  131 T=0.0
  135 RLAB=XMAX-(T*XIJ)/5.0
      W=TC
      K=52
      DO 31 L=M,LL
      K=K-1
      I=MOD(L,5)
      IF(I-1)2,3,2
    3 W=TP
      WRITE (6,101)RLAB,W,(XY(K,N),N=1,101),W,RLAB
      RLAB=RLAB-XIJ
      W=TC
      GO TO 31
    2 WRITE (6,102)W,(XY(K,N),N=1,101),W
   31 CONTINUE
C
   52 KL=1
      GO TO 230
C
   14 NCC=1
   15 RETURN
      END
C             SUBROUTINE SCALE FOR PLOTR              JUNE 21, 1966
C
C     SUBROUTINE 'SCALE' CALCULATES THE SCALING FOR 'PLOTR'
C
      SUBROUTINE SCALE(YMIN,YMAX,YINT,JY,TYMIN,TYMAX,YIJ)
      DIMENSION C(10)
      DATA C           /1.0,1.5,2.0,3.0,4.0,5.0,7.5,10.0,15.0,20.0/
        DATA TEST / 0.76293945E-05/
   50 YR=YMAX-YMIN
      TT=YR/YINT
      J = ALOG10(TT)+TEST
      E=10.0**J
      TT=TT/E
      I=0
      IF(TT-1.0+TEST)205,201,201
  205 TT=TT*10.0
      E=E/10.0
 201  I=I+1
      IF(9-I)1,2,2
    1 E=E*10.0
      I=1
    2 IF(TT-C(I))233,202,201
  233 YIJ=C(I)*E
      GO TO 203
  202 Y=YMIN/C(I)
      J=Y
      T=J
      IF(0.0001-ABS(T-Y))204,233,233
  204 YIJ=C(I+1)*E
  203 X=((YMAX+YMIN)/YIJ-YINT )/2.0+.00001
      K=X
      IF(K)235,240,240
  235 Y=K
      IF(X-Y)236,240,236
  236 K=K-1
  240 TYMIN=K
      TYMIN=YIJ*TYMIN
      TYMAX=TYMIN+YINT*YIJ
      IF (YMAX-TYMAX-TEST)11,11,201
   11 YIJJ=C(I)*E
      XT=((YMAX+YMIN)/YIJJ-YINT)/2.0+.00001
      KT=XT
      IF (KT) 1235,1240,1240
 1235 YT=KT
      IF (XT.NE.YT) KT=KT-1
 1240 TYMINT=KT
      TYMINT=YIJJ*TYMINT
      TYMAXT=TYMINT+YINT*YIJJ
      IF (YMAX-TYMAXT.GT.TEST) GO TO 10
      TYMIN=TYMINT
      TYMAX=TYMAXT
      YIJ=YIJJ
      K=KT
   10 TT=YINT/10.0
      JY=TT+.000001
      YIJ=YINT*(YIJ/10.0)
      J=TYMIN/ YIJ
      IF (K)242,241,241
  242 J=J-1
  241 J=J*JY+JY-K
      JY=J
      RETURN
      END
      SUBROUTINE FORM2(SYMB,XY)
C        SUBROUTINE FORM2 FOR PLOTR (IBM 360)         JUNE 21, 1966
      INTEGER XY,SYMB,BLANK,TEST(18)
       DATA BLANK/'  '/,TEST/'2 ','3 ','4 ','5 ','6 ',
     1'7 ','8 ','9 ','A ','B ','C ','D ','E ','F ','G ','H ','I ','/ '/
      IF(XY.EQ.BLANK)GO TO 50
      DO 30 I=1,17
      IF(XY.NE.TEST(I))GO TO 30
C        PUT IN NEXT SYMBOL OF ARRAY FOR MULTIPLE POINTS
      XY=TEST(I+1)
      GO TO 100
   30 CONTINUE
      IF(XY.EQ.TEST(18))GO TO 100
C        IF OTHER THAN CHARACTERS IN ARRAY TEST PUT IN CHARACTER 2.
      XY=TEST(1)
      GO TO 100
C        IF BLANK, PUT IN SYMBOL
   50 XY=SYMB
  100 RETURN
      END













































