C	SPECTRA3.FOR - CROSS PERIODOGRAM SMOOTHING FROM TWO SERIES
C	TO ESTIMATE COHERENCY AND PHASE SPECTRUM.  COHERENCY
C	GIVES SHARED HARMONICS AT SAME FREQUENCY AND PHASE SPECTRUM
C	GIVES THE DIRECTION OF THIS ASSOCIATION.  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 TRANSFORMS OF BOTH SERIES TO BE SMOOTHED AND THE 
C	SPECTRUMS OF EACH ARE READ IN BY DATIN.  THE SMOOTHING OF
C	BOTH SPECTRUMS (K,NK) MUST HAVE BEEN SAME FOR BOTH SERIES
C	AND IDENTICAL TO THAT SPECIFIED HERE.
C
	DIMENSION TR1(513),TI1(513),TR2(513),TI2(513),K(10)
	CHARACTER * 40 DF
	DATA PI /3.141593/
	WRITE(1,2050)
 2050	FORMAT(/" SPECTRA3 - Cross-Spectrum by Smoothed Cross-Periodogram")
 2001	WRITE(1) "Enter UPPERCASE Transform Input File Name for 1st File: "
	READ(1) DF
	IF(IOREAD(5,2,0,DF)) GOTO 2001
	CALL DATIN(TR1,NPGM1,START,STEP,5)
	CALL DATIN(TI1,NPGM1,START,STEP,5)
	I=IOCLOS(5)
 2002	WRITE(1) "Enter UPPERCASE Transform Input File Name for 2nd File: "
	READ(1) DF
	IF(IOREAD(5,2,0,DF)) GOTO 2002
	CALL DATIN(TR2,NPGM2,START,STEP,5)
	CALL DATIN(TI2,NPGM2,START,STEP,5)
	I=IOCLOS(5)
	IF(NPGM1.EQ.NPGM2) GOTO 20
	WRITE(1,2) NPGM1,NPGM2
    2	FORMAT(' *** ERROR - TRANSFORM LENGTHS ',2I5,' DIFFER.')
	STOP
   20	N=NPGM1
	NP2=(N-1)*2
	WRITE(1,7001)
 7001	FORMAT(' Number of Observations in Original Time Series: '$
        READ(1,7002) NOBS
 7002	FORMAT(I0)
	CON=FLOAT(NP2)**2/(PI*FLOAT(NOBS)*2.0)
	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)
	DO 30 I=1,N
	CR=TR1(I)*TR2(I)+TI1(I)*TI2(I)
	CI=TI1(I)*TR2(I)-TR1(I)*TI2(I)
	TR1(I)=CR*CON
   30	TI1(I)=CI*CON
	WRITE(1,7006)
 7006	FORMAT(' Smoothing Pass - '$
	DO 40 I=1,NK
	WRITE(1,7007) I
 7007	FORMAT(' ',I2$
	CALL MODDAN(TR1,TR2,N,K(I),1.0)
   40	CALL MODDAN(TI1,TI2,N,K(I),-1.0)
	WRITE(1,7008)
 7008	FORMAT(1X)
	CALL POLAR(TR1,TI1,N)
 2003	WRITE(1) "Enter UPPERCASE Spectrum Input File Name for 1st File: "
	READ(1) DF
	IF(IOREAD(5,2,0,DF)) GOTO 2003
	CALL DATIN(TR2,N1,START,STEP,5)
	I=IOCLOS(5)
 2004	WRITE(1) "Enter UPPERCASE Spectrum Input File Name for 2nd File: "
	READ(1) DF
	IF(IOREAD(5,2,0,DF)) GOTO 2003
	CALL DATIN(TI2,N2,START,STEP,5)
	I=IOCLOS(5)
	IF((N1.EQ.N).AND.(N2.EQ.N)) GOTO 60
	WRITE(1,3) N,N1,N2
    3	FORMAT(' *** ERROR - TRANSFORM LENGTH ',I5,' NOT EQUAL TO '
     *  'SPECTRUM LENGTHS ',2I5)
	STOP
   60	DO 50 I=1,N
   50	TR1(I)=TR1(I)/SQRT(TR2(I)*TI2(I))
	START=0.0
	STEP=PI/FLOAT(N-1)
 2052	WRITE(1) "Enter UPPERCASE File for Cross-Spectrum Output: "
	READ(1) DF
	IF(IOWRIT(8,2,0,DF)) GOTO 2052
	DF="CROSS-SPECTRUM ESTIMATE"
	CALL DATOUC(TR1,TI1,N,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 POLAR(X,Y,N)
C
C	CONVERTS REAL AND IMAGINARY PARTS OF CROSS SPECTRUM TO
C	MAGNITUDES AND PHASES.  CONVERSION IS DONE IN PLACE GIVEN:
C
C	X	-  REAL ARRAY  -  REAL PART	-  RETURNS MAGNITUDE
C	Y	-  REAL ARRAY  -  IMAG PART	-  RETURNS PHASE
C	N	-  INTEGER     -  SERIES LENGTH	-  (UNCHANGED)
C
	DIMENSION X(N),Y(N)
	DO 10 I=1,N
	R=SQRT(X(I)**2+Y(I)**2)
	PHI=ATAN2(Y(I),X(I))
	X(I)=R
   10	Y(I)=PHI
	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 DATOUC(X,Y,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, (3) COHERENCY
C	(4) SQUARED COHERENCY, (5) PHASE SPECTRUM AND (6) TIME LEAD
C       AT FREQUENCEY ARE OUTPUT IN A FORMAT SO DATA MAY BE
C	EDITED FOR PLOTTING SOFTWARE INPUT.
C
C	X	-  REAL ARRAY  -  COHERENCY
C	X	-  REAL ARRAY  -  PHASE SPECTRUM
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),Y(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/' *CANNOT BE REPROCESSED*'/1X,2F10.5)
	Z=1024.0
	DO 3 I=1,N
        X2=X(I)**2
        T=0.0
        IF(START.GT.0.0) T=Y(I)/START
	WRITE(M,2) START,Z,X(I),X2,Y(I),T
    2	FORMAT(1X,F7.5,1X,F9.4,1X,F7.5,1X,F7.5,1X,F9.5,1X,F10.4)
	START=START+STEP
    3	Z=2.0*PI/START
	RETURN
	END