      SUBROUTINE FOURIE(X,N)
C
C     PURPOSE--THIS SUBROUTINE PERFORMS A FOURIER ANALYSIS
C              OF THE DATA IN THE INPUT VECTOR X. 
C              THE ANALYSIS CONSISTS OF THE FOLLOWING--
C              1) COMPUTING (AND PRINTING)
C                 (FOR EACH OF THE HARMONIC FREQUENCIES
C                 1/N, 2/N, 3/N, ..., 1/2)
C                 THE CORRESPONDING FOURIER COEFICIENTS,
C                 THE AMPLITUDE, THE PHASE,
C                 THE CONTRIBUTION TO THE TOTAL VARIANCE,
C                 AND THE RELATIVE CONTRIBUTION TO THE TOTAL
C                 VARIANCE.
C              2) PLOTTING OUT A FOURIER LINE SPECTRUM =
C                 THE PERIODOGRAM = THE PLOT OF RELATIVE
C                 CONTRIBUTION TO TOTAL VARIANCE
C                 (AT EACH FOURIER FREQUENCY) VERSUS
C                 THE FOURIER FREQUENCY.
C
C              IN ORDER THAT THE RESULTS OF THE FOURIER ANALYSIS
C              BE VALID AND PROPERLY INTERPRETED, THE INPUT DATA
C              IN X SHOULD BE EQUI-SPACED IN TIME 
C              (OR WHATEVER VARIABLE CORRESPONDS TO TIME).
C
C              THE HORIZONTAL AXIS OF THE SPECTRA PRODUCED
C              BY THIS SUBROUTINE IS FREQUENCY.
C              THIS FREQUENCY IS MEASURED IN UNITS OF
C              CYCLES PER 'DATA POINT' OR, MORE PRECISELY, IN
C              CYCLES PER UNIT TIME WHERE
C              'UNIT TIME' IS DEFINED AS THE
C              ELAPSED TIME BETWEEN ADJACENT OBSERVATIONS.
C              THE RANGE OF THE FREQUENCY AXIS IS 0.0 TO 0.5.
C
C     INPUT ARGUMENTS--X      = THE SINGLE PRECISION VECTOR OF
C                               (UNSORTED) OBSERVATIONS.
C                      N      = THE INTEGER NUMBER OF OBSERVATIONS
C                               IN THE VECTOR X.
C     OUTPUT--2 TO 10 PAGES (DEPENDING ON
C             THE INPUT SAMPLE SIZE) OF 
C             AUTOMATIC PRINTOUT--
C             1) A LISTING OF THE AMPLITUDE,
C                PHASE, CONTRIBUTION TO THE
C                TOTAL VARIANCE, AND RELATIVE
C                CONTRIBUTION TO THE TOTAL
C                VARIANCE FOR EACH OF THE
C                FOURIER FREQUENCIES
C                (1/N, 2/N, 3/N, ..., 1/2).
C                THIS LISTING MAY TAKE AS LITTLE AS 1
C                PAGE OR AS MANY AS N/100 PAGES
C                (THE EXACT NUMBER DEPENDING ON
C                THE INPUT SAMPLE SIZE N).
C                THIS LISTING IS TERMINATED
C                AFTER AT MOST 8 COMPUTER PAGES.
C                IF MORE PAGES ARE DESIRED,
C                CHANGE THE VALUE OF THE
C                VARIABLE     MAXPAG
C                WITHIN THIS SUBROUTINE 
C                FROM 8 TO WHATEVER DESIRED.
C             2) A PLOT OF THE RELATIVE 
C                CONTRIBUTION TO THE
C                TOTAL VARIANCE VERSUS FREQUENCY. 
C     PRINTING--YES.
C     RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N
C                   FOR THIS SUBROUTINE IS 15000. 
C                 --THE SAMPLE SIZE N MUST BE GREATER
C                   THAN OR EQUAL TO 3. 
C     OTHER DATAPAC   SUBROUTINES NEEDED--PLOTSP AND CHSPPF.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, SIN, COS, ATAN.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     COMMENT--FOURIER ANALYSIS DIFFERS FROM SPECTRAL ANALYSIS
C              (AS, FOR EXAMPLE, PRODUCED BY THE DATAPAC
C              TIMESE SUBROUTINE) IN THAT A
C              FOURIER ANALYSIS DOES NO SMOOTHING ON
C              THE SPECTRAL ESTIMATES WHEREAS A SPECRRAL
C              ANALYSIS DOES SMOOTH THE SPECTRAL ESTIMATES. 
C              THE NET RESULT IS THAT THE SPECTRAL
C              ESTIMATES OBTAINED FROM A FOURIER
C              ANALYSIS ARE ALMOST ALWAYS MORE
C              VARIABLE THAN THOSE OBTAINED IN A
C              SPECTRAL ANALYSIS.
C              THE PRACTICAL CONCLUSION IS THAT
C              WHEN THE DATA ANALYST HAS A CHOICE 
C              OF WHETHER TO PERFORM A FOURIER
C              ANALYSIS OR A SPECTRAL ANALYSIS,
C              THE SPECTRAL ANALYSIS SHOULD
C              ALMOST ALWAYS BE PREFERRED.
C            --THE MAXIMUM NUMBER OF FOURIER FREQUENCIES
C              FOR WHICH THE FOURIER COEFFICIENTS IS
C              COMPUTED (AND LISTED) IS N/2 WHERE N IS
C              THE SAMPLE SIZE (LENGTH OF THE
C              DATA RECORD IN THE VECTOR X).
C              THIS RULE IS OVERRIDDEN
C              (FOR LISTING PURPOSES ONLY)
C              IN LARGE DATA SETS AND IS REPLACED 
C              BY THE RULE THAT THE MAXIMUM
C              NUMBER OF LAGS LISTED = 800
C              (WHICH CORRESPONDS TO AN 
C              8-PAGE LISTING OF FOURIER COEFFICIENTS.
C              IF MORE PAGES ARE DESIRED,
C              CHANGE THE VALUE OF THE
C              VARIABLE     MAXPAG
C              WITHIN THIS SUBROUTINE
C              FROM 8 TO WHATEVER DESIRED.
C            --IF THE INPUT OBSERVATIONS IN X ARE CONSIDERED
C              TO HAVE BEEN COLLECTED 1 SECOND APART IN TIME,
C              THEN THE FREQUENCY AXIS OF THE RESULTING
C              SPECTRA WOULD BE IN UNITS OF HERTZ 
C              (= CYCLES PER SECOND).
C            --THE FREQUENCY OF 0.0 CORRESPONDS TO A CYCLE
C              IN THE DATA OF INFINITE (= 1/(0.0))
C              LENGTH OR PERIOD.
C              THE FREQUENCY OF 0.5 CORRESPONDS TO A CYCLE
C              IN THE DATA OF LENGTH = 1/(0.5) = 2 DATA POINTS.
C            --ANY EQUI-SPACED FOURIER ANALYSIS IS
C              INTRINSICALLY LIMITED TO DETECTING FREQUENCIES
C              NO LARGER THAN 0.5 CYCLES PER DATA POINT;
C              THIS CORRESPONDS TO THE FACT THAT THE
C              SMALLEST DETECTABLE CYCLE IN THE DATA
C              IS 2 DATA POINTS PER CYCLE.
C     REFERENCES--JENKINS AND WATTS, ESPECIALLY PAGE 290.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING LABORATORY (205.03)
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE--301-921-2315
C     ORIGINAL VERSION--NOVEMBER  1972. 
C     UPDATED         --NOVEMBER  1975. 
C     UPDATED         --FEBRUARY  1976. 
C
C---------------------------------------------------------------------
C
      CHARACTER*4 ALPERC
      DIMENSION X(1)
      DIMENSION A(7500),B(7500)
      COMMON /BLOCK2/ WS(15000)
      EQUIVALENCE (A(1),WS(1)),(B(1),WS(7501))
      DATA PI/3.14159265358979/
      DATA ALPERC/'%'/
C
      IPR=6
      ILOWER=3
      IUPPER=15000
      MAXPAG=8
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.ILOWER.OR.N.GT.IUPPER)GOTO50
      HOLD=X(1)
      DO65I=2,N
      IF(X(I).NE.HOLD)GOTO90
   65 CONTINUE
      WRITE(IPR, 9)HOLD
      RETURN
   50 WRITE(IPR,17)ILOWER,IUPPER
      WRITE(IPR,47)N
      RETURN
   90 CONTINUE
    9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT (A VECTOR) TO THE FOURIE SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6
     1H *****)
   17 FORMAT(1H , 96H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 FOURIE SUBROUTINE IS OUTSIDE THE ALLOWABLE (,I6,1H,,I6,16H) INTER
     1VAL *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      AN=N
C
C     DETERMINE IF N IS ODD OR EVEN
C
      IEVODD=N-2*(N/2)
      DEL=(AN+1.0)/2.0
      IF(IEVODD.EQ.0)DEL=(AN+2.0)/2.0
C
C     COMPUTE THE SAMPLE MEAN 
C
      SUM=0.0
      DO100I=1,N
      SUM=SUM+X(I)
  100 CONTINUE
      XBAR=SUM/AN
C
C     COMPUTE THE BIASED SAMPLE VARIANCE
C
      SUM=0.0
      DO200I=1,N
      SUM=SUM+(X(I)-XBAR)**2
  200 CONTINUE
      VBIAS=SUM/AN
C
C     COMPUTE THE FOURIER COSINE AND SINE COEFFICIENTS--THEY ARE PLACED
C     IN VECTORS A AND B, RESPECTIVELY. 
C
      NHALF=N/2
      DO400I=1,NHALF
      AI=I
      SUMA=0.0
      SUMB=0.0
      DO500J=1,N
      T=J 
      SUMA=SUMA+X(J)*COS(2.0*PI*(AI/AN)*(T-DEL))
      SUMB=SUMB+X(J)*SIN(2.0*PI*(AI/AN)*(T-DEL))
  500 CONTINUE
      A(I)=SUMA/AN
      B(I)=SUMB/AN
  400 CONTINUE
C
C     WRITE OUT THE SAMPLE SIZE, THE SAMPLE MEAN, 
C     AND THE (BIASED) SAMPLE VARIANCE. 
C
      WRITE(IPR,998)
      WRITE(IPR,801)
      WRITE(IPR,999)
      WRITE(IPR,999)
      WRITE(IPR,805)N
      WRITE(IPR,810)XBAR
      WRITE(IPR,815)VBIAS
      WRITE(IPR,999)
C
C     COMPUTE THE HARMONIC CONTRIBUTION 
C     AT EACH OF THE FOURIER FREQUENCIES.
C     THE FUNDAMENTAL FOURIER FREQUENCY 
C     IS 1/N CYCLES PER DATA POINT
C     (WHERE N = THE INPUT SAMPLE SIZE).
C     THE OTHER FOURIER FREQUENCIES
C     ARE MULTIPLES OR HARMONICS
C     (2/N, 3/N, 4/N, ...1/2) OF THE FUNDAMENTAL. 
C     COMPUTE AMPLITUDES, PHASES, AND
C     CONTRIBUTIONS TO THE VARIANCE AT EACH
C     OF THE FOURIER FREQUENCIES.
C     COMPUTE THE PERCENTAGE CONTRIBUTION
C     TO THE TOTAL VARIANCE AT EACH
C     OF THE FOURIER FREQUENCIES.
C     NOTE--TO SAVE STORAGE, ALSO COPY
C     THE PERCENTAGE CONTRIBUTIONS TO THE VARIANCE)
C     (WHICH WILL LATER BE PLOTTED OUT LIKE A SPECTRUM)
C     INTO THE VECTOR A; THIS WILL DESTROY
C     THE PREVIOUS CONTENTS OF THE VECTOR A.
C     WRITE OUT ALL OF THE ABOVE.
C
      NNPAGE=50
      I=0 
      DO600IPAGE=1,MAXPAG
      WRITE(IPR,998)
      WRITE(IPR,820)
      WRITE(IPR,821)
      WRITE(IPR,822)
      DO700J=1,NNPAGE
      I=I+1
      AI=I
      FFREQ   =AI/AN
      PERIOD=1.0/FFREQ
      ANGRAD   =(AI/AN)*2.0*PI
      ANGDEG   =(AI/AN)*360.0 
      AMP   =SQRT(A(I)*A(I)+B(I)*B(I))
      PHASE1   =ATAN(-B(I)/A(I))
      PHASE2   =PHASE1   *360.0/(2.0*PI)
      CONMSQ   =2.0*AMP   *AMP
      IF(I.EQ.NHALF.AND.IEVODD.EQ.0)CONMSQ=CONMSQ/2.0
      PERCON   =100.0*(CONMSQ   /VBIAS) 
      WRITE(IPR,825)I,FFREQ,PERIOD,A(I),B(I),AMP,PHASE1,PHASE2,
     1CONMSQ,PERCON,ALPERC
      A(I)=PERCON
      IF(I.GE.NHALF)GOTO750
      ISKIP=I-10*(I/10)
      IF(ISKIP.EQ.0)WRITE(IPR,999)
  700 CONTINUE
  600 CONTINUE
  750 CONTINUE
