Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/bmd/bmd06r.for
There is 1 other file named bmd06r.for in the archive. Click here to see a list.
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