CBD01T        AMPLITUDE AND PHASE ANALYSIS          AUGUST  9, 1965
C        THIS IS A SIFTED VERSION OF BMD01T ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
          DIMENSION X(350),FMT(36),YSQ(80),Y(80),SYM(15),YP(15),
     1A(350),B(350),XX(350),YY(350)
      COMMON  X      , FMT    , YSQ    , Y      , SYM    , YP
      COMMON  A      , B      , XX     , YY     , AMPL   , PHASE
      COMMON  TG     , FILTER , KPROB  , NDATA  , FTRI   , NTRI
      COMMON  M      , FM     , NPOINT , POINT  , WO     , ANS
      COMMON  PI     , H
C         DIMENSION X(800),FMT(120),YSQ(500),Y(500),SYM(15),YP(15),
C    1A(800),B(800),XX(800),YY(800)
C     COMMON X,FMT,YSQ,Y,SYM,YP,A,B,XX,YY,AMPL,PHASE,TG,FILTER,KPROB,
C    1NDATA,FTRI,NTRI,M,FM,NPOINT,POINT,WO,ANS,PI,H
C
  305 FORMAT(52H1BMD01T - AMPLITUDE AND PHASE ANALYSIS - VERSION OF ,
     115HAUGUST  9, 1965/
     2 40H HEALTH SCIENCES COMPUTING FACILITY,UCLA/
     314H PROBLEM CODE A6,/
     423H NUMBER OF DATA POINTS I5,///)
C
	DOUBLE PRECISION A123,B123,TODE,CODE
      DATA A123/'FINISH  '/,B123/'PROBLM  '/
      DATA Q002HL/4HYES /
      DATA Q003HL/4H*000/
      YES =(+Q002HL)
      NTAPE=5
	CALL USAGEB('BMD01T')
  304 FORMAT('0CONTROL CARD WAS INCORRECTLY ORDERED OR PUNCHED. IT IS '
     1 'PRINTED BELOW AND PROGRAM HAS BEEN TERMINATED.')
 5001 FORMAT(2A6,4A3,I3,I4,I3,I4,I3,F5.0,22X,2I2)
 5002 FORMAT (1H0,2A6,4A3,I3,I4,I3,I4,I3,F5.0,22X,2I2)
  500 READ (5,5001)TODE,CODE,XG,XDATA,ELTER,XOOL,KPROB,NDATA,NTRI,M,
     1NPOINT,WO,INFORM,KVR
      IF(TODE.EQ.B123) GO TO 20
      IF(TODE.EQ.A123) GO TO 303
  302 CONTINUE
      WRITE (6,3040)
 3040 FORMAT ( '0A ''PROBLM'' OR ''FINISH'' CARD WAS EXPECTED.' )
      WRITE (6,304)
      WRITE (6,5002) TODE,CODE,XG,XDATA,ELTER,XOOL,KPROB,NDATA,NTRI,M,
     1 NPOINT,WO,INFORM,KVR
  303 IF(NTAPE-5)372,372,371
  371 REWIND NTAPE
  372 WRITE(6,3722)
 3722 FORMAT (1H0 / / '0 FINISH CARD ENCOUNTERED.' )
      STOP
   20 CALL TPWD(INFORM,NTAPE)
      WRITE (6,305)CODE,NDATA
      IF(-NTRI)25,302,302
   25 FTRI=NTRI
      IF(M-350)255,255,302
  255 FM=M
      PI=3.14159265
      CRIT=(FM-2.0)/2.0
      IF(CRIT-FTRI)  11,11,12
   11 DO 90 I=1,6
   90 WRITE (6,2000)
      WRITE (6,1011)KPROB
 1011 FORMAT(10X,5H* * */10X,28HBANDWIDTH CHOSEN IN PROBLEM I3,
     1 50H LARGER THAN PI, PROGRAM GOES TO NEXT PROBLEM CARD/10X,5H* * *
     2)
      PI=-4000
   12 IF(KVR.GT.0.AND.KVR.LE.2) GO TO 210
      WRITE (6,4000)
 4000 FORMAT(1H0,23X,71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPEC
     XIFIED, ASSUMED TO BE 1.)
      KVR = 1
  210 KVR=KVR*18
      READ (5,5003)(FMT(I),I=1,KVR)
 5003 FORMAT(18A4)
      WRITE (6,5005) (FMT(I),I=1,KVR)
 5005 FORMAT ('0VARIABLE FORMAT USED IS PRINTED BELOW.'/(1H0, 18A4))
      IF((NDATA-49)*(NDATA-351))215,302,302
  215 MISTAK=0
      READ (NTAPE,FMT)(X(I),I=1,NDATA)
      IF(-PI)40,500,500
   40 IF(XG.NE.YES) GO TO 30
   45 CALL TRANS( MISTAK)
      IF( MISTAK) 500,30,500
   30 IF(XDATA.NE.YES) GO TO 60
   65 WRITE (6,2000)
      WRITE (6,1017)KPROB
 1017 FORMAT(10X,25HORIGINAL DATA OF PROBLEM I3//)
      WRITE (6,1018)(X(I),I=1,NDATA)
 1018 FORMAT(10F12.6)
   60 IF(ELTER.NE.YES) GO TO 126
   95 POINT=NPOINT
      IF(50-NPOINT)951,951,302
  951 IF(WO*(WO-PI))955,955,96
  955 H=PI/POINT
      G=(FTRI+1.0)*2.0*PI/FM
      BONE=WO-.5*G
      IF(BONE) 96,96,97
   96 WRITE (6,1019)KPROB
 1019 FORMAT(10X,49H*****THE MAIN LOBE OF THE FILTER OF THIS PROBLEM(I3,
     166H) IS NOT BETWEEN 0 AND PI,PROGRAM GOES TO THE NEXT PROBLEM, IF 
     2ANY)
      GO TO 500
   97 NONE=(BONE/H)+0.5
      WRITE (6,800)
  800 FORMAT(1H1)
      CALL LANAI (NONE,0.0,0)
      AA=ANS
  100 WRITE (6,2000)
 2000 FORMAT(1H0)
      ZERO=0.0
      WRITE (6,1002)ZERO,BONE,AA
 1002 FORMAT(10X,38HVALUE OF INTEGRAL OF F(X) SQUARE FROM F10.5,
     14H TO F10.5,3H = F10.5//)
      BTWO=WO+.5*G
      IF(BTWO-PI) 101,101,102
  102 WRITE (6,1019)KPROB
      GO TO 500
  101 NTWO=(G/H)+0.5
      CALL LANAI (NTWO,BONE,NONE)
      BB=ANS
  105 WRITE (6,2000)
      WRITE (6,1002)BONE,BTWO,BB
      NTHREE=((PI-BTWO)/H)+0.5
      KK=NONE+NTWO
      CALL LANAI(NTHREE,BTWO,KK)
      CC=ANS
  110 WRITE (6,2000)
      WRITE (6,1002)BTWO,PI,CC
      PC=(AA+CC)*100.0/(AA+BB+CC)
      WRITE (6,1006)
 1006 FORMAT(1H0,9X,10HRESOLUTION7X,16HNO. OF TRIANGLES9X,11HCENTERED AT
     18X,9HBANDWIDTH,10X,19HLEAKAGE IN PER CENT//)
      WRITE (6,1008)M,NTRI,WO,G,PC
 1008 FORMAT(13X,I3,18X,I3,15X,F11.8,8X,F10.7,12X,F12.8)
  115 WRITE (6,2000)
      WRITE (6,1005)
 1005 FORMAT(1H0,9X,34HTHESE ARE THE VALUES OF THE FILTER//)
      NPOINT=NONE+NTWO+NTHREE
      POINT=NPOINT-1
      WRITE (6,1007)(Y(I),I=1,NPOINT)
 1007 FORMAT(7F12.6)
      WRITE (6,1010)
 1010 FORMAT(1H1,20X,19HGRAPH OF THE FILTER//)
      YMAX=-10**10
      YMIN=10**10
      DO 120 I=1,NPOINT
      YMAX=AMAX1(Y(I),YMAX)
  120 YMIN=AMIN1(Y(I),YMIN)
      SYM(1)=(+Q003HL)
      DO 125 I=1,NPOINT
      FI=I-1
      XX(1)=FI*PI/POINT
      YP(1)=Y(I)
  125 CALL PLOTR (XX,0.0,PI,YP,SYM,YMIN,YMAX,1,-1)
      CALL PLOTR (XX,0.0,PI,YP,SYM,YMIN,YMAX,-1,-1)
  126 IF(XOOL.NE.YES) GO TO 131
  127 CALL COEVV
      BZERO=2.0*FTRI/FM
      WRITE (6,1003)
 1003 FORMAT(1H1,9X47HB IS THE COEFFICIENT OF COS J.OMEGA IN C(OMEGA)/
     1 10X,75HAND A IS THE COEFFICIENT OF SIN J.OMEGA IN S(OMEGA) IN SEC
     2TION FOUR STEP 2 ///)
      WRITE (6,1009)BZERO
 1009 FORMAT(10X16HA(   0) =    0.024X9HB(   0) =F11.5)
      DO 130 I=1,M
  130 WRITE (6,1004)I,A(I),I,B(I)
 1004 FORMAT (10X,2HA(I4,4H) = F10.5,20X,2HB(I4,4H) = F10.5)
  131 DO 140 K=1,NDATA
      XX(K)=0.0
      YY(K)=0.0
      DO 155 J=1,M
      JS=K+M+J
      KS=K+M-J+1
      IF((JS+1)-NDATA) 145,145,150
  145 YY(K)=YY(K)+A(J)*(X(JS)-X(KS))
  155 XX(K)=XX(K)+B(J)*(X(JS+1)+X(KS-1))
      IS=M+K+1
      XX(K)=BZERO*(X(IS)+X(IS-1))+XX(K)
  140 CONTINUE
  150 WRITE (6,1020)
 1020 FORMAT(1H1,31X,56HGRAPH OF THE FILTERED INPUT AT THE INDICATED DAT
     1A POINTS//)
      K=K-1
      YMAX=-10.0**10
      YMIN= 10.0**10
      DO 157 I=1,K
      YMAX=AMAX1(XX(I),YMAX)
  157 YMIN=AMIN1(XX(I),YMIN)
      IF(YMAX-YMIN-1.0)1575,1575,158
 1575 DIV=(YMAX-YMIN)/3.0
      YMAX=YMAX+DIV
      YMIN=YMIN-DIV
  158 FK=M
      FMAX=NDATA-M
      DO 159 I=1,K
      FK=FK+1.0
      YP(1)=XX(I)
  159 CALL PLOTR(FK,FM,FMAX,YP,SYM,YMIN,YMAX, 1,-1)
      CALL PLOTR(FK,FM,FMAX,YP,SYM,YMIN,YMAX,-1,-1)
      WRITE (6,1023)KPROB
 1023 FORMAT(1H1,20X,72HAVERAGE FREQUENCY, AMPLITUDE, PHASE AND FINITE M
     1OVING AVERAGES OF SERIES I4/32X,53H(FREQUENCY IN RADIANS/UNIT TIME
     2 AND PHASE IN RADIANS)//)
      WRITE (6,1025)
 1025 FORMAT(2X,4HDATA,7X,9HFREQUENCY,11X,9HAMPLITUDE,13X,5HPHASE,7X,38H
     1FINITE MOVING AVERAGES USING CONSTANTS/1X,5HPOINT,68X,4HB(J),16X,
     24HA(J)//)
      CALL PHAZE(XX(1),YY(1),Z1)
      J=M+1
      DO 165 I=2,K
      J=J+1
      AMPL =SQRT(XX(I)**2+YY(I)**2)
      CALL PHAZE(XX(I),YY(I),Z2)
      FREQ=Z2-Z1
      IF(FREQ)160,164,164
  160 DIV= FREQ+(2.0*PI)
      FREQ=-FREQ
      FREQ=AMIN1(DIV,FREQ)
      IF(FREQ-DIV)163,164,164
  163 FREQ=-FREQ
  164 WRITE (6,1024)J,FREQ,AMPL,Z2,XX(I),YY(I)
 1024 FORMAT(1H I5,4XF12.6,4(8X,F12.6))
  165 Z1=Z2
      GO TO 500
      END
CCOEVV        SUBROUTINE COEVV FOR BMD01T           AUGUST  9, 1965
      SUBROUTINE COEVV
          DIMENSION X(350),FMT(36),YSQ(80),Y(80),SYM(15),YP(15),
     1A(350),B(350),XX(350),YY(350)
      DIMENSION S(350)
      COMMON  X      , FMT    , YSQ    , Y      , SYM    , YP
      COMMON  A      , B      , XX     , YY     , AMPL   , PHASE
      COMMON  TG     , FILTER , KPROB  , NDATA  , FTRI   , NTRI
      COMMON  M      , FM     , NPOINT , POINT  , WO     , ANS
      COMMON  PI     , H
C     DIMENSION X(800),FMT(120),YSQ(500),Y(500),SYM(15),YP(15),
C    1A(800),B(800),XX(800),YY(800)
C     DIMENSION S(800)
C     COMMON X,FMT,YSQ,Y,SYM,YP,A,B,XX,YY,AMPL,PHASE,TG,FILTER,KPROB,
C    1NDATA,FTRI,NTRI,M,FM,NPOINT,POINT,WO,ANS,PI,H
      MMO=M-1
      ARG=FM*WO
      AA=2.0*FTRI*(.54+.46*COS(PI))*SIN(ARG)/FM
      BB=2.0*FTRI*(.54+.46*COS(PI))*COS(ARG)/FM
      KT=NTRI/2
      IF(NTRI-KT*2)20,10,20
   20 A(M)=-AA
      B(M)=BB
      GO TO 30
   10 A(M)=AA
      B(M)=-BB
   30 DO 105 J=1,MMO
      FJ=J
      ARGA=FJ*PI/FM
      ARGB=FJ*FTRI*PI/FM
  105 S(J)=4.0*(.54+.46*COS(ARGA))*SIN(ARGB)/(FM*SIN(ARGA))
      DO 120 J=1,MMO
      FJ=J
      ARG=FJ*WO
      A(J)=-S(J)*SIN(ARG)
  120 B(J)=S(J)*COS(ARG)
      RETURN
      END
CDINT    SUBROUTINE DINT FOR BMD01T                  MARCH 13, 1964
      SUBROUTINE DINT(YSQ,NN,H,ANS)
      DIMENSION YSQ(80),CONS(4)
C     DIMENSION YSQ(200),CONS(4)
C
      CONS(1)=0.34861111
      CONS(2)=1.24583333
      CONS(3)=0.87916667
      CONS(4)=1.02638889
      X=0.0
      J=NN
      DO 10 I=1,4
      X=X+(YSQ(I)+YSQ(J))*CONS(I)
      J=J-1
 10   CONTINUE
      DO 20 I=5,J
      X=X+YSQ(I)
 20   CONTINUE
      ANS=X*H
      RETURN
      END
CLANAI        SUBROUTINE LANAI FOR BMD01T           AUGUST  9, 1965
      SUBROUTINE LANAI (NN,BASE,KIT)
          DIMENSION X(350),FMT(36),YSQ(80),Y(80),SYM(15),YP(15),
     1A(350),B(350),XX(350),YY(350)
      COMMON  X      , FMT    , YSQ    , Y      , SYM    , YP
      COMMON  A      , B      , XX     , YY     , AMPL   , PHASE
      COMMON  TG     , FILTER , KPROB  , NDATA  , FTRI   , NTRI
      COMMON  M      , FM     , NPOINT , POINT  , WO     , ANS
      COMMON  PI     , H
C     DIMENSION X(800),FMT(120),YSQ(500),Y(500),SYM(15),YP(15),
C    1A(800),B(800),XX(800),YY(800)
C     COMMON X,FMT,YSQ,Y,SYM,YP,A,B,XX,YY,AMPL,PHASE,TG,FILTER,KPROB,
C    1NDATA,FTRI,NTRI,M,FM,NPOINT,POINT,WO,ANS,PI,H
      DO 100 I=1,NN
      FI=I
      OMEGA=FI*H+BASE
      CALL TABBY(OMEGA)
      L=I+KIT
      Y(L)=YY(1)
  100 YSQ(I)=YY(1)**2
      WRITE (6,800)
  800 FORMAT(1H0,21HVALUES OF F(X) SQUARE//)
      WRITE (6,1001)(YSQ(I),I=1,NN)
 1001 FORMAT(7F12.6)
      CALL DINT (YSQ,NN,H,ANS)
      RETURN
      END
CTABBY   SUBROUTINE TABBY FOR BMD01T               AUGUST  2, 1963
      SUBROUTINE TABBY(OMEGA)
          DIMENSION X(350),FMT(36),YSQ(80),Y(80),SYM(15),YP(15),
     1A(350),B(350),XX(350),YY(350)
      COMMON  X      , FMT    , YSQ    , Y      , SYM    , YP
      COMMON  A      , B      , XX     , YY     , AMPL   , PHASE
      COMMON  TG     , FILTER , KPROB  , NDATA  , FTRI   , NTRI
      COMMON  M      , FM     , NPOINT , POINT  , WO     , ANS
      COMMON  PI     , H
C     DIMENSION X(800),FMT(120),YSQ(500),Y(500),SYM(15),YP(15),
C    1A(800),B(800),XX(800),YY(800)
C     COMMON X,FMT,YSQ,Y,SYM,YP,A,B,XX,YY,AMPL,PHASE,TG,FILTER,KPROB,
C    1NDATA,FTRI,NTRI,M,FM,NPOINT,POINT,WO,ANS,PI,H
      SUM=0.0
      MMO=M-1
      DO 100 J=1,MMO
      FJ=J
      P=FJ*PI/FM
      Q=FJ*WO
      R=FTRI*FJ*PI/FM
      ARG=FJ*OMEGA
  100 SUM=SUM+(.54+.46*COS(P))*COS(Q)*SIN(R)*COS(ARG)/SIN(P)
      SUM=4.0*SUM/FM
      S=FM*WO
      T=FM*OMEGA
      NIT=NTRI/2
      IF(NTRI-NIT*2) 200,300,200
  200 YY(1)= (2.0*FTRI/FM)+SUM+2.0*FTRI*(.54+.46*COS(PI))*COS(S)*COS(T)/
     1FM
      GO TO 400
  300 YY(1)= (2.0*FTRI/FM)+SUM-2.0*FTRI*(.54+.46*COS(PI))*COS(S)*COS(T)/
     1FM
  400 RETURN
      END
CTPWD    SUBROUTINE TPWD FOR BMD01T         VERSION OF SEPT. 26, 1963
      SUBROUTINE TPWD(NT1,NT2)
      IF(NT1)40,10,12
 10   NT1=5
 12   IF(NT1-NT2)14,19,14
   14 IF(NT2.EQ.5) GO TO 19
      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 BMD01T             JUNE 2, 1964
      SUBROUTINE TRANS (MISTAK)
          DIMENSION X(350),FMT(36),YSQ(80),Y(80),SYM(15),YP(15),
     1A(350),B(350),XX(350),YY(350)
      DIMENSION IBIN (10),CON(10)
      COMMON  X      , FMT    , YSQ    , Y      , SYM    , YP
      COMMON  A      , B      , XX     , YY     , AMPL   , PHASE
      COMMON  TG     , FILTER , KPROB  , NDATA  , FTRI   , NTRI
      COMMON  M      , FM     , NPOINT , POINT  , WO     , ANS
      COMMON  PI     , H
C     DIMENSION X(800),FMT(120),YSQ(500),Y(500),SYM(15),YP(15),
C    1A(800),B(800),XX(800),YY(800)
C     DIMENSION IBIN (10),CON(10)
C     COMMON X,FMT,YSQ,Y,SYM,YP,A,B,XX,YY,AMPL,PHASE,TG,FILTER,KPROB,
C    1NDATA,FTRI,NTRI,M,FM,NPOINT,POINT,WO,ANS,PI,H
      EQUIVALENCE(ITG,SPECTG)
      ASN(Q000FL)=ATAN(Q000FL)/SQRT(1.0-Q000FL**2)
	DOUBLE PRECISION SPECTG,KODE
      DATA SPECTG/'SPECTG  '/
      READ (5,1002)KODE,NTRAN,(IBIN(I),CON(I),I=1,NTRAN)
 1002 FORMAT(A6,I1,8(I2,F6.0))
      IF(KODE-ITG)900,400,900
  400 DO 500 I=1,NTRAN
      IF(IBIN(I)-17) 605,600,605
  605 JESUS=IBIN(I)
      IF(JESUS-6)610,905,610
  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
   14 X(K)=SQRT(X(K))
      GO TO 150
   20 IF(X(K))200,22,23
   22 X(K)=1.0
      GO TO 150
   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
   40 X(K)=EXP(X(K))
      GO TO 150
   50 IF(X(K))200,150,53
   53 IF (X(K)-1.0) 54,55,200
  54  ARG =SQRT(X(K))
      X(K)=ASN(ARG)
      GO TO 150
   55 X(K)=PI/2.0
      GO TO 150
   60 IF(X(K))200,200,61
   61 X(K)=ALOG(X(K))
      GO TO 150
   70 IF(X(K))71,200,71
   71 X(K)=1.0/X(K)
      GO TO 150
   80 X(K)=X(K)+CC
      GO TO 150
   90 X(K)=X(K)*CC
      GO TO 150
  100 IF(X(K))200,200,101
  101 X(K)=X(K)**CC
      GO TO 150
  200 WRITE (6,1001)K,IBIN(I)
 1001 FORMAT(11H0DATA POINT I5,57H VIOLATES THE RESTRICTION FOR TRANSGEN
     1ERATION OF THE TYPEI3,52H. THE PROGRAM CONTINUES LEAVING THE VALUE
     2 UNCHANGED.)
  150 CONTINUE
  500 CONTINUE
  300 RETURN
  900 WRITE (6,1003)KODE
      WRITE (6,1006) KODE,NTRAN,(IBIN(I),CON(I),I=1,NTRAN)
 1006 FORMAT(1H0,A6,I1,8(I2,F6.0))
  901 WRITE (6,1004)
      MISTAK=17
      GO TO 300
 1003 FORMAT(58H0CONTROL CARD ERROR. PROGRAM EXPECTED A SPECTG CARD BUT 
     1A A6,16H CARD WAS FOUND.)
 1004 FORMAT(42H0PROGRAM WILL GO TO THE NEXT PROBLEM CARD.)
  905 WRITE(6,1005) KODE, NTRAN, (IBIN(I), CON(I), I=1,NTRAN)
 1005 FORMAT('0ILLEGAL TRANSGENERATION CODE SPECIFIED. CARD IS PRINTED'
     1 ' BELOW.' / (1H0, A6, I1, 8(I2,F6.0)))
      GO TO 901
      END
CPHAZE        SUBROUTINE PHAZE FOR BMD01T           AUGUST  9, 1965
      SUBROUTINE PHAZE(X,Y,Z)
      PI=3.14159265
      PI2=2.0*PI
      AB=ABS(Y/X)
      PHI=ATAN(AB)
      IF(X)11,12,13
   11 IF(Y)17,30,18
   17 Z=PI+PHI
      GO TO 50
   30 Z=PI
      GO TO 50
   18 Z=PI-PHI
      GO TO 50
   12 IF(Y)35,15,40
   35 Z=0.75*PI2
      GO TO 50
   15 Z=0.0
      GO TO 50
   40 Z=PI/2.0
      GO TO 50
   13 IF(Y)14,15,16
   14 Z=PI2-PHI
      GO TO 50
   16 Z=PHI
   50 RETURN
      END
C             SUBROUTINE PLOTR (IBM 360)              AUGUST 13, 1966
      SUBROUTINE PLOTR(X,ZMIN,ZMAX,Y,SYM,WMIN,WMAX,NC,NP)
C
      DIMENSION Y(15),CLAB(12),SYM(15),GF(10),FMT(12),XY(51,101)
	INTEGER XY, BLANKS
      INTEGER SYM,SYMB
      DATA NCC/2/
	DATA TC,TP,BLANKS/'.','+','     '/
      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)
  800 FORMAT(1H  14X,101A1)
 1001 FORMAT(15X,20(5H+....),1H+)
C
C        NCC ON THE INITIAL ENTRY TO PLOTR IS ASSUMED TO BE SOMETHING
C        OTHER THAN ZERO.
C
      IF(NCC) 50,48,50
   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
  201 IF(KL)220,220,231
  231 WRITE (6,1001)
      IF(KL)250,250,220
  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
  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
  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(JJJ),JJJ=1,9,2)
      IF(KL)227,227,14
  227 IF(JY-5)208,209,208
 208  WRITE (6,800)(TC,I=1,J),(TP,(TC ,I=1,4),K=1,19),TP,(TC,I=1,JYT)
      IF(KL) 250,250,225
  209 WRITE (6,1001)
      IF(KL) 250,250,225
  250 CONTINUE
      NCC=0
      IC=0
      IF(NP)80,11,11
   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
C      ENTRY PLOTS
   48 IF(NC)52,13,49
   49 IF(NP)80,10,10
   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
   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
   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
   52 KL=1
      GO TO 230
   14 NCC=1
   15 RETURN
      END
C             SUBROUTINE SCALE FOR PLOTR              JUNE 21, 1966
      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/
      TEST=1.0/(2.0**20)
   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)10,10,201
   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
	DIMENSION TEST(18)
       INTEGER XY,SYMB,BLANK,TEST
	DATA BLANK/'     '/
	DATA TEST/'2','3','4','5','6','7','8','9','A',
     1 '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
























