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