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