Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/bmd/bmdx92.for
There is 1 other file named bmdx92.for in the archive. Click here to see a list.
DIMENSION Q(5000)
DIMENSION NAME(32),TSV(19)
REAL NAME
DATA NAME /' 1',' 2',' 3',' 4',' 5',' 6',' 7',' 8',' 9',
C '10','11','12','13','14','15','16','17','18','19','20','21','22',
C '23','24','25','26','27','28','29','30','31','32'/
DOUBLE PRECISION VIN,VOUT
DATA VIN/' INPUT'/,VOUT/'OUTPUT'/
DATA LOW/'LOW '/,ORM/'ORM '/
DOUBLE PRECISION NAMCHN
DOUBLE PRECISION PROB,FIN,TIT,CHAN,PAR,XVA,YVA,PRE,LABL,LIBL,
1CODE,LBTAPE,LSV
DATA PROB/'PROBLM '/,FIN/'FINISH '/,TIT/'TITLES '/,CHAN/ '
1NAMCHN '/,PAR/'PARAM '/,XVA/'XVALUE '/,YVA/'YVALUE '/,PRE/ '
2PREFLT '/
INTEGER BCODE,HCODE,BPULSE,HPULSE,EXPNO,PSTART,PEND
COMMON/BRI/FILE(2),LONG,BCODE,HCODE,LENGTH,BPULSE,HPULSE,EXPNO,
1 PSTART,PEND,TSTART,TEND,IDAY,IHR,IMIN,ISEC,FS,ISKIP
COMMON/BRIID/LBTAPE,LBUNIT,IDDPL(29),NID,TITLER(68),ID1BRI,ID2BRI
COMMON/XDATA/ICODE,COM1(770)
COMMON/TABLE/COM2(12)
COMMON/PTABLE/COM3(402)
COMMON/CTABLE/COM4(302)
COMMON/SCODE/COM5(7)
COMMON/SCOUNT/COM6(5)
COMMON/LFILE/COM7(126)
COMMON/GETIME/COM8(6)
COMMON/PSAVEC/COM9(17)
COMMON/PREF/X(11),Y(11),NPOINT
COMMON /LIST/ NSC, LEN4, MSCANS, NV, NCHANS, NB, SR, FILTER, ID,
1 L, PLTFLT, SPEC, PLTWND, IT, DPLTAP, IBEGIN, NF, NWORDS, NWTS,
2 SAMPRT, LOG2NS, PLTPOW, PLTLOG, NS, NSCANS,NFF,LFT,LCT
COMMON /FRITOS/ NTOT,JSCANS,ISTART,IWTS,ITEMP,IRR,LOGIRR,IFT2,
X IROWSP,IAMPL,ILOGA,ISPEC,IBRI,NWO,ISIZE2
DOUBLE PRECISION CADE,STIM,BLANK,NULL
DATA CADE/'CODE '/,STIM/'STIM '/,BLANK/' '/,NULL/'N
1ULL '/
DATA YES, BNO, ONO / 'YES ',' NO ','NO '/
COMMON/CHANEL/NAMCHN(32),NVU,NCODE,NPULS,NNV(32),IPRNT
COMMON/OVERLY/INITAL
DOUBLE PRECISION CALIB,NAMES
DATA CALIB/'CALIB '/,NAMES/'NAMES '/
COMMON/CALI/CAL(32)
C NNV(I)=1 IF CHAN USED, 0 IF NOT USED
C NV NUMBER OF CHANNELS ON DIGITIZED TAPE
C NVU NUMBER OF CHANNELS USED NVU.LE.NV
C
C USE Q AS TEMPORARY STORAGE IN READING
C READ PROBLEM CARD
C PROB - FIRST SIX CHARACTERS OF CONTROL CARD. IT SHOULD
C BE PROBLM OR FINISH.
C CODE - ALPHAMERIC PROBLEM CODE
C NSC - NUMBER OF SCANS
C NV - NUMBER OF VARIABLES, SERIES, OR CHANNELS
C NB - NUMBER OF FREQUENCY BANDS
C SAMPRT - SAMPLING RATE
C FILTER - INDICATES WHETHER PREFILTERING IS TO BE DONE
C ID - DECIMATION RATE
C L - HALF LENGTH OF PREFILTER
C PLTFLT - INDICATES WHETHER GRAPH OF PREFILTER IS DESIRED
C PLTPOW - INDICATES WHETHER GRAPHS OF POWER VERSUS FREQUENCY
C ARE DESIRED
C PLTLOG - INDICATES WHETHER GRAPHS OF LOGARITHM OF POWER
C VERSUS FREQUENCY ARE DESIRED
C SPEC - INDICATES WHETHER COMPLEX CROSS SPECTRA MATRIX IS
C TO BE WRITTEN ON FORTRAN TAPE 1 WITH FORMAT
C (8E14.7)
C PLTWND - INDICATES WHETHER A GRAPH OF THE FREQUENCY WINDOW
C OF THE POWER SPECTRUM ESTIMATES IS DESIRED
C IT - FORTRAN NUMBER FOR INPUT TAPE
C DPLTAP - INDICATES WHETHER INPUT DATA IS IN DPL/BRI FORMAT
C REWND - INDICATES WHETHER INPUT TAPE IS TO BE REWOUND BEFOR
C READING
C NF - NUMBER OF VARIABLE FORMAT CARDS
C IBEGIN - NUMBER OF SCANS TO BE SKIPPED BEFORE ANALYSIS
C BEGINS
C NF - NUMBER OF VARIABLE FORMAT CARDS FOR INPUT DATA
C LFT - LOGICAL TAPE FOR OUTPUT OF FILTERED DATA
C NFF - NUMBER OF VARIABLE FORMAT CARDS FOR LFT (OUTPUT)
C LCT - LOGICAL TAPE FOR OUTPUT OF COMPLEX CROSS SPECTRA MA
C
DATA IPROB,JORMSX,JFILT,JFORM/0,0,0,0/
DOUBLE PRECISION POSITN
DATA POSITN/'POSIT '/
LOGICAL IPRNT
DATA JTIT/0/,JBRI/0/,JORMSY/0/,JCHAN/0/
NTOT=5000
CALL USAGEB('BMDX92')
DATA JCALIB/0/
IBRI = 1
NWO=0
INITAL=0
C
C
C
IPRNT=.TRUE.
NWORDS=0
ICODE=0
DO 45 I=1,32
45 CAL(I)=1.0
DO 50I=1,29
50 IDDPL(I)=0
53 NID=29
IF (LABL.NE.FIN)GO TO 59
C FINISH CARD HAS BEENREAD
530 WRITE (6,1050)
1050 FORMAT('1FINISH CARD ENCOUNTERED. JOB TERMINATED'/ 1X)
511 IF(DPLTAP.EQ.YES)CALL POSIT
STOP
59 IF (IPROB.LE.0) GO TO 60
REWIND 1
WRITE(1,1060)LABL,TSV
1060 FORMAT(A6,18A4,A2)
REWIND 1
READ(1,1061)LABL,CODE,NSC,NV,NB,SAMPRT,FILT,PLTPOW,PLTLOG,
1PLTWND,IT,DPLTAP,REWND,NF,LCT,LBUNIT,IBEGIN,ISKIP
GO TO 61
60 READ (5,1061) LABL,CODE,NSC,NV,NB,SAMPRT,FILT,PLTPOW,PLTLOG,
1PLTWND,IT,DPLTAP,REWND,NF,LCT,LBUNIT,IBEGIN,ISKIP
1061 FORMAT(2A6,I6,2I3,F6.0,4A3,I2,A3,A2,I1,2I2,I6,I4)
IF (LABL .EQ. FIN) GO TO 530
61 IF(LABL.EQ.PROB) GO TO 62
WRITE (6,1062)LABL
1062 FORMAT('0',A6,' CARD ENCOUNTERED WHEN A PROBLEM CARD WAS EXPECTED
1, CARD WILL BE IGNORED.')
GO TO 60
62 MFLAG=0
IF (IT.LE.0)IT=5
FILTER=BNO
L=LZZ
ID=IDZZ
IF(ISKIP.EQ.0)ISKIP=2
IF(DPLTAP.EQ.YES)IT=13
IPROB=IPROB+1
IF(FILT.EQ.ORM.OR.FILT.EQ.LOW) FILTER=YES
IF (FILTER.EQ.YES)GO TO 70
NFF=0
ID=1
L=0
FILT=BNO
PLTFLT=BNO
70 READ(5,1060,END=79)LABL,TSV
REWIND 1
WRITE(1,1060)LABL,TSV
REWIND 1
IF(LABL.EQ.PROB)GO TO 200
IF(LABL.EQ.PRE)GO TO 80
IF(LABL.EQ.TIT)GO TO 90
IF(LABL.EQ.CHAN)GO TO 100
IF (LABL.EQ.NAMES) GO TO 1991
IF(LABL.EQ.CALIB)GO TO 141
IF(LABL.EQ.PAR)GO TO 110
IF(LABL.EQ.XVA)GO TO 120
IF(LABL.EQ.YVA)GO TO 130
IF(LABL.EQ.FIN)GO TO 200
IF(LABL.EQ.POSITN)GO TO 8140
C CARDS ASSUMED TO BE FORMAT CARDS
IF (NF.NE.0.OR.NFF.NE.0)GO TO 72
C
WRITE(6,1072)LABL,TSV
1072 FORMAT('0ERROR ON CARD FOUND'/5X,A6,18A4,A2/
1 ' WE CONTINUE WITH THE NEXT PROBLEM')
MFLAG=1
GO TO 70
72 JFORM=1
NWORDS=18*NF
IF(NF.LE.0) GO TO 74
C READ INPUT FORMAT
READ(1,1075) (Q(I),I=1,NWORDS)
IRR = 5
GO TO 75
74 IRR = 1
1075 FORMAT(18A4)
75 NWO=18*NFF
IF(NFF.LE.0)GO TO 78
I1=NWORDS+1
I2=NWORDS+NWO
READ(IRR ,1075) (Q(I),I=I1,I2)
78 IF(IT.NE.5)GO TO 70
IPROB=0
GO TO 200
79 LABL=FIN
GO TO 200
C IF PREFILTER CARD READ
80 JFILT=1
READ(1,1080) LABL,ID,L,NPOINT,PLTFLT,LFT,NFF
1080 FORMAT(A6,I3,I6,I2,A3,I2,I1)
IF(FILT.EQ.ORM)GO TO 82
IF(L.LE.0)L=12*ID
GO TO 84
82 IF(L.LE.0)L=50
84 LZZ=L
IDZZ=ID
GO TO 70
C IF TITLE CARD IS READ
90 IF(JTIT.LT.4)GO TO 92
1220 WRITE(6,1090)
1090 FORMAT('0MORE THAN FOUR ''TITLES'' WERE FOUND. PROGRAM USES THE '
1 'FIRST FOUR AND IGNORES THE OTHERS.')
GO TO 70
92 IF(JTIT.GT.0)GO TO 94
DO 93 I=1,68
93 TITLER(I)=BLANK
JTIT=0
94 I1=JTIT*17+1
I2=I1+16
JTIT=JTIT+1
READ(1,1070) LABL,(TITLER(I),I=I1,I2)
1070 FORMAT(A6,16A4,A2)
GO TO 70
C IF NAMCHN CARD IS READ
100 JCHAN=1
READ(1,1100) (NAMCHN(I),I=1,16)
IF (NV.LE.16) GO TO 70
READ(5,1100) (NAMCHN(I),I=17,NV)
GO TO 70
1991 JCHAN = 1
1100 FORMAT(6X,16A4)
READ(1,1101) (NAMCHN(I),I=1,8)
1101 FORMAT(6X,8A8)
IF(NV.LE.8) GO TO 70
READ (5,1101) (NAMCHN(I),I=9,NV)
GO TO 70
C IF BRI PARAM CARD IS READ
110 JBRI=1
READ(1,1110) LABL,LONG,BCODE,HCODE,LENGTH,BPULSE,HPULSE,EXPNO,
1 PSTART,PEND,TSTART,TEND,IDAY,IHR,IMIN,ISEC,LBTAPE,ID1BRI,ID2BRI
2,FILE
1110 FORMAT(A6,I3,2I5,2I3,I5,3I3,2F7.0,I3,3I2,A6,2I3,A4,A2)
FS=SAMPRT
IF (LONG .EQ. 0) LONG = 30
IF (LENGTH .EQ. 0) LENGTH = 2
GO TO 70
C IF XVALUE CARD IS READ
120 JORMSX=1
READ(1,1120) LABL,(X(I),I=1,11)
1120 FORMAT(A6,11F6.0)
GO TO 70
C IF YVALUE CARD IS READ
130 JORMSY=1
READ(1,1120) LABL,(Y(I),I=1,11)
GO TO 70
C IF BRI POSIT PARAM IS READ
8140 READ(1,1110)LABL,LONG,BCODE,HCODE,LENGTH,BPULSE,HPULSE,EXPNO,
1PSTART,PEND,TSTART,TEND,IDAY,IHR,IMIN,ISEC,LBTAPE,ID1BRI,ID2BRI
2,FILE
FS=SAMPRT
IF (LONG .EQ. 0) LONG = 30
IF (LENGTH .EQ. 0) LENGTH = 2
CALL POSIT
CCCC WHAT TO DO ABOUT ERROR RETURNS
GO TO 70
141 READ(1,1140) (CAL(I),I=1,NV)
JCALIB=1
DO 142 I=1,NV
IF(CAL(I).EQ.0.0)CAL(I)=1.0
142 CONTINUE
1140 FORMAT(6X,16F4.0)
GO TO 70
C CHECKING FOR CONTROL CARDS AND CHECKING PARAMETERS
200 IF (MFLAG.EQ.1)GO TO 53
IF (FILTER.EQ.YES.AND.JFILT.EQ.0)GO TO 219
IF(FILT.NE.ORM)GO TO 216
IF(JORMSX.EQ.0.OR.JORMSY.EQ.0)GO TO 218
IF(JORMSX.EQ.2.AND.JORMSY.EQ.1)GO TO 218
IF(JORMSY.EQ.2)GO TO 8212
IF(NPOINT.LE.0.OR.NPOINT.GT.11)NPOINT=11
210 NPOINT=NPOINT-1
IF(Y(NPOINT).EQ.0.)GO TO 210
NPOINT=NPOINT+1
IF(NPOINT.LE.1)GO TO 9920
JORMSY=2
8212 IF(JORMSX.EQ.2)GO TO 216
DO 8214 I=1,NPOINT
8214 X(I)=X(I)/SAMPRT
IF(ID.LE.0)ID=MAX1(0.5/X(NPOINT)+0.001,1.0)
JORMSX=2
216 GO TO 220
218 WRITE(6,1218)
1218 FORMAT('1BOTH AN XVALUE AND YVALUE CARD WERE NOT PRESENT FOR A
1CASE SPECIFYING AN ORMSBY FILTER'/' THIS PROBLEM WILL BE SKIPPED')
GO TO 53
219 WRITE(6,1219)
1219 FORMAT('1NO PREFILTER CARD WAS FOUND FOR A PROBLEM SPECIFYING
1EITHER AN ORMSBY OR LOW PASS FILTER'/' THIS PROBLEM WILL BE SKIPPE
2D')
GO TO 53
220 NTIT=IABS(JTIT)
JTIT=-NTIT
IF(JCHAN.GT.0)GO TO 226
IF(DPLTAP.NE.YES)GO TO 222
WRITE (6,8220)
8220 FORMAT ('1NO NAMCHN CARD WAS PRESENT WHEN DPLTAPE WAS SPECIFIED'/'
1 THIS PROBLEM WILL BE SKIPPED')
GO TO 53
226 NVU=MIN0(NVU,NV)
IF(JCHAN.EQ.2)GO TO 8227
JCHAN=2
NVU=0
NCODE=0
NPULS=0
DO 224 I=1,NV
IF(NAMCHN(I).EQ.CADE.OR.NAMCHN(I).EQ.STIM.OR.NAMCHN(I).EQ.NULL.OR.
1NAMCHN(I).EQ.BLANK) GO TO 223
NNV(I)=1
NVU=NVU+1
NAMCHN(NVU)=NAMCHN(I)
GO TO 224
223 NNV(I)=0
IF (NAMCHN(I) .NE. CADE ) GO TO8226
NCODE=I
GO TO 224
8226 IF(NAMCHN(I).NE.STIM) GO TO 224
NPULS=I
224 CONTINUE
IF(NPULS.EQ.0) NPULS=NCODE
GO TO 8227
222 DO8224 I=1,NV
NNV(I)=I
8224 NAMCHN(I)=NAME(I)
NVU=NV
8227 IF (DPLTAP.NE.YES.OR.JBRI.EQ.1) GO TO 230
WRITE (6,1230)
1230 FORMAT('1NO PARAM CARD FOUND WHEN DPLTAPE WAS SPECIFIED'/' PROBLE
1M WILL BE SKIPPED')
230 IF(JFORM.GT.0.OR.DPLTAP.EQ.YES)GO TO 232
WRITE (6,1232)
1232 FORMAT('1NO FORMAT CARDS FOUND WHEN BCD INPUT FROM CARDS OR TAPE
1WAS EXPECTED'/' THIS PROBLEM WILL BE SKIPPED')
GO TO 53
232 CONTINUE
250 NBT=2
252 NBT=2*NBT
IF (NB.GT.NBT)GO TO 252
251 NB=NBT
IF (DPLTAP.NE.YES) GO TO 2516
IF(PEND.EQ.0.AND.TEND.EQ.0.)WRITE(6,9912)
9912 FORMAT('0ON PARAMR CARD BOTH PEND AND TEND ARE EQUAL ZERO. ' /
1' THIS PROBLEM POSITIONS THE BRI TAPE ONLY.')
NSC=0
IF(PEND.EQ.0)NSC=TEND*SAMPRT
2516 IF(FILT.EQ.LOW)ID=MAX0(ID,2)
IF(PLTFLT.NE.YES) PLTFLT=BNO
IF(NSC.EQ.0)GO TO 2540
NS=(NSC-(2*L+1))/ID+1
LOG2NS=0
NSCANS=1
2544 LOG2NS=LOG2NS+1
NSCANS=NSCANS*2
IF(NS.GT.NSCANS) GO TO 2544
NB=MIN0(NSCANS/4,NB)
C PRINT TITLES AND SUMMARY OF PROBLEM CARD
2540 WRITE(6,145)
145 FORMAT('1')
WRITE(6,101)CODE,NV,NVU,NB
101 FORMAT
1(' BMDX92 - JET SPEED SPECTRUM ESTIMATION - REVISED ',
2'JULY 29, 1968'
X / 1X, 5X, 40HHEALTH SCIENCES COMPUTING FACILITY, UCLA
X // 1H0, 12HPROBLEM CODE, 36X, A6,
X / 1X, 5X, 16HNUMBER OF SERIES, 27X, I6,
X/6X,'NUMBER OF SERIES TO BE ANALYZED' 12X,I6,
X / 1X, 5X, 25HNUMBER OF FREQUENCY BANDS, 18X, I6 )
IF(NSC.EQ.0) GO TO 1010
WRITE(6,1013)NSC,NS
1013 FORMAT('0',5X,'NUMBER OF SCANS READ', 23X,I6/6X,
1 'NUMBER OF SCANS AFTER DECIMATION' ,11X,I6)
1010 IF(NTIT.LE.0)GO TO 1420
NTITA=17*NTIT
WRITE(6,1011)(TITLER(I),I=1,NTITA)
1011 FORMAT(1H0 / / 4(15X,17A4 /))
WRITE(6,1012)
1012 FORMAT (1H0 / / )
1420 IF (SAMPRT .LT. .00001) SAMPRT = 1.0
WRITE (6, 143) SAMPRT
143 FORMAT (1X, 5X, 13HSAMPLING RATE, 24X, F12.5)
WRITE(6,131)FILT
131 FORMAT (1H0, 5X, 27HPREFILTERING OF EACH SERIES, 16X, 3X, A3)
IF (FILTER .NE. YES) GO TO 211
WRITE (6, 132) ID, L, PLTFLT
132 FORMAT
X ( 1X, 10X, 15HDECIMATION RATE, 23X, I6
X / 1X, 10X, 28HHALF LENGTH OF THE PREFILTER, 10X, I6
X / 1X, 10X, 31HGRAPH OF THE FREQUENCY RESPONSE
X / 1X, 10X, 27H FUNCTION OF THE PREFILTER, 11X, 3X, A3 )
IF(FILT.NE.ORM)GO TO 211
WRITE(6,9211)(Y(I),I=1,NPOINT)
WRITE(6,9212)(X(I),I=1,NPOINT)
9211 FORMAT('0THE ORMSBY FILTER SPECIFIED WAS'/5X,
1 ' Y VALUES' 11F8.4)
9212 FORMAT(5X,' X VALUES' 11F8.4)
WRITE(6,9213)
9213 FORMAT(10X,'THE X VALUES HAVE BEEN DIVIDED BY THE SAMPLING RATE')
211 CONTINUE
SPEC=BNO
IF(LCT.GT.0)SPEC=YES
IF (PLTPOW .NE. YES) PLTPOW = BNO
IF (PLTLOG .NE. YES) PLTLOG = BNO
IF (PLTWND .NE. YES) PLTWND = BNO
WRITE(6,133) PLTPOW,PLTLOG,LCT,SPEC,PLTWND
133 FORMAT
X ( 1H0, 5X, 36HGRAPHS OF POWER VERSUS FREQUENCY FOR
X / 1X, 5X, 13H AUTOSPECTRA, 30X, 3X, A3
X / 1X, 5X, 35HGRAPHS OF LOGARITHM OF POWER VERSUS
X / 1X, 5X, 27H FREQUENCY FOR AUTOSPECTRA, 16X, 3X, A3
X / 1X, 5X, 36HMATRICES OF CROSS SPECTRA IN COMPLEX
X / 1X, 5X, ' FORM WRITTEN ON FORTRAN TAPE',I3,
X / 1X, 5X, ' USING FORMAT (30A4)',23X,A3,
X / 1X, 5X, 32HGRAPH OF THE FREQUENCY WINDOW OF
X / 1X, 5X, 30H THE POWER SPECTRUM ESTIMATES, 13X, 3X, A3)
IBEGIN = MAX0(IBEGIN, 0)
IF(JCHAN.GT.0)WRITE(6,9134)(NAMCHN(I),I=1,NVU)
9134 FORMAT('0 THE ANALYZED SERIES NAMES ARE'/(15X,8A12))
IF(JCALIB.GT.0)WRITE(6,9135)(CAL(I),I=1,NV)
9135 FORMAT('0 THE SCALING FACTORS FOR EACH SERIES ARE'/(15X,12F7.3
1))
IF (DPLTAP .EQ. YES) GO TO 212
IF (IT .LE. 0) IT = 5
NF = MAX0(NF, 1)
IF (IT .NE. 5) GO TO 213
C IF DATA INPUT IS FROM CARDS THE NEXT 9 LINES APPLY
WRITE (6, 134)
134 FORMAT (1H0, 5X, 19HINPUT IS FROM CARDS)
WRITE (6, 135) IBEGIN
135 FORMAT
X ( 1X, 10X, 29HNUMBER OF SCANS TO BE SKIPPED
X / 1X, 10X, 24H BEFORE ANALYSIS BEGINS, 14X, I6 )
WRITE (6, 137) NF
137 FORMAT(11X,'NUMBER OF VARIABLE FORMAT CARDS FOR INPUT',7X,I6)
GO TO 214
C IF DATA INPUT IS FROM BCD TAPE THE NEXT 12 LINES APPLY
213 CONTINUE
IF (REWND .NE. ONO) REWND = YES
WRITE (6, 136) IT, REWND
136 FORMAT
X ( 1H0, 5X, 22HINPUT IS FROM BCD TAPE
X / 1X, 10X, 27HINPUT TAPE (FORTRAN NUMBER), 11X, I6
X / 1X, 10X, 33HINPUT TAPE REWOUND BEFORE READING, 5X, 3X, A3 )
IF (REWND .EQ. YES) REWIND IT
WRITE (6, 135) IBEGIN
WRITE (6, 137) NF
GO TO 214
212 CONTINUE
C IF DATA INPUT IS FROM BRI/DPL DIGITIZED TAPE THE NEXT 10
C LINES APPLY
NF = 0
WRITE (6, 138)
138 FORMAT
X ( 1H0, 5X, 35HINPUT IS FROM TAPE DIGITIZED AT DPL )
WRITE (6, 139)
139 FORMAT
X ( 1X, 10X, 40HDPL TAPE FORMAT IS ASSUMED FOR THE INPUT
X / 1X, 10X, 35H TAPE, SO NO FORMAT CARDS ARE READ)
WRITE(6,140)LONG,BCODE,HCODE,LENGTH,BPULSE,HPULSE,EXPNO,FILE,
1 PSTART,PEND,TSTART,TEND,IDAY,IHR,IMIN,ISEC,LBUNIT,LBTAPE
140 FORMAT
X(11X,'LONG = ',I3,' BCODE = ',I5,' HCODE = ',I5,' LENGTH = ',I3
X/11X,'BPULSE = ',I3,' HPULSE = ',I5,' EXPNO = ',I3,' FILE = ',
XA4,A2
X/11X,'PSTART = ',I3,' PEND = ',I3,' TSTART = ',F7.3,' TEND = '
X,F7.3/11X,'IDAY = ',I3,' IHR = ',I2,' IMIN = ',I2,' ISEC = ',
XI2,' LBUNIT = ',I2,' LBTAPE = ',A6)
214 CONTINUE
IF(NF.LE.0) GO TO 2140
WRITE(6,1111) VIN, (Q(I),I=1,NWORDS)
1111 FORMAT('0VARIABLE',A6,' FORMAT IS' / (5X,18A4))
2140 IF(NFF.LE.0) GO TO 2148
WRITE(6,2141) LFT
2141 FORMAT(1X,10X,'LOGICAL TAPE FOR OUTPUT OF FILTERED DATA',7X,I3)
WRITE(6,2142) NFF
2142 FORMAT(1X,10X,'NUMBER OF VARIABLE FORMAT CARDS FOR OUTPUT OF '
1 'FILTERED DATA',5X,I3)
I1=NWORDS+1
I2=NWORDS+NWO
WRITE(6,1111) VOUT, (Q(I),I=I1,I2)
2148 CONTINUE
CALL BEANS
IF( ISIZE2 .GT. NTOT ) GO TO 58
IF((ISTART+NCHANS*NSCANS).GT.NTOT)GO TO 58
CALL INPUT (JSCANS,Q, Q(ISTART), Q, Q(IWTS), Q(ITEMP), Q(ISTART),
1 IRR,LOGIRR,Q(ITEMP),Q(IFT2))
IF(ICODE.LT.0) GO TO 53
IF(ICODE.GT.0)GO TO 511
CALL OUTPUT (Q, IROWSP, Q(IAMPL), Q(ILOGA), Q(ISPEC), Q(IBRI))
IF (PLTWND .NE. YES) GO TO 708
C IF A GRAPH OF THE FREQUENCY WINDOW OF THE POWER SPECTRUM
C ESTIMATES IS REQUESTED SUBROUTINE WNDPLT PRODUCES THIS
C GRAPH
CALL WNDPLT ( NSCANS, Q )
708 CONTINUE
GO TO 53
58 CONTINUE
WRITE(6,112)NTOT,ISTART,NCHANS,NSCANS
112 FORMAT('1THE SPACE TO STORE THE DATA IS GREATER THAN ' I8,'WORDS'/
1 '0THE FIRST STORAGE LOCATION FOR THE DATA IS' I8,/
2 '0THE NUMBER OF CHANNELS IS 'I4,5X,' AND THE NUMBER OF SCANS IS'
2 I6/'0THE NEXT PROBLEM WILL BE TRIED' )
MFLAG=1
GO TO 53
9920 WRITE(6,9921)
9921 FORMAT('0ORMSBY FILTER SPECIFIED ON PROBLM CARD.BUT NUMBER OF '
1 'POINTS TO SPECIFY FILTER IS ZERO.'/' PROBLEM OMITTED.')
MFLAG=1
GO TO 70
END
SUBROUTINE INPUT ( JSCANS, P, SCANS, FMT, WTS, TEMP, FILT, IRR,
1 LOGIRR,ITEMP,FMT2)
COMMON /LIST/ NSC, LEN4, MSCANS, NV, NCHANS, NB, SR, FILTER, ID,
1 L, PLTFLT, SPEC, PLTWND, IT, DPLTAP, IBEGIN, NF, NWORDS, NWTS,
2 SAMPRT, LOG2NS, PLTPOW, PLTLOG, NS, NSCANS,NFF,LFT,LCT
COMMON/PREF/X(11),Y(11),NPOINT
DOUBLE PRECISION NAMCHN
COMMON/CHANEL/NAMCHN(32),NVU,NCODE,NPULS,NNV(32)
COMMON/XDATA/ICODE
C
COMMON/BRI/FILE(2),LONG,BCODE,HCODE,LENGTH,BPULSE,HPULSE,EXPNO,
1 PSTART,PEND,TSTART,TEND,IDAY,IHR,IMIN,ISEC,FS,ISKIP
INTEGER BCODE,HCODE,BPULSE,HPULSE,EXPNO,PSTART,PEND
DIMENSION SCANS(JSCANS,1),P(1),FMT(1),WTS(1),TEMP(1),FMT2(1)
DIMENSION FILT(IRR, 1)
DIMENSION ITEMP(1)
DIMENSION ITEMPO(20)
DATA YES / 3HYES /
DATA PI / 3.1415927 /
C MSCANS - MAXIMUM NO OF SCANS THAT CAN BE READ IN
C JSCANS - MAXIMUM NO OF SCANS THAT CAN BE STORED - POWER OF 2
IF(DPLTAP.EQ.YES) GO TO 20
J=0
DO 10 I=1,NV
IF(NNV(I).LE.0) GO TO 10
J=J+1
ITEMPO(J)=I
10 CONTINUE
20 IF(FILTER.NE.YES) GO TO 603
IF(NPOINT.GT.0) GO TO 2611
C IF PREFILTERING IS DESIRED, SUBROUTINE GENWTS GENERATES THE
C WEIGHTS USED FOR PREFILTERING
CALL GENWTS(ID, L, WTS)
GO TO 701
2611 CONTINUE
CALL GENFLT (WTS(L+1), L, X, Y, NPOINT, ITEMP)
LL = 2 * L + 2
DO 700 I=1,L
700 WTS(I) = WTS(LL-I)
C IF A GRAPH OF THE FREQUENCY RESPONSE FUNCTION OF THE
C PREFILTER IS REQUESTED THE FOLLOWING STATEMENTS UP TO
C ' 707 CONTINUE' PRODUCE THAT GRAPH.
701 IF(PLTFLT .NE. YES) GO TO 707
DO 703 I1 = 1, IRR
FILT(I1, 1) = 0.0
703 FILT(I1, 2) = 0.0
LP1 = L + 1
DO 704 I1 = 1, LP1
I3 = L + I1
704 FILT(I1, 1) = WTS(I3)
DO 705 I1 = 1, L
I3 = IRR - L + I1
705 FILT(I3, 1) = WTS(I1)
C SUBROUTINE FAST OBTAINS THE FINITE FOURIER TRANSFORM OF
C THE FILTER WEIGHTS
CALL FAST ( FILT(1, 1), FILT(1, 2), LOGIRR)
IRR2 = IRR / 2 + 1
IFLAG = 3
C SUBROUTINE PLOT1 GRAPHS THE TRANSFORMED FILTER WEIGHTS TO
C PRODUCE THE GRAPH OF THE PREFILTER. IFLAG = 3 PRODUCES
C THE APPROPRIATE TITLE ON THE GRAPH AND THE CORRECT BASE
C LINE FOR THE GRAPH
CALL PLOT1( FILT(1, 1), I, IRR2, SAMPRT, IFLAG)
707 CONTINUE
603 CONTINUE
IF(DPLTAP.EQ.YES) GO TO 261
IF (IBEGIN .EQ. 0) GO TO 84
C READ SCANS TO BE SKIPPED BEFORE PREFILTERING (OR ANALYSIS,
C IF NO PREFILTERING IS DONE) STARTS
IF (DPLTAP .EQ. YES) GO TO 261
DO 72 I1 = 1, IBEGIN
READ(IT,FMT,END=996,ERR=992) (SCANS(1,I2),I2=1,NV)
72 CONTINUE
GO TO 262
261 CONTINUE
CALL POSIT
C ITEMPO AND TEMP MAY NOT BE EQUIV FOR POSITP
IF(IERROR.GT.0) RETURN
IF(PEND.EQ.0.AND.TEND.EQ.0)GO TO 994
262 CONTINUE
84 CONTINUE
C
C IF PREFILTERING IS DONE THE FOLLOWING STATEMENTS DOWN TO
C ' 206 CONTINUE' ARE USED TO READ THE DATA AND PREFILTER AT
C THE SAME TIME
DO 201 I1=1,JSCANS
DO 201 I2 = 1, NCHANS
201 SCANS(I1, I2) = 0.0
IF(FILTER.NE.YES) GO TO 601
KSCAN = 1
INDEX1 = 1
DO 206 I3=1,MSCANS
IF (DPLTAP .EQ. YES) GO TO 267
READ(IT,FMT,END=210,ERR=996) (TEMP(I2),I2=1,NV)
DO 2010 I=1,NVU
J=ITEMPO(I)
2010 TEMP(I)=TEMP(J)
GO TO 268
267 CALL POSIT
IF (ICODE .NE. 0) GO TO 99
IF(IFSTOP.GT.0) GO TO 210
268 CONTINUE
I1 = KSCAN
INDEX2 = INDEX1
203 DO 204 I2=1,NVU
204 SCANS(I1, I2) = SCANS(I1, I2) + TEMP(I2)*WTS(INDEX2)
I1 = I1 +1
IF (I1 .GT. NS) GO TO 205
INDEX2 = INDEX2 - ID
IF (INDEX2 .GT. 0) GO TO 203
205 INDEX1 = INDEX1 + 1
IF (INDEX1 .LE. NWTS) GO TO 206
INDEX1 = INDEX1 - ID
KSCAN = KSCAN + 1
206 CONTINUE
NV=NVU
IF(NSC.EQ.0) GO TO 2100
WRITE(6,207) NSC
207 FORMAT('0NUMBER OF SCANS READ FOR THIS PROBLEM WAS',I5)
GO TO 55
2100 WRITE(6,2101) MSCANS
2101 FORMAT('0MAXIMUM NUMBER OF SCANS (',I5,') HAVE BEEN READ FOR THIS
1PROBLEM BEFORE STOP CONDITION ENCOUNTERED.'/' REST OF DATA ARE IGN
2ORED.')
GO TO 600
210 IF(NSC.NE.0) GO TO 212
NS=KSCAN-1
212 CONTINUE
NV=NVU
IF(NSC)600,600,55
601 CONTINUE
IF (FILTER .EQ. YES) GO TO 600
C IF PREFILTERING IS NOT DONE THE NEXT STATEMENTS ARE USED
C TO READ THE DATA
263 CONTINUE
DO 265 I1=1,MSCANS
IF(DPLTAP.EQ.YES) GO TO 2610
READ(IT,FMT,END=264,ERR=996)(TEMP(I2),I2=1,NV)
DO 2609 I2=1,NVU
J=ITEMPO(I2)
2609 SCANS(I1,I2)=TEMP(J)
GO TO 265
C
2610 CALL POSIT
IF (ICODE .NE. 0) GO TO 99
IF(IFSTOP.GT.0) GO TO 264
DO 266 I2 = 1, NVU
266 SCANS(I1,I2)=TEMP(I2)
265 CONTINUE
NV=NVU
C
IF (NSC.EQ.0) GO TO 2652
WRITE(6,207) NSC
GO TO 55
2652 WRITE (6,2101) JSCANS
C
GO TO 600
264 NS=I1
NV = NVU
IF(NSC.GT.0) GO TO 55
600 CONTINUE
C NS - EFFECTIVE NO OF SCANS READ IN.
C NSCANS IS THE LEAST POWER OF 2 GREATER THAN OR EQUAL TO NS,
C (NSCANS IS LOWER CASE M IN THE WRITE UP)
LOG2NS = 0
NSCANS = 1
54 IF (NS .LE. NSCANS) GO TO 55
LOG2NS = LOG2NS + 1
NSCANS = 2 * NSCANS
GO TO 54
55 CONTINUE
C NB SHOULD BE GREATER THAN OR EQUAL TO 4 AND LESS THAN OR
C EQUAL TO NSCANS/4
NSD4 = NSCANS / 4
IF(NB.GT.NSD4)WRITE(6,2102)NB,NSCANS
2102 FORMAT('1THE NO. OF BANDS SPECIFIED',I5,5X,'FOR THIS PROBLEM IS GR
1EATER THAN',I5,'/4. THE RESULTS MAY BE IN ERROR')
NB = MIN0(NB, NSD4)
ITEST = (2*NV - 1) * (NB + 1)
IF (ITEST .LE. NSCANS) GO TO 56
LEN1 = ( (2*NV - 1) * (NB + 1) - NSCANS)**2
LEN1 = LEN1 / (4 * (NB + 1) )
LEN1 = LEN1 + 2 * NV * (NB + 1) + 1
IF(LEN1.LE.LEN4) GO TO 56
WRITE(6,999)
999 FORMAT(35H0PROBLEM EXCEEDS AVAILABLE STORAGE.)
STOP
56 CONTINUE
IF(NFF.EQ.0) GO TO 602
WRITE(LFT,FMT2)((SCANS(I1,I2),I2=1,NV),I1=1,NS)
WRITE (6,560) LFT, NV, NS
560 FORMAT('0FILTERED DATA WRITTEN ON UNIT',I3 /
1 ' NUMBER OF SERIES ', I2, ', NUMBER OF SCANS ',I5)
602 IF (NS .EQ. NSCANS) GO TO 74
C ZERO OUT THE LAST PART OF EACH SERIES IF THE NUMBER OF SCANS
C IS NOT A POWER OF 2
I1BEGN = NS + 1
DO 75 I1 = I1BEGN, NSCANS
DO 75 I2 = 1, NV
SCANS(I1, I2) = 0.0
75 CONTINUE
74 CONTINUE
C DETREND EACH SERIES BY SUBTRACTING FROM EACH SERIES ITS
C LEAST SQUARES LINEAR REGRESSION LINE
FNS = NS
TBAR = 0.5 * (FNS + 1.0)
TSUMSQ = (FNS * (FNS + 1.0) * (2.0*FNS + 1.0) ) / 6.0
DO 76 I2 = 1, NV
SUM = 0.0
CRSPRO = 0.0
DO 77 I1 = 1, NS
SUM = SUM + SCANS(I1, I2)
CRSPRO = CRSPRO + FLOAT(I1) * SCANS(I1, I2)
77 CONTINUE
FMEAN = SUM / FNS
BETA = (CRSPRO - FNS*TBAR*FMEAN) / (TSUMSQ - FNS*TBAR*TBAR)
DO 78 I1 = 1, NS
FREG = FMEAN + BETA*(FLOAT(I1) - TBAR)
SCANS(I1, I2) = SCANS(I1, I2) - FREG
78 CONTINUE
76 CONTINUE
C WINDOW EACH SERIES WITH A COSINE TAPER
IR = NS / 10
R = IR
DO 79 I1 = 1, IR
FI1 = I1
FINT = FI1 - 0.5
WINDOW = 0.5 * (1.0 - COS(PI * FINT / R ) )
I3 = NS+ 1 - I1
DO 80 I2 = 1, NV
SCANS(I1, I2) = WINDOW * SCANS(I1, I2)
SCANS(I3, I2) = WINDOW * SCANS(I3, I2)
80 CONTINUE
79 CONTINUE
IF (NCHANS .EQ. NV) GO TO 82
C IF THE NUMBER OF SERIES IS ODD, FILL A DUMMY SERIES WITH
C ZEROS
DO 83 I1 = 1, NSCANS
SCANS(I1, NCHANS) = 0.0
83 CONTINUE
82 CONTINUE
CALL TRANS ( P, SCANS, NSCANS, JSCANS, NV, NB, LOG2NS )
C THE CROSS SPECTRUM ESTIMATES IN ARRAY P ARE SCALED BY
C MULTIPLYING BY C
WNDPWR = FNS - 1.25 * R
FSCANS = NSCANS
FNB = NB
FD = FSCANS / (2.0 * FNB)
C = 0.25 / (SR * (FD + 1.0) * WNDPWR)
IROWSP = 2 * (NB + 1)
ICOLSP = (NV * (NV + 1) ) / 2
ISIZEP = IROWSP * ICOLSP
DO 85 I1 = 1, ISIZEP
P(I1) = C * P(I1)
85 CONTINUE
RETURN
99 WRITE(6,990) ICODE
990 FORMAT ('0READING TERMINATED. PICKUP SET ICODE = ' , I8)
RETURN
992 WRITE (6,993)
993 FORMAT('0END FILE FOUND ON POSITIONING INPUT TAPE. PROBLEM OMITTED
1.')
994 ICODE=-1
RETURN
996 WRITE (6,997)
997 FORMAT('0ERROR ON INPUT TAPE. PROBLEM OMITTED.')
GO TO 994
END
SUBROUTINE OUTPUT(P,IROWSP,AMPL,ALOGAM,COMSPC,DPDATA)
COMMON /LIST/ NSC, LEN4, MSCANS, NV, NCHANS, NB, SR, FILTER, ID,
1 L, PLTFLT, SPEC, PLTWND, IT, DPLTAP, IBEGIN, NF, NWORDS, NWTS,
2 SAMPRT, LOG2NS, PLTPOW, PLTLOG, NS, NSCANS,NFF,LFT,LCT
DOUBLE PRECISION LBTAPE
COMMON/BRIID/LBTAPE,LBUNIT,IDDPL(29),NID,TITLER(68),ID1BRI,ID2BRI
COMMON/BRI/FILE(2),LONG,BCODE,HCODE,LENGTH,BPULSE,HPULSE,EXPNO,
1 PSTART,PEND,TSTART,TEND,IDAY,IHR,IMIN,ISEC,FS,ISKIP
INTEGER BCODE,HCODE,BPULSE,HPULSE,EXPNO,PSTART,PEND
DIMENSION DPDATA(1)
DIMENSION P(IROWSP, 1), AMPL(1), ALOGAM(1)
DIMENSION COMSPC(1)
DIMENSION AM(13),CO(13),PH(5),II(13),JJ(5)
DOUBLE PRECISION NAMCHN,INAME,JNAME
DIMENSION PHASE(2)
DIMENSION INAME(13),JNAME(13)
DIMENSION XAMCHN(32)
COMMON/CALI/CAL(32)
COMMON/CHANEL/NAMCHN(32),NVU,NCODE,NPULS,NNV(32)
EQUIVALENCE(NAMCHN,XAMCHN)
DATA PHASE/4HPHA ,4HSE /
DATA FREQ,AMP,COH/4HFREQ,4HAMP ,4HCOH /
DATA YES/4HYES /
DATA BLANK/4H /
C DATA NCASE/0/
ICOLSP = (NV * (NV + 1) ) / 2
DATA PI / 3.14159265 /
C DF = DEGREES OF FREEDOM
DF=FLOAT(NS)/FLOAT(NB)
NBAND=MIN0(NB,35)
LMAX = NB + 1
IF(LBUNIT.LE.0) GO TO 50
NID=29
180 FORMAT('0CASE WILL BE WRITTEN IN MERGE FORMAT ON TAPE ' A6,
1' ON LOGICAL',I3,'.'/' ID VECTOR IS'/' ',22I4,I7,7I4)
IDDPL(1)=1
IDDPL(2)=1
IDDPL(3)=ID1BRI
IDDPL(4)=ID2BRI
IDDPL(5)=DF
C NCASE=NCASE+1
IDDPL(6)=IDDPL(6)+1
IDDPL(7)=SAMPRT
IDDPL(8)=NV
IDDPL(9) =-1
IDDPL(10)=EXPNO
IDDPL(11)=PSTART
IDDPL(13)=TSTART
IDDPL(14)=PEND
IDDPL(16)=TEND*100.0
IDDPL(17)=(SR/FLOAT(NB*2))*(LMAX-1)
IDDPL(18)=(SR/FLOAT(NB*2))*1000.
IDDPL(19)=LBUNIT
IDDPL(20)=0
IDDPL(21)=0
IDDPL(22)=IDAY
IDDPL(23)=IHR*10000+IMIN*100+ISEC
IDDPL(24)=2
IDDPL(25)=ID
WRITE(6,180)LBTAPE,LBUNIT,IDDPL
C SPECTRA FROM X92 * 2
NDP=456
DO 199 I=1,NDP
199 DPDATA(I)=BLANK
K=95
DO 200 I=1,NV
DPDATA(K+1)=XAMCHN(2*I)
DPDATA(K)=XAMCHN(2*I-1)
200 K=K+2
C
DO 22 I=1,17
22 DPDATA(I)=TITLER(I)
C
K=433
DO 24 I=18,34
DPDATA(K)=TITLER(I)
24 K=K+1
WRITE(LBUNIT)NID,IDDPL,NDP,(DPDATA(I),I=1,NDP)
NDP=0
50 CONTINUE
C THE FOLLOWING STATEMENTS DOWN TO ' 10 CONTINUE' PRODUCE
C THE TABLE OF THE POWER AT EACH FREQUENCY BAND FOR EACH
C SERIES
K = 1
WRITE (6, 101) DF
101 FORMAT (1H1, 20X, 27HAUTOSPECTRA FOR EACH SERIES, 10X, F9.2,
X 1X, 18HDEGREES OF FREEDOM )
WRITE (6, 113)
113 FORMAT (1H0, / 1X, 57HTABLE OF THE POWER AT EACH FREQUENCY BAND FO
XR EACH SERIES )
DO 10 I = 1, NV
II(K) = I
INAME(K)=NAMCHN(I)
IF((K.EQ.13) .OR. (I.EQ.NV)) GO TO 11
K = K + 1
GO TO 10
11 CONTINUE
WRITE(6,114) (INAME(L),L=1,K)
114 FORMAT('0'/' FREQ',3X,13A9)
DO 15 L = 1, LMAX
FR = SR * FLOAT(L - 1) / FLOAT(NB * 2)
IF(LBUNIT.LE.0.OR.L.GE.52)GO TO 115
IDDPL(21)=IDDPL(21)+1
NDP=NDP+1
DPDATA(NDP)=FR
115 CONTINUE
DO 14 M = 1, K
NIJ = (II(M) - 1) * NV + II(M) - ( (II(M) - 1) * II(M) ) / 2
C = P(2 * L - 1, NIJ)
Q = P(2 * L, NIJ)
AM(M)=SQRT(C**2+Q**2)*CAL(II(M))**2
IF(LBUNIT.LE.0.OR.L.GE.52)GO TO 14
NDP=NDP+1
DPDATA(NDP)=AM(M)*2.0
14 CONTINUE
WRITE (6, 104) FR, (AM(M1), M1 = 1, K)
104 FORMAT(1X,F8.3,13F9.2)
15 CONTINUE
K = 1
10 CONTINUE
IF(LBUNIT.LE.0) GO TO 130
IDDPL(1)=2
C IDDPL(21)=LMAX
WRITE(LBUNIT)NID,IDDPL,NDP,(DPDATA(I),I=1,NDP)
130 CONTINUE
IF ( (PLTPOW .NE. YES) .AND. (PLTLOG .NE. YES) ) GO TO 19
C THE FOLLOWING STATEMENTS DOWN TO ' 19 CONTINUE' PRODUCE THE
C GRAPHS OF POWER AND LOGARITHM OF POWER VERSUS FREQUENCY
C FOR EACH SERIES
DO 16 I = 1, NV
DO 18 L = 1, LMAX
NIJ = (I - 1) * NV + I - ((I - 1) * I) / 2
C = P(2*L - 1, NIJ)
Q = P(2*L, NIJ)
AMPL(L) = SQRT(C**2 + Q**2)
18 CONTINUE
IF (PLTLOG .NE. YES) GO TO 31
DO 32 L = 1, LMAX
32 ALOGAM(L) = ALOG(AMPL(L))
31 CONTINUE
IF (PLTPOW .NE. YES) GO TO 20
C SUBROUTINE PLOT1 PRODUCES THE GRAPH OF POWER VERSUS
C FREQUENCY FOR SERIES I. IFLAG = 1 PRODUCES THE
C APPROPRIATE TITLE AND CORRECT BASE LINE
IFLAG = 1
CALL PLOT1(AMPL, I, LMAX, SR, IFLAG)
20 CONTINUE
IF (PLTLOG .NE. YES) GO TO 21
C SUBROUTINE PLOT1 PRODUCES THE GRAPH OF LOGARITHM OF POWER
C VERSUS FREQUENCY FOR SERIES I. IFLAG = 2 PRODUCES THE
C APPROPRIATE TITLE AND THE CORRECT BASE LINE
IFLAG = 2
CALL PLOT1(ALOGAM, I, LMAX, SR, IFLAG)
21 CONTINUE
16 CONTINUE
19 CONTINUE
IF ( NV .EQ. 1 ) GO TO 29
C THE FOLLOWING STATEMENTS DOWN TO ' 1 CONTINUE' PRODUCE THE
C CROSS SPECTRA FOR ALL PAIRS OF SERIES
C SUBROUTINE WNDPLT PRODUCES A GRAPH OF THE FREQUENCY WINDOW
C OF THE POWER SPECTRUM ESTIMATES
IDDPL(2)=0
K=1
NVM1 = NV - 1
DO 1 I = 1, NVM1
IP1 = I + 1
DO 1 J = IP1, NV
II(K)=I
INAME(K)=NAMCHN(I)
JJ(K)=J
JNAME(K)=NAMCHN(J)
IF ( (K .EQ. 5) .OR. ( (J .EQ. NV) .AND. (I .EQ. NVM1))) GO TO 2
K=K+1
GO TO 1
2 WRITE(6,8) DF
8 FORMAT (1H1,20X,32HCROSS SPECTRA BETWEEN THE SERIES,10X,F9.2,1X,
1 18HDEGREES OF FREEDOM)
WRITE(6,111) (INAME(L),L=1,K)
111 FORMAT(/10X,5(8X,A8,'/',6X))
WRITE(6,112) (JNAME(L),L=1,K)
112 FORMAT(10X,5(11X,A8,4X))
WRITE(6,4) FREQ, (AMP,PHASE,COH ,L=1,K)
4 FORMAT ( 3X, A4, 2X, 5(5X, A3, 3X, A3, A2, 3X, A3, 1X) )
NDP=0
IDDPL(21)=0
DO 5 L=1,LMAX
FR=SR*FLOAT(L-1)/FLOAT(NB*2)
IF(LBUNIT.LE.0.OR.L.EQ.1.OR.L.GE.51)GO TO 215
NDP=NDP+1
DPDATA(NDP)=FR
IDDPL(21)=IDDPL(21)+1
215 CONTINUE
DO 6 M=1,K
NIJ=(II(M)-1)*NV+JJ(M)-((II(M)-1)*II(M))/2
NII=(II(M)-1)*NV+II(M)-((II(M)-1)*II(M))/2
NJJ=(JJ(M)-1)*NV+JJ(M)-((JJ(M)-1)*JJ(M))/2
C=P(2*L-1,NIJ)
Q=P(2*L,NIJ)
AM(M)=SQRT(C**2+Q**2)
IF ( (Q .EQ. 0.0) .AND. (C .EQ. 0.0) ) GO TO 25
PH(M) = 180.0 * ATAN2(Q, C) / PI
GO TO 26
25 PH(M) = 0.0
26 CONTINUE
PP = SQRT(P(2*L-1,NII)*P(2*L-1,NJJ))
IF ( (PP .EQ. 0.0) .AND. (AM(M) .EQ. 0.0) ) GO TO 27
CO(M) = AM(M) / PP
GO TO 28
27 CO(M) = 0.0
28 CONTINUE
AM(M)=AM(M)*CAL(JJ(M))*CAL(II(M))
IF(LBUNIT.LE.0.OR.L.EQ.1.OR.L.GE.51)GO TO 6
DPDATA(NDP+1)=AM(M)*2.0
DPDATA(NDP+2)=PH(M)
DPDATA(NDP+3)=CO(M)
NDP=NDP+3
6 CONTINUE
5 WRITE(6,7) FR, (AM(M),PH(M),CO(M) ,M=1,K)
7 FORMAT (1X, F8.3, 5(F10.3, 1X, F6.1, 1X, F5.3) )
IDDPL(20)=K
K=1
IF(LBUNIT.LE.0) GO TO 1
IDDPL(1)=4
IDDPL(2)=IDDPL(2)+1
C IDDPL(21)=LMAX-1
WRITE(LBUNIT)NID,IDDPL,NDP,(DPDATA(L),L=1,NDP)
1 CONTINUE
29 CONTINUE
IF (SPEC .NE. YES) GO TO 151
LLMAX = 2*(NB + 1)
NVP1 = NV + 1
NV2 = 2 * NV
DO 82 LL = 1, LLMAX, 2
DO 82 I = 1, NV
N = 0
IM1 = I - 1
K = I - NV
IF (I .EQ. 1) GO TO 85
DO 86 J = 1, IM1
K = K + NVP1 - J
N = N +1
COMSPC(N) = P(LL, K)
N = N + 1
86 COMSPC(N) = - P(LL + 1, K)
85 K = K + NVP1 - I
DO 80 J = I, NV
N = N + 1
COMSPC(N) = P(LL, K)
N = N + 1
COMSPC(N) = P(LL + 1, K)
80 K = K + 1
WRITE(LCT,122) (COMSPC(J),J=1,NV2)
122 FORMAT ( 30A4 )
82 CONTINUE
151 CONTINUE
RETURN
END
SUBROUTINE TRANS (P, Y, N, JSCANS, NV, NB, LOG2NS )
DIMENSION P(1), Y(JSCANS,1)
M = LOG2NS
DO 1 I=1,NV,2
1 CALL FAST(Y(1,I),Y(1,I+1),M)
NT=N/2
DO 8 I=1,NV,2
DO 8 J=2,NT
K=N-J+2
T1=(Y(J,I+1)-Y(K,I+1))
T2=(Y(K,I)-Y(J,I))
Y(J,I)=(Y(J,I)+Y(K,I))
Y(J,I+1)=(Y(J,I+1)+Y(K,I+1))
Y(K,I)=T1
8 Y(K,I+1)=T2
NE=N/(4*NB)
IIMIN=NE+1
IIMAX=NT-NE
IISTEP=2*NE
L=1
DO 2 I=1,NV
DO 2 J=I,NV
P(L)=Y(1,I)*Y(1,J)
P(L+1)=0.0
DO 3 K=2,IIMIN
M=N-K+2
P(L)=P(L)+2.0*(Y(K,I)*Y(K,J)+Y(M,I)*Y(M,J))
3 CONTINUE
L=L+2
DO 4 II=IIMIN,IIMAX,IISTEP
P(L)=0.0
P(L+1)=0.0
KMAX=II+IISTEP
DO 5 K=II,KMAX
M=N-K+2
P(L)=P(L)+Y(K,I)*Y(K,J)+Y(M,I)*Y(M,J)
5 P(L+1)=P(L+1)+Y(K,I)*Y(M,J)-Y(M,I)*Y(K,J)
4 L=L+2
P(L)=Y(NT+1,I)*Y(NT+1,J)
P(L+1)=0.0
KMIN=NT+1-NE
DO 6 K=KMIN,NT
M=N-K+2
P(L)=P(L)+2.0*(Y(K,I)*Y(K,J)+Y(M,I)*Y(M,J))
6 CONTINUE
2 L=L+2
RETURN
END
SUBROUTINE FAST(X,Y,M)
C SUBROUTINE FAST OBTAINS FINITE FOURIER TRANSFORMS OF THE
C SERIES IN PAIRS
DIMENSION X(1),Y(1)
N=2**M
DO 1 L=1,M
IMAX=2**(M-L)
JDELT=2*IMAX
ARG=6.2831853/FLOAT(JDELT)
C=COS(ARG)
S=SIN(ARG)
U=1.0
V=0.0
DO 1 I=1,IMAX
DO 2 J=I,N,JDELT
K=J+IMAX
XJ=X(J)+X(K)
YJ=Y(J)+Y(K)
XK=X(J)-X(K)
YK=Y(J)-Y(K)
X(K)=U*XK-V*YK
Y(K)=U*YK+V*XK
X(J)=XJ
2 Y(J)=YJ
UT=C*U-S*V
V=C*V+S*U
1 U=UT
J=1
NT=N/2
IMAX=N-1
DO 3 I=1,IMAX
IF(I.GE.J) GO TO 5
XT=X(J)
X(J)=X(I)
X(I)=XT
YT=Y(J)
Y(J)=Y(I)
Y(I)=YT
5 K=NT
4 IF(K.GE.J) GO TO 3
J=J-K
K=K/2
GO TO 4
3 J=J+K
RETURN
END
SUBROUTINE GENWTS(ID, L, WTS)
DIMENSION WTS(1)
DATA PI / 3.1415927 /
IF (L .EQ. 0) GO TO 21
FK = L / ID - 2
FKK = FK + 0.5
FL = L
DO 22 IS = 1, L
S = IS
INDEX1 = (L + 1) + IS
HS = (1.0 + COS(PI * S / FL) ) / (4.0 * FL)
22 WTS(INDEX1) = HS * SIN(PI * FKK * S / FL) / SIN(PI*S/ (2.0*FL) )
DO 23 IS = 1, L
INDEX1 = L + 1 - IS
INDEX2 = L + 1 + IS
23 WTS(INDEX1) = WTS(INDEX2)
WTS(L + 1) = FKK / FL
RETURN
21 WTS(1) = 1.0
RETURN
END
SUBROUTINE GENFLT ( H, NW, X, Y, L, Q)
DIMENSION H(2),Q(100),X(2),Y(2)
DATA P/6.2831853/
NW1=NW+1
L1=L+1
DO 1 I=2,L
1 Q(I)=(Y(I)-Y(I-1))/(X(I)-X(I-1))
Q(1)=0.
Q(L1)=0.
DO 2 I=2,NW1
T=P*FLOAT(I-1)
TT=T*T
H(I)=0.
DO 3 J=2,L1
3 H(I)=H(I)+(Q(J-1)-Q(J))*COS(T*X(J-1))/TT
2 H(I)=(H(I)+(Y(L)*SIN(T*X(L))-Y(1)*SIN(T*X(1)))/T)*2.
T=2.*(Y(L)*X(L)-Y(1)*X(1))
DO 4 J=2,L
4 T=T-(Y(J)-Y(J-1))*(X(J)+X(J-1))
H(1)=T
RETURN
END
SUBROUTINE PLOT1(AMPL, I, LMAX, SR, IFLAG)
DIMENSION AMPL(LMAX)
DOUBLE PRECISION NAMCHN
COMMON/CHANEL/NAMCHN(32),NVU,NCODE,NPULS,NNV(32)
AMPMAX = AMPL(1)
AMPMIN = AMPL(1)
DO 1 I1 = 2, LMAX
IF (AMPL(I1) .GT. AMPMAX) AMPMAX = AMPL(I1)
IF (AMPL(I1) .LT. AMPMIN) AMPMIN = AMPL(I1)
1 CONTINUE
IF (IFLAG .EQ. 2 .OR. IFLAG .EQ. 3) GO TO 2
CALL ROUND(AMPMAX, AMPLO, IAMP)
WIDTH = 0.01 * FLOAT(IAMP)
BASE = 0.0
IF(IFLAG.EQ.1) WRITE(6,110) NAMCHN(I)
110 FORMAT (1H1, 47HGRAPH OF THE POWER AT EACH FREQUENCY FOR SERIES,
X 1X,A4// 1H0,9HFREQUENCY,51X,5HPOWER)
IF (IFLAG .EQ. 4) WRITE (6, 124)
124 FORMAT (1H1, 61HGRAPH OF THE FREQUENCY WINDOW OF THE POWER SPECTRU
XM ESTIMATES / / 1H0, 9HFREQUENCY, 51X, 5HPOWER )
IF (IFLAG .EQ. 4) IFLAG = 1
CALL GRAPH(AMPL, LMAX, SR, BASE, WIDTH, AMPLO, IFLAG)
RETURN
2 CONTINUE
RANGE = 1.15 * (AMPMAX - AMPMIN)
CALL ROUND (RANGE, SIZELO, ISIZE)
FSIZE = ISIZE
UNITS = .1 * FSIZE * SIZELO
ISTEPS = AMPMIN / UNITS
FSTEPS = ISTEPS
TBASE = FSTEPS * UNITS
13 CONTINUE
IF (TBASE .LT. AMPMIN) GO TO 12
FSTEPS = FSTEPS - 1.0
TBASE = FSTEPS * UNITS
GO TO 13
12 CONTINUE
BASE = FSTEPS * 0.1 * FSIZE
WIDTH = 0.01 * FSIZE
17 CONTINUE
IF (ABS(BASE) .LE. 10.5) GO TO 18
BASE = .1 * BASE
WIDTH = .1 * WIDTH
SIZELO = 10. * SIZELO
GO TO 17
18 CONTINUE
120 FORMAT (1H1, 64HGRAPH OF THE LOGARITHM OF THE POWER AT EACH FREQUE
XNCY FOR SERIES,1X,A4//1H0,'FREQUENCY',41X,'LOGARITHM OF THE POWER'
X)
IF (IFLAG .EQ. 2) WRITE (6,120) NAMCHN(I)
IF (IFLAG .EQ. 3) WRITE (6, 123)
123 FORMAT (1H1, 57HGRAPH OF THE FREQUENCY RESPONSE FUNCTION OF THE PR
XEFILTER / / 1H0, 9HFREQUENCY, 51X, 5HPOWER )
IF (IFLAG .EQ. 3) IFLAG = 1
CALL GRAPH(AMPL, LMAX, SR, BASE, WIDTH, SIZELO, IFLAG)
RETURN
END
SUBROUTINE ROUND(RANGE, SIZELO, ISIZE)
SIZEHI = 10.0
SIZELO = 1.0
21 CONTINUE
IF ( (SIZELO .LT. RANGE) .AND. (RANGE .LE. SIZEHI) ) GO TO 22
IF (RANGE .GE. SIZEHI) GO TO 23
SIZEHI = .1 * SIZEHI
SIZELO = .1 * SIZELO
GO TO 21
23 SIZEHI = 10. * SIZEHI
SIZELO = 10. * SIZELO
GO TO 21
22 CONTINUE
DO 25 I1 = 2, 10
FI1 = I1
SIZE = FI1 * SIZELO
ISIZE = I1
IF (RANGE .LT. SIZE) GO TO 24
25 CONTINUE
24 CONTINUE
RETURN
END
SUBROUTINE GRAPH(AMPL, LMAX, SR, BASE, WIDTH, SIZE, IFLAG)
DIMENSION AMPL(LMAX)
DIMENSION ARRAY(100), SCALE(10)
DIMENSION TITLE(11, 2)
DATA TITLE /
X 4H , 4H , 4H , 4H , 4H MU, 4HLTIP, 4HLY E, 4HACH ,
X 4HPOWE, 4HR BY, 4H ,
X 4H MUL, 4HTIPL, 4HY EA, 4HCH L, 4HOGAR, 4HITHM, 4H OF , 4HTHE ,
X 4HPOWE, 4HR BY, 4H /
DATA BLANK, STAR / 1H , 1H* /
NB = LMAX - 1
DO 31 I1 = 1, 100
ARRAY(I1) = BLANK
31 CONTINUE
ARRAY0 = BLANK
TBASE = BASE * SIZE
TWIDTH = WIDTH * SIZE
IF (SIZE .GT. 500000.) GO TO 25
IF (SIZE .GT. .5) GO TO 26
IF (SIZE .GT. .00005) GO TO 27
25 WRITE (6, 111) (TITLE(I2, IFLAG), I2= 1, 11), SIZE
111 FORMAT (1H0, 45X, 11A4, 1X, 1PE8.1 )
GO TO 29
26 WRITE (6, 112) (TITLE(I2, IFLAG), I2= 1, 11), SIZE
112 FORMAT (1H0, 45X, 11A4, 1X, F7.0 )
GO TO 29
27 WRITE (6, 113) (TITLE(I2, IFLAG), I2= 1, 11), SIZE
113 FORMAT (1H0, 45X, 11A4, 1X, F7.4 )
29 CONTINUE
DO 24 I1 = 1, 10
SCALE(I1) = BASE + 10. * FLOAT(I1) * WIDTH
24 CONTINUE
WRITE (6, 114) BASE, SCALE
114 FORMAT (1H0, 3X, 11F10.4 )
WRITE (6, 115)
115 FORMAT(1H0,9X,10(10H+.........),1H+)
WRITE (6, 121)
121 FORMAT (1X, 9X, 1H., 99X, 1H. / 1X, 9X, 1H., 99X, 1H. )
WRITE (6, 122)
122 FORMAT (1X)
DO 32 I1 = 1, LMAX
FR = SR * FLOAT(I1 - 1) / FLOAT(NB * 2)
INDEX = (AMPL(I1) - TBASE) / TWIDTH + 0.5
IF ( (INDEX .GE. 1) .AND. (INDEX .LE. 100) ) ARRAY(INDEX) = STAR
IF (INDEX .EQ. 0) ARRAY0 = STAR
WRITE (6, 119) FR, ARRAY0, ARRAY, AMPL(I1)
119 FORMAT (1X, F8.3, 1X, A1, 100A1, 2X, 1PE10.2 )
IF ( (INDEX .GE. 1) .AND. (INDEX .LE. 100) ) ARRAY(INDEX) =BLANK
ARRAY0 = BLANK
32 CONTINUE
WRITE (6, 122)
WRITE (6, 121)
WRITE (6, 115)
WRITE (6, 122)
WRITE (6, 114) BASE, SCALE
IF (SIZE .GT. 500000. ) GO TO 33
IF (SIZE .GT. .5 ) GO TO 34
IF (SIZE .GT. .00005 ) GO TO 35
33 WRITE (6, 111) (TITLE(I2, IFLAG), I2= 1, 11), SIZE
GO TO 36
34 WRITE (6, 112) (TITLE(I2, IFLAG), I2= 1, 11), SIZE
GO TO 36
35 WRITE (6, 113) (TITLE(I2, IFLAG), I2= 1, 11), SIZE
36 CONTINUE
RETURN
END
SUBROUTINE WNDPLT ( NSCARS, WINDOW )
DIMENSION WINDOW ( NSCARS, 2)
COMMON /LIST/ NSC, LEN4, MSCANS, NV, NCHANS, NB, SR, FILTER, ID,
1 L, PLTFLT, SPEC, PLTWND, IT, DPLTAP, IBEGIN, NF, NWORDS, NWTS,
2 SAMPRT, LOG2NS, PLTPOW, PLTLOG, NS, NSCANS,NFF,LFT,LCT
DATA PI / 3.1425927 /
C SET ALL ELEMENTS OF SECOND COLUMN OF WINDOW EQUAL TO 0
DO 403 I1 = 1, NSCANS
403 WINDOW(I1, 2) = 0.0
C SET UP WINDOW WEIGHTS IN FIRST COLUMN OF WINDOW
IR = NS / 10
R = IR
DO 405 I1 = 1, IR
I3 = NS + 1 - I1
FI1 = I1
FINT = FI1 - 0.5
WINDOW(I1, 1) = 0.5 * (1.0 - COS(PI * FINT / R) )
WINDOW(I3, 1) = WINDOW(I1, 1)
405 CONTINUE
IF (NS .EQ. NSCANS) GO TO 406
I1BEGN = NS + 1
DO 407 I1 = I1BEGN, NSCANS
407 WINDOW(I1, 1) = 0.0
406 CONTINUE
IRP1 = IR + 1
I1END = NS - IR
DO 408 I1 = IRP1, I1END
408 WINDOW(I1, 1) = 1.0
C SUBROUTINE FAST OBTAINS THE FINITE FOURIER TRANSFORM OF THE
C WEIGHTS USED FOR WINDOWING
CALL FAST(WINDOW(1, 1), WINDOW(1, 2), LOG2NS)
C REPLACE EACH ELEMENT OF FIRST COLUMN OF WINDOW WITH SQUARE
C OF AMPLITUDE OF CORRESPONDING ELEMENT OF TRANSFORMED SERIES
DO 409 I1 = 1, NSCANS
409 WINDOW(I1, 1) = WINDOW(I1, 1)**2 + WINDOW(I1, 2)**2
C ZERO OUT COLUMN 2 OF WINDOW
DO 410 I1 = 1, NSCANS
410 WINDOW(I1, 2) = 0.0
ID = NSCANS / (2 * NB)
I3D = 3 * ID
NSD2 = NSCANS / 2
LMAX = MIN0(I3D, NSD2)
C FORM SUMS OF ELEMENTS OF COLUMN 1 OF WINDOW AND PUT THEM IN
C IN COLUMN 2 OF WINDOW
IID = NSCANS / (2 * NB)
IHALFD = IID / 2
IIDP1 = IID + 1
I5 = NSCANS - IHALFD + 1
DO 411 I1 = 1, LMAX
I4 = I5
DO 412 I3 = 1, IIDP1
WINDOW(I1, 2) = WINDOW(I1, 2) + WINDOW(I4, 1)
I4 = I4 + 1
IF (I4 .EQ. NSCANS + 1) I4 = 1
412 CONTINUE
I5 = I5 + 1
IF (I5 .EQ. NSCANS + 1) I5 = 1
411 CONTINUE
C SCALE ELEMENTS OF COLUMN 2 BY MULTIPLYING BY C
FNS = NS
WNDPWR = FNS - 1.25 * R
FSCANS = NSCANS
FNB = NB
FD = FSCANS / (2.0 * FNB)
C = 1.0 / (SR * (FD + 1.0) * WNDPWR)
DO 413 I1 = 1, LMAX
413 WINDOW(I1 , 2) = C * WINDOW(I1, 2)
SR1 = SR * FLOAT(LMAX) / FLOAT(NSD2)
C SUBROUTINE PLOT1 GRAPHS THE ELEMENTS OF COLUMN 2 OF WINDOW
C TO PRODUCE THE FREQUENCY WINDOW OF THE POWER SPECTRUM
C ESTIMATES. IFLAG = 2 PRODUCES THE APPROPRIATE TITLE ON
C THE GRAPH AND PRODUCES A BASE LINE OF ZERO
IFLAG = 4
CALL PLOT1(WINDOW(1, 2), I, LMAX, SR1, IFLAG)
RETURN
END
SUBROUTINE POSIT
C
WRITE (6,10)
10 FORMAT ('0INPUT TAPE WAS SPECIFIED AS UCLA BRI/DPL DIGITIZED TAPE.
2'/' FOR INFORMATION ON THIS KIND OF INPUT, PLEASE WRITE TO THE UCL
3A HEALTH SCIENCES COMPUTING FACILITY.'
4 ' THE PROGRAM HAS BEEN TERMINATED.')
STOP
END
SUBROUTINE BEANS
DOUBLE PRECISION NAMCHN,LBTAPE
COMMON /FRITOS/ NTOT,JSCANS,ISTART,IWTS,ITEMP,IRR,LOGIRR,IFT2,
X IROWSP,IAMPL,ILOGA,ISPEC,IBRI,NWO,ISIZE2
COMMON/CHANEL/NAMCHN(32),NVU,NCODE,NPULS,NNV(32)
COMMON/BRIID/LBTAPE,LBUNIT,IDDPL(29),NID,TITLER(68),ID1BRI,ID2BRI
COMMON /LIST/ NSC, LEN4, MSCANS, NV, NCHANS, NB, SR, FILTER, ID,
1 L, PLTFLT, SPEC, PLTWND, IT, DPLTAP, IBEGIN, NF, NWORDS, NWTS,
2 SAMPRT, LOG2NS, PLTPOW, PLTLOG, NS, NSCANS,NFF,LFT,LCT
DATA YES/'YES '/
C SAMPRT IS THE INITIAL SAMPLING RATE (LOWER CASE S' IN THE
C WRITE UP). SR IS THE SAMPLING RATE AFTER PREFILTERING.
C SR = SAMPRT IF NO PREFILTERING IS DONE. (SR IS LOWER
C CASE S IN THE WRITE UP)
SR = SAMPRT / FLOAT(ID)
C THE ARRAY Q IS USED FOR SEVERAL PURPOSES IN THE PROGRAM.
C THESE ARE
C 1) AN ARRAY, CALLED WTS IN SUBROUTINE INPUT, USED FOR
C STORING PREFILTER WEIGHTS
C 2) AN ARRAY, CALLED FMT IN SUBROUTINE INPUT, USED FOR
C STORING VARIABLE FORMAT CARDS
C 3) AN ARRAY OF LENGTH NV, CALLED TEMP IN SUBROUTINE
C INPUT, USED FOR TEMPORARY STORAGE OF INPUT DATA
C 4) AN ARRAY, CALLED SCANS IN SUBROUTINE INPUT AND
C Y IN SUBROUTINE TRANS, WHICH IS USED FOR STORING THE
C SERIES BEING ANALYZED
C 5) AN ARRAY, CALLED FILT IN SUBROUTINE INPUT, USED FOR
C PRODUCING THE GRAPH OF THE FREQUENCY RESPONSE
C FUNCTION OF THE LOW PASS FILTER
C 6) AN ARRAY, CALLED FILT IN SUBROUTINE INPUT AND
C WINDOW IN SUBROUTINE WNDPLT, USED FOR PRODUCING THE
C GRAPH OF THE FREQUENCY WINDOW OF THE POWER SPECTRUM
C ESTIMATES
C 7) AN ARRAY, CALLED P IN SUBROUTINES INPUT, TRANS, AND
C OUTPUT, USED FOR STORING THE MATRIX OF CROSS
C SPECTRUM ESTIMATES
C 8) TWO ARRAYS, CALLED AMPL AND ALOGAM IN SUBROUTINE
C OUTPUT, USED IN PRODUCING THE GRAPHS OF POWER AND
C LOGARITHM OF POWER FOR THE SERIES
C 9) AN ARRAY, CALLED COMSPC IN SUBROUTINE OUTPUT, USED
C FOR WRITING THE MATRICES OF CROSS SPECTRUM ESTIMATES
C ON FORTRAN UNIT 1.
C THESE VARIOUS ARRAYS ARE OVERLAPPED IN THE FOLLOWING
C MANNER. FIRST ITEMS (1), (2), (3), AND (6) OCCUPY Q. THEN
C ITEM (5) REPLACES ITEM (6). THEN ITEM (4) REPLACES ITEM
C (5). THEN ITEM (4) IS CONSUMED IN PIECES TO PRODUCE ITEM
C (7) WHICH BEGINS AT THE START OF Q AND REPLACES (1), (2),
C AND (3). AS PARTS OF ITEM (4) ARE USED ITEM (7) REPLACES
C THEM. ITEM (8) FOLLOWS ITEM (7) IN Q. THEN ITEM (9)
C REPLACES ITEM (8) IN Q.
C THE SUBROUTINE THAT DOES FINITE FOURIER TRANSFORMS DOES 2
C SERIES AT A TIME. IF THERE IS AN ODD NUMBER OF SERIES, IT
C IS NECESSARY TO PROVIDE A DUMMY SERIES FILLED WITH ZEROS
C TO BE USED IN TRANSFORMING THE LAST SERIES. NVPRIM IS THE
C SMALLEST EVEN NUMBER GREATER THAN OR EQUAL TO NV.
NCHANS=((NVU+1)/2)*2
IF(NSC.GT.0) GO TO 699
LEN1 = 2 * NVU * ( NB + 1 )
JSCAN1=(NTOT-LEN1)/NCHANS
JSCANS=1
54 JSCANS=JSCANS*2
IF(JSCANS.LT.JSCAN1)GO TO 54
55 JSCANS=JSCANS/2
LEN2=JSCANS*NCHANS
5501 MSCANS=JSCANS*ID
5502 LEN1=NTOT-LEN2
NS=JSCANS
NSCANS=JSCANS
GO TO 5503
699 JSCANS=NSCANS
MSCANS=NSC
IF(DPLTAP.EQ.YES) MSCANS=NSC+1
LEN2=JSCANS*NCHANS
C ADD 1 SCAN TO MSCANS SO THAT READ WILL STOPPED BY PICKUP
C IF NUMBER OF SCANS KNOWN,JSCANS=NO OF EFF SCANS (NEAREST POWER OF
C AND MSCANS = TOTAL SCANS READ.
C IF NUMBER OF SCANS UNKNOWN,JSCANS=MAX THAT CAN BE STORED
C AND MSCANS = MAX THAT CAN BE READ.
LEN1 = 0
IF ((2*NVU-1)*(NB+1) .LT. NSCANS) GO TO 700
LEN1=((2*NVU-1)*(NB+1)-NSCANS)**2
LEN1=LEN1/(4*(NB+1))
700 LEN1 = LEN1 + 2 * NVU * (NB+1) + 1
C
C NWORDS IS THE AMOUNT OF SPACE THAT IS USED TO HOLD THE
C VARIABLE FORMAT CARDS
5503 CONTINUE
C IWTS IS WHERE THE ARRAY WTS BEGINS IN Q
IFT2=NWORDS+1
IWTS=NWORDS+NWO+1
C NWTS IS THE AMOUNT OF SPACE USED TO HOLD THE PREFILTER
C WEIGHTS
NWTS = 2 * L + 1
C ITEMP IS WHERE THE ARRAY TEMP BEGINS IN Q
ITEMP = IWTS + NWTS
C LEN3 IS THE AMOUNT OF SPACE OCCUPIED BY ITEMS (1), (2), AND
C (3)
LEN3=NWORDS+NWO+NVU+NWTS
C WHILE ITEM (4) IS IN Q THERE MUST BE SPACE FOR ITEMS (1),
C (2), AND (3) AND LATER FOR THE FIRST PART OF ITEM (7).
C THEREFORE LEN4, THE MAXIMUM OF LEN1 AND LEN3, IS FOUND.
LEN4 = MAX0(LEN1, LEN3)
C ISTART IS WHERE ITEMS (4), (5), AND (6) START IN Q
ISTART=LEN4
C IRR IS THE SMALLEST POWER OF 2 GREATER THAN OR EQUAL TO 4*L
C 2**LOGIRR = IRR
IRR = 1
IF (PLTFLT .NE. YES) GO TO 707
LENSER = 4 * L
LOGIRR = 0
702 IRR=2*IRR
LOGIRR = LOGIRR + 1
IF(LENSER.GT.IRR)GO TO 702
707 CONTINUE
IROWSP = 2 * (NB + 1)
ICOLSP=(NVU*(NVU+1))/2
ISIZEP = IROWSP * ICOLSP
IAMPL = ISIZEP + 1
ILOGA = ISIZEP + 1
IF (PLTPOW .EQ. YES) ILOGA = ILOGA + (NB + 1)
IF(LBUNIT.GT.0) IBRI = ILOGA + (NB +1)
ISPEC = ISIZEP + 1
LEN5 = 0
LEN6 = 0
IF (PLTPOW .EQ. YES) LEN5 = LEN5 + (NB + 1)
IF (PLTLOG .EQ. YES) LEN5 = LEN5 + (NB + 1)
IF(LBUNIT.GT.0) LEN5 = LEN5 + MAX0(456,16*MAX0(35,NB))
IF (SPEC .EQ. YES) LEN6 = 2 * NV
ISIZE2 = ISIZEP + MAX0(LEN5, LEN6)
RETURN
END