C
C     PLOT OUT THE PERCENTAGE CONTRIBUTIONS
C     TO THE TOTAL VARIANCE AT
C     EACH OF THE FOURIER FREQUENCIES
C     (1/N, 2/N, 3/N, ..., 1/2).
C     THIS WILL CORRESPOND TO A SPECTRAL
C     PLOT IN SPECTRAL ANALYSIS.
C
      CALL PLOTSP(A,NHALF,0)
      WRITE(IPR,855)
C
  801 FORMAT(1H ,44X,16HFOURIER ANALYSIS)
  805 FORMAT(1H ,40X,41HTHE SAMPLE SIZE N                      = ,I8) 
  810 FORMAT(1H ,40X,41HTHE SAMPLE MEAN                        = ,F20.8)
  815 FORMAT(1H ,40X,41HTHE SAMPLE VARIANCE (WITH DIVISOR N-1) = ,F20.8)
  820 FORMAT(1H ,40H     I   FOURIER   PERIOD      FOURIER  ,
     1 30H      FOURIER       AMPLITUDE ,
     1  46H      PHASE          PHASE         VARIANCE   ,
     1   10H  RELATIVE)
  821 FORMAT(1H ,44H        FREQUENCY            COEFFICIENT    ,
     1 11HCOEFFICIENT              ,
     1  59H                    RADIANS        DEGREES        COMPONENT,
     1   12H    VARIANCE)
  822 FORMAT(1H ,43H     (CYCLES/POINT)             A(I)       ,
     1 14H    B(I)      ,
     1  50H                                                  ,
     1   22H         COMPONENT (%))
  825 FORMAT(1H ,I6,2X,F8.6,1X,F8.2,6(1X,E14.7),2X,F6.2,A1) 
  855 FORMAT(1H ,40X,56HPERIODOGRAM = FOURIER LINE SPECTRUM OF THE ORIGI
     1NAL DATA)
  998 FORMAT(1H1)
  999 FORMAT(1H )
C
      RETURN
      END 
