C	SPECTRAL.FOR - COMPUTATION OF THE DISCRETE FOURIER TRANSFORM
C	AND THE PERIODOGRAM OF A TIME SERIES.  THE PERIODOGRAM IS
C	COMPUTED AS THE SQUARED AMPLITUDE R(W)**2 RATHER THAN THE
C	ACTUAL PERIODOGRAM (N/(8*PI))R(W)**2 FOR COMPUTATIONAL USE
C	IN SUBSEQUENT PROCESSING.
C
	DIMENSION X(1024),Y(1024)
	CHARACTER * 40 DF
	WRITE(1,2050)
 2050	FORMAT(/" SPECTRAL - Fourier Transform And Periodogram")
 2001	WRITE(1) "Enter UPPERCASE Input File Name: "
	READ(1) DF
	IF(IOREAD(5,2,0,DF)) GOTO 2001
	M=5
	DO 10 I=1,1024
	X(I)=0.0
   10	Y(I)=0.0
	CALL DATIN(X,N,START,STEP,M)
	WRITE(1,7001) 
 7001	FORMAT(' Order of Trend to Remove (0=Mean, 1=Linear): '$
        READ(1,7002) K
 7002	FORMAT(I0)
	CALL DETRND(X,N,K)
	WRITE(1,7003)
 7003	FORMAT(' Proportion of Data to Taper: '$
        READ(1,7004) P
 7004	FORMAT(F0.0)
	CALL TAPER(X,N,P)
	NP2=1
   20	NP2=NP2*2
	IF(NP2.LT.N) GOTO 20
	INV=0
	CALL FT01A(X,Y,NP2,INV)
	IF(INV.EQ.0) GOTO 30
	WRITE(1) "*** ERROR IN SUBROUTINE FT01A ***"
	GOTO 99
   30	NPGM=NP2/2+1
 2052	WRITE(1) "Enter UPPERCASE File for Transform Output: "
	READ(1) DF
	IF(IOWRIT(8,2,0,DF)) GOTO 2052
	DF="REAL PART OF TRANSFORM"
	CALL DATOUT(X,NPGM,0.0,1.0,8,DF)
	DF="IMAGINARY PART OF TRANSFORM"
	CALL DATOUT(Y,NPGM,0.0,1.0,8,DF)
        INV=IOCLOS(8)
	CON=(2.0*FLOAT(NP2)/FLOAT(N))**2
	DO 40 I=1,NPGM
   40	X(I)=(X(I)**2+Y(I)**2)*CON
 2053	WRITE(1) "Enter UPPERCASE File for Periodogram Output: "
	READ(1) DF
	IF(IOWRIT(8,2,0,DF)) GOTO 2053
	DF="SQUARED AMPLITUDE"
	CALL DATOUT(X,NPGM,0.0,1.0,8,DF)
        INV=IOCLOS(8)
   99	STOP
	END
C
	SUBROUTINE FT01A(XR,XI,N,INV)
C
C	THIS SUBROUTINE IMPLEMENTS THE SANDE-TUKEY RADIX-2 FAST FOURIER
C	TRANSFORM.  EITHER THE DIRECT OR INVERSE TRANSFORM MAY BE COMP-
C	UTED.  PARAMETERS ARE:
C
C	XR  -  REAL ARRAY  -  REAL PART OF THE  -  RETURNS REAL PART
C			      SERIES		   OF TRANSFORM
C	XI  -  REAL ARRAY  -  IMAGINARY PART OF -  RETURNS IMAGINARY
C			      THE SERIES	   PART OF TRANSFORM
C	N   -  INTEGER     -  LENGTH OF SERIES  -  (SAME)
C	INV -  INTEGER     -  0 FOR DIRECT      -  RETURNS ERROR CODE
C			      1 FOR INVERSE        (-1 IF ERROR)
C
C	DIRECT IS:
C	       N
C	(1/N) SUM X(T)*EXP(-2*PI*I*(T-1)*(J-1)/N), J=1,...,N
C             T=1
C
C	INVERSE IS:
C	       N
C	      SUM X(T)*EXP( 2*PI*I*(T-1)*(J-1)/N), J=1,...,N
C             T=1
C
	DIMENSION XR(N),XI(N),UR(15),UI(15)
	LOGICAL FIRST
	DATA FIRST /.TRUE./
	IF(.NOT.FIRST) GOTO 120
	UR(1)=0.0
	UI(1)=1.0
	DO 110 I=2,15
	UR(I)=SQRT((1.0+UR(I-1))/2.0)
  110	UI(I)=UI(I-1)/(2.0*UR(I))
	FIRST=.FALSE.
  120	IF((N.GT.0).AND.(FLOAT(N).LE.2.0**16)) GOTO 130
	INV=-1
	RETURN
  130	N0=1
	II=0
  140	N0=N0+N0
	II=II+1
	IF(N0.LT.N) GOTO 140
	I1=N0/2
	I3=1
	I0=II
	WRITE(1,7001)
 7001	FORMAT(' FFTing data pass - '$
	DO 260 I4=1,II
	WRITE(1,7002) I4
 7002	FORMAT(' ',I2$
	DO 250 K=1,I1
	WR=1.0
	WI=0.0
	KK=K-1
	DO 230 I=1,I0
	IF(KK.EQ.0) GOTO 240
	IF(MOD(KK,2).EQ.0) GOTO 230
	J0=I0-I
	WS=WR*UR(J0)-WI*UI(J0)
	WI=WR*UI(J0)+WI*UR(J0)
	WR=WS
  230	KK=KK/2
  240	IF(INV.EQ.0) WI=-WI
	L=K
	DO 250 J=1,I3
	L1=L+I1
	ZR=XR(L)+XR(L1)
	ZI=XI(L)+XI(L1)
	Z=WR*(XR(L)-XR(L1))-WI*(XI(L)-XI(L1))
	XI(L1)=WR*(XI(L)-XI(L1))+WI*(XR(L)-XR(L1))
	XR(L1)=Z
	XR(L)=ZR
	XI(L)=ZI
  250	L=L1+I1
	I0=I0-1
	I3=I3+I3
  260	I1=I1/2
	WRITE(1,7003)
 7003	FORMAT('  *')
	UM=1.0
	IF(INV.EQ.0) UM=1.0/FLOAT(N0)
	DO 310 J=1,N0
	K=0
	J1=J-1
	DO 320 I=1,II
	K=2*K+MOD(J1,2)
  320	J1=J1/2
	K=K+1
	IF(K.LT.J) GOTO 310
	ZR=XR(J)
	ZI=XI(J)
	XR(J)=XR(K)*UM
	XI(J)=XI(K)*UM
	XR(K)=ZR*UM
	XI(K)=ZI*UM
  310	CONTINUE
	RETURN
	END
C
	SUBROUTINE DETRND(X,N,K)
C
C	COMPUTES AND SUBTRACTS EITHER THE SERIES MEAN OR A LINEAR
C	LEAST SQUARES TREND FROM THE SERIES DATA, E.G. POLYMONIAL
C	OF DEGREE K
C
C	X  -  REAL ARRAY  -  THE INPUT SERIES -  RETURNS DETRENDED
C	N  -  INTEGER     -  SERIES LENGTH    -  (UNCHANGED)
C	K  -  INTEGER     -  DEGREE OF POLYN  -  (UNCHANGED)
C
	DIMENSION X(N)
	WRITE(1,7001)
 7001	FORMAT(' Detrending data...')
	SUMX=0.0
	DO 20 I=1,N
   20	SUMX=SUMX+X(I)
	XBAR=SUMX/FLOAT(N)
	WRITE(1,7002) XBAR
 7002	FORMAT(5X,'Mean Removed = ',F13.5)
	DO 30 I=1,N
   30	X(I)=X(I)-XBAR
	IF(K.LE.0) RETURN
	TBAR=FLOAT(N+1)/2.0
	SUMTT=(FLOAT(N)*(FLOAT(N)**2-1))/12.0
	SUMTX=0.0
	DO 40 I=1,N
   40	SUMTX=SUMTX+X(I)*(FLOAT(I)-TBAR)
	BETA=SUMTX/SUMTT
	WRITE(1,7003) BETA
 7003	FORMAT(5X,'Slope Removed = ',F13.5)
	DO 50 I=1,N
   50	X(I)=X(I)-BETA*(FLOAT(I)-TBAR)
	RETURN
	END
C
	SUBROUTINE TAPER(X,N,P)
C
C	APPLIES SPLIT-COSINE-BELL TAPERING TO THE TIME SERIES.  THE
C	ARGUMENT P IS THE TOTAL PROPORTION OF THE SERIES TO TAPER.
C	TUKEY SUGGESTS .1 TO .2, .25 IS THE DEFAULT IN MANY PACKAGES.
C
C	X  -  REAL ARRAY  -  THE INPUT SERIES -  RETURNS TAPERED
C	N  -  INTEGER     -  SERIES LENGTH    -  (UNCHANGED)
C	P  -  REAL        -  PROP. TO  TAPER  -  (UNCHANGED)
C
	DATA PI /3.141593/
	DIMENSION X(N)
	IF((P.LE.0.0).OR.(P.GE.1.0)) RETURN
	WRITE(1,7001)
 7001	FORMAT(' Tapering data...')
	M=INT(P*FLOAT(N)+0.5)/2
	DO 10 I=1,M
	WEIGHT=0.5-0.5*COS(PI*(FLOAT(I)-0.5)/FLOAT(M))
	X(I)=X(I)*WEIGHT
   10	X(N+1-I)=X(N+1-I)*WEIGHT
	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(1024)
	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 DATOUT(X,N,START,STEP,M,HEAD)
C
C	OUTPUT A SERIES OF VALUES IN FORMAT REQUIRED FOR DATIN.
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
	WRITE(1,7001)
 7001	FORMAT(' Writing data...')
	WRITE(M,1) HEAD,N,START,STEP
    1	FORMAT(1X,A0/1X,I5/' (5E15.8)'/1X,2F10.5)
	WRITE(M,2) (X(I),I=1,N)
    2	FORMAT(1X,5E15.8)
	RETURN
	END
D,FMT
	WRITE(1,7001)
 7001	FORMAT(' Reading data...')
	READ(M,2) HEAD,N,FMT,START,STEP
    
