CID                     BMD06R
C        ASYMPTOTIC REGRESSION - MAIN PROGRAM          MAY  9, 1966
C        THIS IS A SIFTED VERSION OF BMD06R ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
C        IT WAS THEN CONVERTED TO 360 FORTRAN IV (H-LEVEL)
C
      COMMON  X1
      COMMON  NXY    , TXY    , NY     , D123   , TRANS  , LTRANX
      COMMON  CONST  , CONSX  , KVT    , SWITCH , Y1     , N1
      COMMON  NS     , N2     , NTL    , CONT   , IPASS  , IX
      COMMON  IY     , YV     , R0     , C      , C1     , TY0
      COMMON  TY1    , TY2    , R      , Y012   , RS     , AS
      COMMON  BS     , SUMSQS , AVE    , X2     , Y2     , Y3
      COMMON  I5     , I6
      DIMENSION X1(400),X2(400),CF(9),   Y2(1700),                NTL(
     110),CONT(10),Y3(1700),C(3,3),C1(3,3),Y012(3),RS(51),U(400),AS(50)
     2,BS(50),SUMSQS(51),AVE(400),PRE(400)
     3,N1(400),N2(400)
      EQUIVALENCE(X1,PRE,U)
      DOUBLE PRECISION A123,B123,TODE,PROB
      LOGICAL WRIT
      WRIT=.TRUE.
  502 FORMAT('1BMD06R - ASYMPTOTIC REGRESSION - REVISED ',
     1'OCTOBER 29, 1968'/
     2 40H HEALTH SCIENCES COMPUTING FACILITY,UCLA//)
C
      I5=5
	CALL USAGEB('BMD06R')
      I6=6
      DATA A123,B123/'FINISH','PROBLM'/
      DATA GOOG/3HYES/
      D123 =  GOOG
  100 WRITE(I6,502)
      READ(5,501)TODE,PROB,NXY,TXY,RES,SPARE,LTRANX,CONSX,XSORT,GPLOT,
     1STEP,XRED,MORE,OUT,RHO1,NTAPE,KVT
C
      IF(TODE.EQ.B123)GO TO 303
  301 IF(TODE.EQ.A123)WRIT=.FALSE.
C
C
C
  302 IF(WRIT)WRITE(I6,304)
      IF(I5.NE.5)REWIND I5
      STOP
  303 CALL TPWD(NTAPE,I5)
 332  CONTINUE
C
      WRITE (I6,701)
      WRITE (I6,503)
      WRITE (I6,504)PROB
      WRITE (I6,530)NXY,TXY
      WRITE (I6,532)XRED,XSORT
      WRITE (I6,535)RES
      WRITE (I6,537)LTRANX,OUT
      WRITE (I6,538)CONSX,KVT
      IF((NXY-3)*(NXY-400))305,302,302
  305 IF(KVT.GT.0.AND.KVT.LE.10)GO TO 3
      KVT=1
      WRITE(I6,4000)
 3    CONTINUE
      SUMSQS(1)=0.0
      SWITCH=0
      LTRANX=LTRANX+1
      CALL READ
      NCN=NS
C
C CHECK INPUT DATA
      DO 70 I=1,NXY
   70 N2(I)=0
      N2(1)=N1(1)
      DO 71 I=2,NXY
   71 N2(I)=N2(I-1)+N1(I)
      IF(OUT .NE. D123) GO TO 6
C PRINT INPUT DATA
    4 WRITE (I6,505)
      WRITE (I6,506)
      DO 5 I=1,NXY
      NE=N2(I)
      NS=NE-N1(I)+1
      DO 72 L=NS,NE
   72 WRITE (I6,507)I,X1(I),Y2(L)
    5 CONTINUE
    6 DO 7 I=1,NXY
      IF(N1(I))8,8,7
    7 CONTINUE
      GO TO 10
    8 WRITE (I6,508)
      IF(TRANS .NE. D123) GO TO 100
    9 DO 200 I=1,MORE
  200 READ (5,509)
      GO TO 100
C TRANSFORM X AND SORT X AND Y
   10 IF(LTRANX-1)    302,73,74
   74 CALL TRANX
      GO TO 75
   73 DO 76 I=1,NXY
   76 X2(I)=X1(I)
 75   IF(XSORT .NE. D123) GO TO 15
   11 CALL SORTX
      GO TO 20
   15 NF=0
      DO 16 I=1,NXY
      NE=N2(I)
      NS=NE-N1(I)+1
      DO 17 J=NS,NE
      NF=NF+1
   17 Y3(NF)=Y2(J)
   16 CONTINUE
 20   IF(XRED .NE. D123) GO TO 600
C SCALE X
  601 RNXY=NXY-1
      SMAL=X2(1)
      RANGE=X2(NXY)-SMAL
      DO 602 I=1,NXY
  602 X2(I)=((X2(I)-SMAL)/RANGE)*RNXY
 600  IF(MORE) 6001,26,27
 6001 WRITE(6,6002)
 6002 FORMAT('0ILLEGAL''TRANSGENERATION CARD COUNT'' SPECIFIED, MUST BE
     1 .GE. 1 AND .LE. 99.   PROGRAM TERMINATED.')
      GO TO 302
   26 MORE=1
C BEGIN REGRESSION LOOP
   27 DO 1000 IMM=1,MORE
      WRITE (I6,513)IMM
      NF=0
C SET UP Y FOR REGRESSION AND TRANSFORMATION
      DO 29 I=1,NXY
      NE=N2(I)
      NS=NE-N1(I)+1
      DO 28 L=NS,NE
      NF=NF+1
   28 Y2(NF)=Y3(L)
   29 CONTINUE
   21 CALL TRANSY(JQ)
      IF(JQ)25,301,251
  251 REWIND 1
      READ(1,501)TODE,PROB,NXY,TXY,RES,SPARE,LTRANX,CONSX,XSORT,GPLOT,
     1STEP,XRED,MORE,OUT,RHO1,NTAPE,KVT
      GO TO 303
C GET START VALUE FOR R
 25   CALL FINDR
      XN = NXY
      XB = 0.0
      DO 3993 I=1,NXY
 3993 XB = XB + X2(I)
      XB = XB/XN
      I = XB
      XB = I
      DO 3994 I=1,NXY
 3994 X2(I) = X2(I)-XB
      IF(RHO1)2511,2511,2510
 2510 R0 = RHO1
 2511 CONTINUE
      WRITE (I6,520)R0
      PRINT = 0.0
      IF(STEP .EQ. D123) GO TO 4877
 4777 CONTINUE
      WRITE (I6,522)
      PRINT = 1.0
 4877 CONTINUE
C  INITIALIZE LOOP
      ICOUNT=0
      NIT = 100
      RS(1)=R0
      TY0=0.0
      DO 445 I=1,NXY
      TY0=TY0+AVE(I)*(FLOAT(N1(I)))
  445 CONTINUE
      Y012(1)=TY0
      R=R0
      SUMSQ = 0.0
C BEGIN  ITERATION FOR ALPHA, BETA AND R.
 996  CONTINUE
      ICOUNT = ICOUNT + 1
      R0 = R
      CALL SETUP
      CALL INVERT(C1,3)
      A=C1(1,1)*TY0+C1(1,2)*TY1+C1(1,3)*TY2
      B=C1(1,2)*TY0+C1(2,2)*TY1+C1(2,3)*TY2
      R=R0+(C1(1,3)*TY0+C1(2,3)*TY1+C1(3,3)*TY2)/B
      R = ABS(R)
 40   CONTINUE
      SUMSQ=0.0
      DO 45 I=1,NXY
      YEST=A+B*(R**X2(I))
      DIF=AVE(I)-YEST
      SUMSQ=SUMSQ+DIF**2
   45 CONTINUE
      B = B*R**(-XB)
C PRINT ITERATION RESULTS IF DESIRED
      IF (PRINT )477,481,477
 477  CONTINUE
      WRITE (I6,523)ICOUNT,A, B, R, SUMSQ
      GO TO 47
 481  CONTINUE
      WRITE (I6,923)ICOUNT,A, B, R, SUMSQ
      WRITE (I6,518)
      DO 39 I=1,3
      WRITE (I6,516)
      K=I-1
      WRITE (I6,517)(C(I,J),J=1,3),K,Y012(I)
   39 WRITE (I6,516)
      WRITE (I6,519)
      DO 49 I=1,3
      WRITE (I6,516)
      K=I-1
      WRITE (I6,517)(C1(I,J),J=1,3),K,Y012(I)
   49 WRITE (I6,516)
 47   IF(ABS(R-R0)-ABS(R)*1.E-5)60,60,308
 308  IF (ICOUNT-NIT+1)99,99,309
 99   CONTINUE
      GO TO 996
C END OF ITERATION FOR ALPHA, BETA, AND R.
C PRINT RESULTS
 309  CONTINUE
      WRITE (I6,9309)NIT
 60   CONTINUE
      DO 6111 I=1,NXY
 6111 X2(I) = X2(I) +XB
      CALL SETUP
      CALL INVERT(C1,3)
      WRITE (I6,518)
      DO 62 I=1,3
      WRITE (I6,516)
      K=I-1
      WRITE (I6,517)(C(I,J),J=1,3),K,Y012(I)
   62 WRITE (I6,516)
      WRITE (I6,519)
      DO 63 I=1,3
      WRITE (I6,516)
      K=I-1
      WRITE (I6,517)(C1(I,J),J=1,3),K,Y012(I)
   63 WRITE (I6,516)
      WRITE (I6,524)
      WRITE (I6,525)
      N=0
      SE=0.0
      DO101 I=1,NXY
      NE=N2(I)
      NS=NE-N1(I)+1
      PRE(I)=A+B*(R**X2(I))
      N=N+N1(I)
      DO 101 J=NS,NE
      SE=SE+(Y2(J)-PRE(I))**2
  101 CONTINUE
      FN=N-3
      AS(1)= A
      BS(1)= B
      RS(1)= R
      SE=SE/FN
      SA=SQRT(C1(1,1)*SE)
      SB=SQRT(C1(2,2)*SE)
      SR=(SQRT(C1(3,3)*SE))/ABS(B)
      WRITE (I6,526)A,SA
      WRITE (I6,527)B,SB
      WRITE (I6,528)R, SR
      CALL ANOVAR
      IF(RES .NE. D123) GO TO 641
 640  CONTINUE
      WRITE (I6,529)
      WRITE (I6,630)
      DO 64 I=1,NXY
      NE=N2(I)
      NS=NE-N1(I)+1
      DO 65 J=NS,NE
      DIFF=Y2(J)-PRE(I)
   65 WRITE (I6,631)I,X2(I),Y2(J),PRE(I),DIFF
   64 CONTINUE
 641  CONTINUE
      IF(GPLOT .NE. D123) GO TO 66
   67 CALL PLOT  (ICOUNT)
   66 TRANS=D123
 1000 CONTINUE
      GO TO 100
C
  304 FORMAT(45H0CONTROL CARDS INCORRECTLY ORDERED OR PUNCHED)
 501  FORMAT(2A6,I3,A2,A3,A3,I2,F6.0,4A3,I2,A3,F6.0, 14X, 2I2)
  503 FORMAT(13H PROBLEM CARD//)
  504 FORMAT(13H PROBLEM CODE6X,A6)
  505 FORMAT(14H0ORIGINAL DATA//)
  506 FORMAT(5H  NO.6X,7HX VALUE6X,7HY VALUE//)
  507 FORMAT(1H I4,2F13.4)
  508 FORMAT(61H0NO Y VALUE FOR SOME X,PROGRAM WILL GO TO NEXT PROBLEM,I
     1F ANY)
  509 FORMAT(I2)
  513 FORMAT(9H0FIT NO. I3)
  515 FORMAT(15H0ITERATION NO. I2)
  516 FORMAT(2H *62X,1H*)
  517 FORMAT(2H *3F20.8,7H  *  Y(I1,2H)=F20.8)
  518 FORMAT(19H0INFORMATION MATRIX//)
  519 FORMAT(30H0INVERSE OF INFORMATION MATRIX//)
  520 FORMAT(/23H0INITIAL ESTIMATE OF R=F10.4//)
  521 FORMAT(/1H 25X,19HPARAMETER ESTIMATES)
  522 FORMAT(14H ITERATION NO.8X,1HA14X,1HB14X,1HR9X,20HSUM(E(Y)-MEAN(Y)
     1)**2//)
  523 FORMAT(1H I9,F19.6,2F15.6,1X,F16.6)
  524 FORMAT(/1H04X,5HFINAL6X,8HSTANDARD)
  525 FORMAT(1H 4X,9HESTIMATES2X,10HDEVIATIONS//)
  526 FORMAT(3H A=F11.4,2X,F10.4//)
  527 FORMAT(3H B=F11.4,2X,F10.4//)
  528 FORMAT(3H R=F11.4,2X,F10.4)
  529 FORMAT(20H1 TABLE OF RESIDUALS//)
  530 FORMAT(16H NO. OF X VALUES3X,I6,8X,14H INPUT PATTERN9X,A6)
 532  FORMAT(12H X RE-SCALED 10X,A6,8X,5H SORT17X,A6)
 535  FORMAT(     16H PRINT RESIDUALS 6X, A3)
  537 FORMAT(13H X TRANS CODE6X,I6,8X,12H OUTPUT DATA10X,A3)
  538 FORMAT(11H X CONSTANT8X,F11.5,3X,16H VARIABLE FORMAT3X,I6)
  630 FORMAT(69H0 NO.         X VALUE         Y VALUE     Y PREDICTED   
     1     RESIDUAL//)
  631 FORMAT(1H I4,4F16.4)
  701 FORMAT(20H REGRESSION EQUATION/1H015X,1HX/1H 4X,11HY = A + B.R//)
 923  FORMAT(11H0ITERATION=I5, 3X, 2HA=F12.6,3X,2HB=F12.6,3X,2HR=F12.6
     1,3X, 23HSUM(E(Y)-MEAN(Y))**2  =F12.6)
 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
     1IED, ASSUMED TO BE 1.)
 9309 FORMAT( 36H0CONVERGENCE CRITERION NOT MET AFTER  I5, 11H ITERATION
     1S  )
      END
      SUBROUTINE ANOVAR
C             SUBROUTINE ANOVAR FOR BMD06R             MAY  9, 1966
      COMMON  X1
      COMMON  NXY    , TXY    , NY     , D123   , TRANS  , LTRANX
      COMMON  CONST  , CONSX  , KVT    , SWITCH , Y1     , N1
      COMMON  NS     , N2     , NTL    , CONT   , IPASS  , IX
      COMMON  IY     , YV     , R0     , C      , C1     , TY0
      COMMON  TY1    , TY2    , R      , Y012   , RS     , AS
      COMMON  BS     , SUMSQS , AVE    , X2     , Y2     , Y3
      COMMON  I5     , I6
      DIMENSION X1(400),X2(400),         Y2(1700),                NTL(
     110),CONT(10),Y3(1700),C(3,3),C1(3,3),Y012(3),RS(51),U(400),AS(50)
     2,BS(50),SUMSQS(51),AVE(400),PRE(400)
     3,N1(400),N2(400)
      EQUIVALENCE(X1,PRE,U)
C
      ADD=0.0
      TOT=0.0
      TMEAN=0.0
      YHAT=0.0
C N2(NXY) CONTAINS THE TOTAL NUMBER OF DEPENDENT VARIABLE VALUES AND
C TY0 CONTAINS THE SUM OF ALL DEPENDENT VARIABLE VALUES.
C SUM, THEREFORE, CONTAINS THE GRAND MEAN.
      DIV=N2(NXY)
      SUM=TY0/DIV
      DO 10 I=1,NXY
      NE=N2(I)
      NS=NE-N1(I)+1
      YHAT=YHAT+((AVE(I)-PRE(I))**2)*FLOAT(N1(I))
      DO 7 J=NS,NE
      ADD=ADD+(Y2(J)-SUM)**2
      TOT=TOT+(Y2(J)-AVE(I))**2
    7 TMEAN=TMEAN+(Y2(J)-PRE(I))**2
   10 CONTINUE
      WRITE (I6,100)
      WRITE (I6,101)
      WRITE (I6,106)
      WRITE (I6,102)ADD
      WRITE (I6,103)TOT
      WRITE (I6,104)TMEAN
      WRITE (I6,105)YHAT
      WRITE (I6,106)
  100 FORMAT(1H020X,20HANALYSIS OF VARIANCE//)
  101 FORMAT(11H DEVIATIONS20X,14HSUM OF SQUARES)
  102 FORMAT(11H FROM MEAN 20X,F13.4)
  103 FORMAT(16H FROM X(I) MEANS15X,F13.4//)
  104 FORMAT(11H FROM CURVE20X,F13.4)
  105 FORMAT(25H OF X(I) MEANS FROM CURVE6X,F13.4)
  106 FORMAT(1H 44(1H*))
      RETURN
      END
      SUBROUTINE FINDR
C           SUBROUTINE FINDR FOR BMD06R                MAY  9, 1966
      COMMON  X1
      COMMON  NXY    , TXY    , NY     , D123   , TRANS  , LTRANX
      COMMON  CONST  , CONSX  , KVT    , SWITCH , Y1     , N1
      COMMON  NS     , N2     , NTL    , CONT   , IPASS  , IX
      COMMON  IY     , YV     , R0     , C      , C1     , TY0
      COMMON  TY1    , TY2    , R      , Y012   , RS     , AS
      COMMON  BS     , SUMSQS , AVE    , X2     , Y2     , Y3
      COMMON  I5     , I6
      DIMENSION X1(400),X2(400),         Y2(1700),                NTL(
     110),CONT(10),Y3(1700),C(3,3),C1(3,3),Y012(3),RS(51),U(400),AS(50)
     2,BS(50),SUMSQS(51),AVE(400),PRE(400)
     3,N1(400),N2(400)
      DIMENSIONCF(9)
      EQUIVALENCE(X1,PRE,U)
C
      DO 1 I=1,NXY
      SUM=0.0
      NE=N2(I)
      NS=NE-N1(I)+1
      DO 2 J=NS,NE
 2    SUM = SUM + Y2(J)
      AVE(I) = SUM/FLOAT(N1(I))
 1    CONTINUE
      CALL REGRES( X2, AVE,NXY, CF)
      IF(ABS(CF(2))-.01)50,50,55
 50   R0 = .9999
      GO TO 100
 55   CONTINUE
      IF(CF(5))60,65,60
 65   CONTINUE
      GO TO 50
 60   CONTINUE
      R0 = (CF(4)*(X2(2)-X2(NXY))+CF(5)*(X2(2)**2 -X2(NXY)**2))/
     1(CF(4)*(X2(1)-X2(NXY-1))+CF(5)*(X2(1)**2-X2(NXY-1)**2))
      IF(R0)50, 50,100
 100  RETURN
      END
      SUBROUTINE INVERT(A,N)
C            SUBROUTINE INVERT FOR BMD06R              MAY  9, 1966
      COMMON  X1
      COMMON  NXY    , TXY    , NY     , D123   , TRANS  , LTRANX
      COMMON  CONST  , CONSX  , KVT    , SWITCH , Y1     , N1
      COMMON  NS     , N2     , NTL    , CONT   , IPASS  , IX
      COMMON  IY     , YV     , R0     , C      , C1     , TY0
      COMMON  TY1    , TY2    , R      , Y012   , RS     , AS
      COMMON  BS     , SUMSQS , AVE    , X2     , Y2     , Y3
      COMMON  I5     , I6
      DIMENSION X1(400),X2(400),         Y2(1700),                NTL(
     110),CONT(10),Y3(1700),C(3,3),C1(3,3),Y012(3),RS(51),U(400),AS(50)
     2,BS(50),SUMSQS(51),AVE(400),PRE(400)
     3,N1(400),N2(400)
      DIMENSION L(3),M(3),A(3,3)
      EQUIVALENCE(X1,PRE,U)
C
      D=1.0
      DO80 K=1,N
      L(K)=K
      M(K)=K
      BIGA=A(K,K)
      DO20 I=K,N
      DO20 J=K,N
      IF(ABS(BIGA)-ABS(A(I,J))) 10,20,20
   10 BIGA=A(I,J)
      L(K)=I
      M(K)=J
   20 CONTINUE
C     INTERCHANGE ROWS
      J=L(K)
      IF(L(K)-K) 35,35,25
   25 DO30 I=1,N
      HOLD=-A(K,I)
      A(K,I)=A(J,I)
   30 A(J,I)=HOLD
C     INTERCHANGE COLUMNS
   35 I=M(K)
      IF(M(K)-K) 45,45,37
   37 DO40 J=1,N
      HOLD=-A(J,K)
      A(J,K)=A(J,I)
   40 A(J,I)=HOLD
C     DIVIDE COLUMN BY MINUS PIVOT
   45 DO55 I=1,N
   46 IF(I-K)50,55,50
   50 A(I,K)=A(I,K)/(-A(K,K))
   55 CONTINUE
C     REDUCE MATRIX
      DO65 I=1,N
      DO65 J=1,N
   56 IF(I-K) 57,65,57
   57 IF(J-K) 60,65,60
   60 A(I,J)=A(I,K)*A(K,J)+A(I,J)
   65 CONTINUE
C     DIVIDE ROW BY PIVOT
      DO75 J=1,N
   68 IF(J-K)70,75,70
   70 A(K,J)=A(K,J)/A(K,K)
   75 CONTINUE
C     CONTINUED PRODUCT OF PIVOTS
      D=D*A(K,K)
C     REPLACE PIVOT BY RECIPROCAL
      A(K,K)=1.0/A(K,K)
   80 CONTINUE
C     FINAL ROW AND COLUMN INTERCHANGE
      K=N
  100 K=(K-1)
      IF(K) 150,150,103
  103 I=L(K)
      IF(I-K) 120,120,105
  105 DO110 J=1,N
      HOLD=A(J,K)
      A(J,K)=-A(J,I)
  110 A(J,I)=HOLD
  120 J=M(K)
      IF(J-K) 100,100,125
  125 DO130 I=1,N
      HOLD=A(K,I)
      A(K,I)=-A(J,I)
  130 A(J,I)=HOLD
      GO TO 100
 150  CONTINUE
      IF(D) 160,999,160
 999  WRITE (6,995)
 995  FORMAT(49H ERROR IN ATTEMPTING TO INVERT A SINGULAR MATRIX. )
 160  RETURN
      END
      SUBROUTINE PLOT (ICOUNT)
C             SUBROUTINE  PLOT FOR BMD06R              MAY  9, 1966
      COMMON  X1
      COMMON  NXY    , TXY    , NY     , D123   , TRANS  , LTRANX
      COMMON  CONST  , CONSX  , KVT    , SWITCH , Y1     , N1
      COMMON  NS     , N2     , NTL    , CONT   , IPASS  , IX
      COMMON  IY     , YV     , R0     , C      , C1     , TY0
      COMMON  TY1    , TY2    , R      , Y012   , RS     , AS
      COMMON  BS     , SUMSQS , AVE    , X2     , Y2     , Y3
      COMMON  I5     , I6
      DIMENSION X1(400),X2(400),         Y2(1700),                NTL(
     110),CONT(10),Y3(1700),C(3,3),C1(3,3),Y012(3),RS(51),U(400),AS(50)
     2,BS(50),SUMSQS(51),AVE(400),PRE(400)
     3,N1(400),N2(400)
       DIMENSION FMT(120)
      DIMENSION SYM(15),YYY(15)
      EQUIVALENCE(X1,PRE,U)
C
      A=AS(ICOUNT)
      B=BS(ICOUNT)
      R=RS(ICOUNT+1)
      SS=10.0**10
      FF=-10.0**10
      XS=SS
      XL=FF
      DO 300 I=1,NXY
      XS=AMIN1(XS,X2(I))
      XL=AMAX1(XL,X2(I))
      SS=AMIN1(SS,AVE(I),PRE(I))
  300 FF=AMAX1(FF,AVE(I),PRE(I))
      WRITE (I6,100)
      WRITE (I6,201)
C     DIMENSION SYM(15),YYY(15)
      DATA SYM(1),SYM(2)/'A','P'/
  302 XX=X2(1)
      YYY(2)=PRE(1)
      YYY(1)=AVE(1)
      NC=2
      CALL PLOTR(XX,XS,XL,YYY,SYM,SS,FF,NC,+1)
      DO 50 I=2,NXY
      XX=X2(I)
      YYY(2)=PRE(I)
      YYY(1)=AVE(I)
 35   CALL PLOTR(XX,XS,XL,YYY,SYM,SS,FF,NC,+1)
   50 CONTINUE
      NC = 0
      CALL PLOTR(XX,XS,XL,YYY,SYM,SS,FF,NC,+1)
 100  FORMAT(12H1GRAPH CODES,9X,13H A=AVERAGED Y,5X,14H P=PREDICTED Y,
     15X,7H B=BOTH)
  201 FORMAT(48H0X AND Y VALUES ARE PLOTTED IN TRANSFORMED UNITS)
  403 RETURN
      END
      SUBROUTINE READ
C             SUBROUTINE READ FOR BMD06R               MAY  9, 1966
      COMMON  X1
      COMMON  NXY    , TXY    , NY     , D123   , TRANS  , LTRANX
      COMMON  CONST  , CONSX  , KVT    , SWITCH , Y1     , N1
      COMMON  NS     , N2     , NTL    , CONT   , IPASS  , IX
      COMMON  IY     , YV     , R0     , C      , C1     , TY0
      COMMON  TY1    , TY2    , R      , Y012   , RS     , AS
      COMMON  BS     , SUMSQS , AVE    , X2     , Y2     , Y3
      COMMON  I5     , I6
      DIMENSION X1(400),X2(400),FMT(180),Y2(1700),Z(400),        NTL(
     110),CONT(10),Y3(1700),C(3,3),C1(3,3),Y012(3),RS(51),U(400),AS(50)
     2,BS(50),SUMSQS(51),AVE(400),PRE(400)
     3,N1(400),N2(400)
      EQUIVALENCE (X1,PRE,U),(IA1,A1),(Z,X2)
      DOUBLE PRECISION A1,IA1,ISAM
      DATA A1,G123/'SAMSIZ','XY'/
C
      NS=0
      NCARDS=(NXY+21)/22
      L2=0
      DO 50 J=1,NCARDS
      L1=L2+1
      L2=L2+22
      IF(NXY-L2)35,36,36
 35   L2=NXY
 36   READ (5,103)ISAM,(N1(I),I=L1,L2)
      IF(ISAM.EQ.IA1)GO TO 50
 40   WRITE (6,4000)A1
      NXY=-22
 50   CONTINUE
      IF(NXY)55,55,56
   55 STOP
 56   NTOT=0
      DO 70 I=1,NXY
      IF(N1(I))80,80,75
 80   N1(I)= 1
 75   NTOT=NTOT+N1(I)
 70   CONTINUE
      IF(NTOT-1700)4,4,300
 4    KVT=KVT*18
      READ (5,101)(FMT(I),I=1,KVT)
      WRITE(6,41) (FMT(I), I=1,KVT)
 41   FORMAT(' THE VARIABLE FORMAT IS  ',20A4,/(1H0,24X,20A4))
      JUMP=1
 5    IF(TXY .EQ. G123) GO TO 7
    6 JUMP=2
    7 DO 25 I=1,NXY
      NY=N1(I)
      GO TO (9,11), JUMP
    9 READ (I5,FMT)D,(Z(J),J=1,NY)
      X1(I)=D
      DO 19 J=1,NY
   16 NS=NS+1
      Y2(NS)=Z(J)
   19 CONTINUE
      GO TO 25
   11 READ (I5,FMT)(Z(J),J=1,NY),D
      X1(I)=D
      DO 24 J=1,NY
   26 NS=NS+1
      Y2(NS)=Z(J)
   24 CONTINUE
   25 CONTINUE
 101  FORMAT(18A4)
  102 FORMAT(12F6.0)
 103  FORMAT(A6,22I3)
 4000 FORMAT(1H0,A6,33H CARD MISPUNCHED OR OUT OF ORDER.)
 4010 FORMAT(84H0TOTAL NUMBER OF VALUES FOR DEPENDENT VARIABLE EXCEEDS 5
     1000. PROGRAM CANNOT PROCEED.)
      RETURN
 300  WRITE (6,4010)
      GO TO 55
      END
      SUBROUTINE REGRES (X,  F,  N,  C)
C             SUBROUTINE REGRES FOR BMD06R             MAY  9, 1966
      DIMENSIONA(9),B(3),C(9),F(1),X(1)
C
      DO 1 I=1,9
 1    A(I)=0.0
      B(1) = 0.0
      B(2) = 0.0
      B(3) = 0.0
      DO 100 I=1,N
      X2 = X(I)**2
      A(2) = A(2)+X(I)
      A(3) = A(3) +X2
      A(6) = A(6) + X2*X(I)
      A(9) = A(9) + X2**2
      B(1) = B(1) + F(I)
      B(2) = B(2) +F(I)*X(I)
      B(3) = B(3)+ F(I)*X2
 100  CONTINUE
      A(1) = N
      A(4) = A(2)
      A(5) = A(3)
      A(7) = A(3)
      A(8) = A(6)
      D = A(1)*A(5)-A(2)**2
      C(1) = (A(5)*B(1)-A(2)*B(2))/D
      C(2) = (-A(2)*B(1)+A(1)*B(2))/D
      CALL INVERT(A,3)
      C(3) = A(1)*B(1)+A(2)*B(2) +A(3)*B(3)
      C(4) = A(4) *B(1) + A(5)*B(2) + A(6)*B(3)
      C(5) = A(7) *B(1) + A(8)*B(2) + A(9)*B(3)
      RETURN
C
      END
      SUBROUTINE SETUP
C        SUBROUTINE SETUP FOR BMD06R                   MAY  9, 1966
      COMMON  X1
      COMMON  NXY    , TXY    , NY     , D123   , TRANS  , LTRANX
      COMMON  CONST  , CONSX  , KVT    , SWITCH , Y1     , N1
      COMMON  NS     , N2     , NTL    , CONT   , IPASS  , IX
      COMMON  IY     , YV     , R0     , C      , C1     , TY0
      COMMON  TY1    , TY2    , R      , Y012   , RS     , AS
      COMMON  BS     , SUMSQS , AVE    , X2     , Y2     , Y3
      COMMON  I5     , I6
      DIMENSION X1(400),X2(400),         Y2(1700),                NTL(
     110),CONT(10),Y3(1700),C(3,3),C1(3,3),Y012(3),RS(51),U(400),AS(50)
     2,BS(50),SUMSQS(51),AVE(400),PRE(400)
     3,N1(400),N2(400)
      EQUIVALENCE(X1,PRE,U)
C
      V=ALOG(R)
      DO 10 I=1,3
      DO 10 J=1,3
   10 C(I,J)=0.0
      TY1=0.0
      TY2=0.0
      DO 20 I=1,NXY
      FN=N1(I)
      C(1,1)=C(1,1)+FN
      U(I)=EXP(X2(I)*V)
      U3=X2(I)*U(I)/R
      U4=U(I)**2
      U5=X2(I)*U4/R
      U6=X2(I)*U5/R
      C(1,2)=C(1,2)+(U(I)*FN)
      C(1,3)=C(1,3)+(U3*FN)
      C(2,2)=C(2,2)+(U4*FN)
      C(2,3)=C(2,3)+(U5*FN)
      C(3,3)=C(3,3)+(U6*FN)
      TY1=TY1+U(I)*AVE(I)*FN
   15 TY2=TY2+U3*AVE(I)*FN
   20 CONTINUE
      C(3,1)=C(1,3)
      C(3,2)=C(2,3)
      C(2,1)=C(1,2)
      Y012(2)=TY1
      Y012(3)=TY2
      DO 30 I=1,3
      DO 30 J=1,3
   30 C1(I,J)=C(I,J)
      RETURN
      END
      SUBROUTINE SORTX
C           SUBROUTINE SORTX FOR BMD06R                MAY  9, 1966
      COMMON  X1
      COMMON  NXY    , TXY    , NY     , D123   , TRANS  , LTRANX
      COMMON  CONST  , CONSX  , KVT    , SWITCH , Y1     , N1
      COMMON  NS     , N2     , NTL    , CONT   , IPASS  , IX
      COMMON  IY     , YV     , R0     , C      , C1     , TY0
      COMMON  TY1    , TY2    , R      , Y012   , RS     , AS
      COMMON  BS     , SUMSQS , AVE    , X2     , Y2     , Y3
      COMMON  I5     , I6
      DIMENSION  N3(400),N4(400)       ,N1(400),N2(400)
      DIMENSION X1(400),X2(400),         Y2(1700),                NTL(
     110),CONT(10),Y3(1700),C(3,3),C1(3,3),Y012(3),RS(51),U(400),AS(50)
     2,BS(50),SUMSQS(51),AVE(400),PRE(400)
      EQUIVALENCE(X1,PRE,U),(AVE,N3)
C
      NF=0
      DO 60 I=1,NXY
   33 N3(I)=I
   60 N4(I)=N1(I)
      DO 10 I=1,NXY
      DO 5 J=I,NXY
      IF(X2(I)-X2(J)) 5,5,4
    4 XD=X2(I)
      X2(I)=X2(J)
      X2(J)=XD
      K1=N3(I)
      N3(I)=N3(J)
      N3(J)=K1
      K2=N1(I)
      N1(I)=N1(J)
      N1(J)=K2
    5 CONTINUE
   10 CONTINUE
      DO 6 I=1,NXY
      K1=N3(I)
      NE=N2(K1)
      NS=NE-N4(K1)+1
      DO 3 K=NS,NE
      NF=NF+1
    3 Y3(NF)=Y2(K)
    6 CONTINUE
      DO 20 I=1,NXY
   20 N2(I)=0
      N2(1)=N1(1)
      DO 30 I=2,NXY
   30 N2(I)=N2(I-1)+N1(I)
      RETURN
      END
      SUBROUTINE TPWD(NT1,NT2)
C        SUBROUTINE TPWD FOR BMD06R                    MAY  9, 1966
C
C
      IF(NT1)40,10,12
 10   NT1=5
 12   IF(NT1-NT2)14,19,14
   14 IF(NT2.EQ.5)GO TO 18
 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
C
C
C
C
   40 WRITE(6,49)
C
C
      STOP
 49   FORMAT(25H ERROR ON TAPE ASSIGNMENT)
      END
      SUBROUTINE TRANX
C             SUBROUTINE TRANX FOR BMD06R              MAY  9, 1966
      COMMON  X1
      COMMON  NXY    , TXY    , NY     , D123   , TRANS  , LTRANX
      COMMON  CONST  , CONSX  , KVT    , SWITCH , Y1     , N1
      COMMON  NS     , N2     , NTL    , CONT   , IPASS  , IX
      COMMON  IY     , YV     , R0     , C      , C1     , TY0
      COMMON  TY1    , TY2    , R      , Y012   , RS     , AS
      COMMON  BS     , SUMSQS , AVE    , X2     , Y2     , Y3
      COMMON  I5     , I6
      DIMENSION X1(400),X2(400),         Y2(1700),                NTL(
     110),CONT(10),Y3(1700),C(3,3),C1(3,3),Y012(3),RS(51),U(400),AS(50)
     2,BS(50),SUMSQS(51),AVE(400),PRE(400)
     3,N1(400),N2(400)
      EQUIVALENCE(X1,PRE,U)
      ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
C
      LTRANX=LTRANX-1
      KTRANX=LTRANX
      IF(LTRANX*(LTRANX-11))5,150,175
    5 DO 100 I=1,NXY
      D=X1(I)
      GO TO (10,20,30,40,50,60,70,80,90,200,170),LTRANX
 10   IF(D)160,11,12
 11   X2(I)=0.0
      GO TO 100
 12   X2(I)=SQRT(D)
      GO TO 100
 20   IF(D)160,21,22
 21   X2(I)=1.0
      GO TO 100
 22   X2(I)=SQRT(D)+SQRT(D+1.0)
      GO TO 100
 30   IF(D)160,160,31
 31   X2(I) = ALOG10(D)
      GO TO 100
   40 X2(I)=EXP(D)
      GO TO 100
 50   IF(D)160,11,51
 51   IF(D-1.0)52,53,160
 52   A=SQRT(D)
      X2(I)=ASN(A)
      GO TO 100
 53   X2(I)=1.57079633
      GO TO 100
   60 FN=NXY
      A=D/(FN+1.0)
      B=A+1.0/(FN+1.0)
      IF(A)160,61,62
 61   IF(B)160,11,64
 62   IF(B)160,65,66
 64   X2(I)=ASN(SQRT(B))
      GO TO 100
 65   X2(I)=ASN(SQRT(A))
      GO TO 100
 66   X2(I)=ASN(SQRT(A))+ASN(SQRT(B))
      GO TO 100
 70   IF(D)71,160,71
 71   X2(I)=1.0/D
      GO TO 100
   80 X2(I)=D+CONSX
      GO TO 100
   90 X2(I)=D*CONSX
      GO TO 100
 200  IF(D)210,11,210
 210  X2(I)=D**CONSX
      GO TO 100
 170  IF(D)155,155,171
 171  X2(I) = ALOG(D)
      GO TO 100
 150  WRITE (6,4000)
 4000 FORMAT(108H0ILLEGAL TRANSGENERATION CODE FOR THE INDEPENDENT VARIA
     1BLE. PROGRAM WILL PROCEED LEAVING VARIABLE UNCHANGED.)
      X2(I)=D
      GO TO 100
 155  KTRANX=17
 160  WRITE (6,4010)I,KTRANX
      X2(I)=D
      GO TO 100
 4010 FORMAT(4H0THEI4,94HTH. VALUE OF THE INDEPENDENT VARIABLE VIOLATES   
     1THE RESTRICTION FOR TRANSGENERATION OF THE TYPEI3,1H./52H THE PROG
     2RAM CONTINUES LEAVING THIS VALUE UNCHANGED.)
 175  IF(LTRANX-17)150,176,150
 176  LTRANX=11
      GO TO 5
  100 CONTINUE
      RETURN
      END
      SUBROUTINE TRANSY(JQ)
C             SUBROUTINE TRANSY FOR BMD06R             MAY  9, 1966
      COMMON  X1
      COMMON  NXY    , TXY    , NY     , D123   , TRANS  , LTRANX
      COMMON  CONST  , CONSX  , KVT    , SWITCH , Y1     , N1
      COMMON  NS     , N2     , NTL    , CONT   , IPASS  , IX
      COMMON  IY     , YV     , R0     , C      , C1     , TY0
      COMMON  TY1    , TY2    , R      , Y012   , RS     , AS
      COMMON  BS     , SUMSQS , AVE    , X2     , Y2     , Y3
      COMMON  I5     , I6
      DIMENSION X1(400),X2(400),         Y2(1700),                NTL(
     110),CONT(10),Y3(1700),C(3,3),C1(3,3),Y012(3),RS(51),U(400),AS(50)
     2,BS(50),SUMSQS(51),AVE(400),PRE(400)
     3,N1(400),N2(400),MT(19)
      EQUIVALENCE(X1,PRE,U)
      DOUBLE PRECISION TODE,B123,C123
      DATA B123,C123/'PROBLM  ','SPECTG  '/
      ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
      JQ=-1
      READ(5,1100)TODE,MT
      REWIND 1
      WRITE(1,1100)TODE,MT
 1100 FORMAT(A6,A2,18A4)
      IF(TODE.EQ.B123)GO TO 2
      IF(TODE.NE.C123)GO TO 1
      REWIND 1
      READ(1,1101)TODE,NT,(NTL(I),CONT(I),I=1,NT)
  300 WRITE (I6,1102)
      WRITE (I6,1103)
      DO 114 I=1,NT
  114 WRITE (I6,1104)NTL(I),CONT(I),I
    9 DO 99 I=1,NT
      JUMP=NTL(I)+1
      FMULY=CONT(I)
      IF(JUMP*(JUMP-19))200,203,203
 200  IF(11-JUMP)204,201,201
 204  IF(JUMP-18)203,202,203
 202  JUMP=12
  201 DO 98 J=1,NXY
      NE=N2(J)
      NS=NE-N1(J)+1
      GO TO(99,10,20,30,40,50,60,70,80,90,100,110),JUMP
   10 DO 11 K=NS,NE
      IF(Y2(K))88,13,14
   13 Y2(K)=0.0
      GO TO 11
   14 Y2(K)=SQRT(Y2(K))
   11 CONTINUE
      GO TO 98
   20 DO 21 K=NS,NE
      IF(Y2(K))88,22,23
   22 Y2(K)=1.0
      GO TO 21
   23 Y2(K)=SQRT(Y2(K))+SQRT(Y2(K)+1.0)
   21 CONTINUE
      GO TO 98
   30 DO 31 K=NS,NE
      IF(Y2(K))88,88,31
 31   Y2(K) = ALOG10(Y2(K))
      GO TO 98
   40 DO 41 K=NS,NE
   41 Y2(K)=EXP(Y2(K))
      GO TO 98
   50 DO 51 K=NS,NE
      IF(Y2(K))88,52,53
   52 Y2(K)=0.0
      GO TO 51
   53 IF(Y2(K)-1.0)54,55,88
   54 Y2(K)=ASN(SQRT(Y2(K)))
      GO TO 51
   55 Y2(K)=3.14159265/2.0
   51 CONTINUE
      GO TO 98
   60 DIV=N1(J)
      DO 61 K=NS,NE
      A=Y2(K)/(DIV+1.0)
      B=A+1.0/(DIV+1.0)
      IF(A)88,62,63
   62 IF(B)88,64,65
   64 Y2(K)=0.0
      GO TO 61
   65 Y2(K)=ASN(SQRT(B))
      GO TO 61
   63 IF(B)88,66,67
   66 Y2(K)=ASN(SQRT(A))
      GO TO 61
   67 Y2(K)=ASN(SQRT(A))+ASN(SQRT(B))
   61 CONTINUE
      GO TO 98
   70 DO 71 K=NS,NE
      IF(Y2(K)) 71,88,71
   71 Y2(K)=1.0/Y2(K)
      GO TO 98
   80 DO 81 K=NS,NE
   81 Y2(K)=Y2(K)+FMULY
      GO TO 98
   90 DO 91 K=NS,NE
   91 Y2(K)=Y2(K)*FMULY
      GO TO 98
  100 DO 101 K=NS,NE
      IF(Y2(K))102,103,102
  103 Y2(K)=0.0
      GO TO 101
  102 Y2(K)=Y2(K)**FMULY
  101 CONTINUE
      GO TO 98
  110 DO 111 K=NS,NE
      IF(Y2(K))88,88,111
  111 Y2(K)=ALOG(Y2(K))
   98 CONTINUE
   99 CONTINUE
      GO TO 1000
 203  WRITE (6,4000)NTL(I)
      GO TO 99
 4000 FORMAT(29H0ILLEGAL TRANSGENERATION CODEI3,93H SPECIFIED FOR THE DE
     1PENDENT VARIABLE. THE PROGRAM WILL PROCEED WITHOUT THIS TRANSGENER
     2ATION.)
 88   N=K-NS+1
      WRITE (6,4010)N,J,NTL(I)
      GO TO 98
 4010 FORMAT(4H0THEI4,56HTH. VALUE OF THE DEPENDENT VARIABLE CORRESPONDI
     1NG TO THEI4,37HTH. VALUE OF THE INDEPENDENT VARIABLE/57H VIOLATES   
     2THE RESTRICTION FOR TRANSGENERATION OF THE TYPEI3,1H.52H THE PROGR
     3AM CONTINUES LEAVING THIS VALUE UNCHANGED.)
 1101 FORMAT(A6,I1,8(I2,F6.0))
 1102 FORMAT(20H0TRANSFORMATION CARD)
 1103 FORMAT(27H0 CODE   CONSTANT  PASS NO.)
 1104 FORMAT(1H0I5,F13.4,I5)
    2 JQ=JQ+1
    1 JQ=JQ+1
 1000 RETURN
      END
      SUBROUTINE PLOTR(X,ZMIN,ZMAX,Y,SYM,WMIN,WMAX,NC,NP)
C
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
      DIMENSION Y(15),CLAB(12),GF(10),FMT(12),XY(51,101),SYM(15)
      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
   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
      SUBROUTINE SCALE(YMIN,YMAX,YINT,JY,TYMIN,TYMAX,YIJ)
C             SUBROUTINE SCALE FOR PLOTR              JUNE 21, 1966
C
C     SUBROUTINE 'SCALE' CALCULATES THE SCALING FOR 'PLOTR'
C
      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**17)
   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
      INTEGER XY,SYMB,BLANK,B
      DATA BLANK/'  '/,B/'B '/
C
      IF(XY.EQ.BLANK)GO TO 50
C
      IF(XY.NE.SYMB)XY=B
C
C
      GO TO 100
C        IF BLANK, PUT IN SYMBOL
   50 XY=SYMB
  100 RETURN
      END
C
