C	SPECTRA2.FOR - PERIODOGRAM SMOOTHING TO OBTAIN SPECTRUM
C	ESTIMATE FROM OUTPUT OF SPECTRAL.FOR.  SPECTRUM ESTIMATE
C	IS COMPUTED BY REPEATED SMOOTHING WITH MODIFIED DANIELL
C	WEIGHTS.  INPUT IS:
C
C	NOBS	-  INTEGER	 - LENGTH OF ORIGINAL TIME SERIES
C	NK	-  INTEGER	 - NUMBER OF SMOOTHING PASSES
C	K	-  INTEGER ARRAY - ARRAY OF SMOOTHING PARAMETERS
C
C	THE PERIODOGRAM TO BE SMOOTHED IS READ IN BY DATIN.
C
	DIMENSION X(513),Y(513),K(10)
	CHARACTER * 40 DF
	DATA PI /3.141593/
	WRITE(1,2050)
 2050	FORMAT(/" SPECTRA2 - Spectrum Estimates by Smoothed Periodogram")
 2001	WRITE(1) "Enter UPPERCASE Periodogram Input File Name: "
	READ(1) DF
	IF(IOREAD(5,2,0,DF)) GOTO 2001
	CALL DATIN(X,NPGM,START,STEP,5)
	NP2=(NPGM-1)*2
	WRITE(1,7001)
 7001	FORMAT(' Number of Observations in Original Time Series: '$
        READ(1,7002) NOBS
 7002	FORMAT(I0)
	WRITE(1,7003)
 7003	FORMAT(' Number of Modified Daniell Passes to Make: '$
        READ(1,7002) NK
	DO 7005 I=1,NK
	WRITE(1,7004) I
 7004	FORMAT('    Weight (1/2 length of moving average) for pass ',I2,': '$
 7005   READ(1,7002) K(I)
	CON=FLOAT(NOBS)/(8.0*PI)
	DO 10 I=1,NPGM
   10	X(I)=X(I)*CON
	WRITE(1,7006)
 7006	FORMAT(' Smoothing Pass - '$
	DO 50 I=1,NK
	WRITE(1,7007) I
 7007	FORMAT(' ',I2$
   50	CALL MODDAN(X,Y,NPGM,K(I),1.0)
	WRITE(1,7008)
 7008	FORMAT(1X)
	START=0.0
	STEP=2.0*PI/FLOAT(NP2)
 2052	WRITE(1) "Enter UPPERCASE File for Spectrum Output: "
	READ(1) DF
	IF(IOWRIT(8,2,0,DF)) GOTO 2052
	DF="SPECTRUM ESTIMATE"
	CALL DATOUP(X,NPGM,START,STEP,8,DF)
	I=IOCLOS(8)
	STOP
	END
C
	FUNCTION EXTEND(X,I,N,SYM)
C
C	RETURNS THE ITH TERM IN SERIES X, EXTENDED IF NECESSARY WITH
C	EVEN OR ODD SYMMETRY ACCORDING TO SIGN OF SYM (+1 OR -1).
C	IF SYM=0 THE EXTENDED VALUE IS ZERO.
C
	DIMENSION X(N)
	IF(N.GT.1) GOTO 10
	WRITE(1,1) N
    1	FORMAT(' *** ERROR VALUE OF N IN EXTEND IS ',I5)
	STOP
   10	J=I
	CON=1
   20	IF(J.GE.1) GOTO 30
	J=2-J
	CON=CON*SYM
   30	IF(J.LE.N) GOTO 40
	J=2*N-J
	CON=CON*SYM
	GOTO 20
   40	EXTEND=X(J)*CON
	RETURN
	END
C
	SUBROUTINE MODDAN(X,Y,N,K,SYM)
C
C	MODIFIED DANIELL SMOOTHING TO SERIES X.  COMPUTES SPECTRUM
C	ESTIMATE FROM PERIODOGRAMS.  ASSUMES SERIES IS SYMMETRIC ABOUT
C	BOTH END VALUES TO PROVIDE END VALUES OF OUTPUT.  PARAMETERS
C	ARE:
C
C	X	-  REAL ARRAY  -  THE SERIES	-  RETURNS SMOOTHED
C	Y	-  REAL ARRAY  -  SCRATCH	-  RETURNS SERIES
C	K	-  INTEGER     -  HALF-LENGTH	-  (UNCHANGED)
C				  OF MOVING AVERAGE
C	SYM	-  REAL        -  +1 OR -1	-  (UNCHANGED)
C				  SYMMETRY INDICATOR
C	N	-  INTEGER     -  SERIES LENGTH	-  (UNCHANGED)
C
	DIMENSION X(N),Y(N)
	DO 10 I=1,N
   10	Y(I)=X(I)
	IF(K.LE.0) RETURN
	LIM=K-1
	CON=1.0/FLOAT(2*K)
	DO 20 I=1,N
	X(I)=Y(I)
	IF(LIM.EQ.0) GOTO 20
	DO 30 J=1,LIM
   30	X(I)=X(I)+EXTEND(Y,I-J,N,SYM)+EXTEND(Y,I+J,N,SYM)
   20	X(I)=(X(I)+(EXTEND(Y,I-K,N,SYM)+EXTEND(Y,I+K,N,SYM))*0.5)*CON
	RETURN
	END
C
	SUBROUTINE DATIN(X,N,START,STEP,M)
C
C	INPUT A SERIES OF VALUES IN RUNTIME FORMAT.  THE FIRST FOUR
C	CARDS MUST BE (FREE FORMAT):
C		1) A TITLE CARD (<72 COLS.)
C		2) NUMBER OF CASES (I5)
C		3) DATA FORMAT (<72 COLS.)
C		4) START TIME VALUE AND STEP (2F10.5)
C
C	X	-  REAL ARRAY  -  RETURN THE SERIES
C	N	-  INTEGER     -  RETURN SERIES LENGTH
C	START	-  REAL	       -  RETURN SERIES TIME VALUE AT 1ST POINT
C	STEP	-  REAL        -  RETURN TIME INCREMENT BETWEEN POINTS
C	M	-  INTEGER     -  UNIT NUMBER (UNCHANGED)
C
	DIMENSION X(513)
	CHARACTER * 72 HEAD,FMT
	WRITE(1,7001)
 7001	FORMAT(' Reading data...')
	READ(M,2) HEAD,N,FMT,START,STEP
    2	FORMAT(A0/I5/A0/2F10.5)
	WRITE(1,3) HEAD,N,FMT,START,STEP
    3	FORMAT(5X,'The Input Data Header Reads:'/5X,A0/
     *	5X,'Series Length Input is = ',I5/5X,'Format Of Input Data:'/
     *	5X,A0/5X,'Time Origin = ',F11.5,'  Time Increment = ',F11.5)
        READ(M,FMT) (X(I),I=1,N)
	RETURN
	END
C
	SUBROUTINE DATOUP(X,N,START,STEP,M,HEAD)
C
C	OUTPUT A SERIES OF VALUES IN FORMAT REQUIRED FOR DATIN.
C	UNLIKE DATOUT (1) FREQ VALUES, (2) PERIOD FREQ,  AND
C	(3) LOG OF SPECTRUM ARE OUTPUT WITH (4) SPECTRUM
C	SO DATA MAY BE EDITED FOR PLOTTING SOFTWARE INPUT.
C
C	X	-  REAL ARRAY  -  THE SERIES
C	N	-  INTEGER     -  SERIES LENGTH
C	START	-  REAL	       -  SERIES TIME VALUE AT 1ST POINT
C	STEP	-  REAL        -  TIME INCREMENT BETWEEN POINTS
C	M	-  INTEGER     -  UNIT NUMBER (UNCHANGED)
C	HEAD    -  CHAR*40     -  HEADER TITLE
C
	DIMENSION X(N)
	CHARACTER * 40 HEAD
	DATA PI /3.141593/
	WRITE(1,7001)
 7001	FORMAT(' Writing data...')
	WRITE(M,1) HEAD,N,START,STEP
    1	FORMAT(1X,A0/1X,I5/' (28X,E15.8)'/1X,2F10.5)
	Y=1024.0
	DO 3 I=1,N
	Z=ALOG(X(I))
	WRITE(M,2) START,Y,Z,X(I)
    2	FORMAT(1X,F7.5,1X,F9.4,1X,F9.5,1X,E15.8)
	START=START+STEP
    3	Y=2.0*PI/START
	RETURN
	END