C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *      
C     SUBROUTINE FOR FOURIER, POWER SPECTRA AND AUTOCORRELATION         
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *      
C                                                                       
C                                             CODED BY Y.OHSAKI         
C                                                                       
C     PURPOSE                                                           
C       TO COMPUTE FOURIER SPECTRUM, POWER SPECTRUM AND/OR AUTOCORRE-   
C       LATION OF A SERIES OF EQUI-SPACED DATA                          
C                                                                       
C     USAGE                                                             
C       CALL FPAC(N,X,ND1,DT,IND,F,G,R,ND2,NFOLD,DF)                    
C                                                                       
C     DESCRIPTION OF AUGUMENTS                                          
C       N      - TOTAL NUMBER OF DATA                                   
C       X(ND1) - EQUI-SPACED DATA                                       
C       ND1    - DIMENSION OF X IN CALLING PROGRAM  ND1.LE.8192         
C       DT     - TIME INCREMENT IN DATA IN SEC                          
C       IND    - 100 FOR FOURIER SPECTRUM                               
C                010 FOR POWER SPECTRUM                                 
C                001 FOR AUTOCORRELATION                                
C       F(ND2) - FOURIER SPECTRUM                                       
C       G(ND2) - POWER SPECTRUM                                         
C       R(ND2) - AUTOCORRELATION                                        
C       ND2    - DIMENSION OF F,G,R IN CALLING PROGRAM                  
C       NFOLD  - TOTAL NUMBER OF DATA IN SPECTRUM AND AUTOCORRELATION   
C       DF     - FREQUENCY INCREMENT IN FOURIER AND POWER SPECTRA IN HZ 
C                                                                       
C     REMARKS                                                           
C      (1) PARAMETER IND IS ADDIBLE. IF, FOR INSTANCE, IND=101(100+001),
C          FOURIER SPECTRUM AND AUTOCORRELATION ARE COMPUTED            
C      (2) ND2.GE.NT/2+1, WHERE NT IS POWER OF 2 EQUAL TO N OR MINIMUM  
C          LARGER THAN N                                                
C      (3) AUTOCORRELATION IS NORMALIZED IN TERMS OF THE MEAN OF DATA   
C          SQUARED                                                      
C      (4) AUTOCORRELATION IS SPACED AT THE SAME INTERVAL AS THE GIVEN  
C          SERIES OF DATA                                               
C                                                                       
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     
C       FAST                                                            
C                                                                       
      SUBROUTINE FPAC(N,X,ND1,DT,IND,F,G,R,ND2,NFOLD,DF)                
C                                                                       
      COMPLEX   C(8192)                                                 
      DIMENSION X(ND1),F(ND2),G(ND2),R(ND2)                             
C                                                                       
C     INITIALIZATION                                                    
C                                                                       
      DO 110 M=1,N                                                      
      C(M)=CMPLX(X(M),0.)                                               
  110 CONTINUE                                                          
      NT=2                                                              
  120 IF(NT.GE.N) GO TO 130                                             
      NT=NT*2                                                           
      GO TO 120                                                         
  130 IF(NT.EQ.N) GO TO 150                                             
      DO 140 M=N+1,NT                                                   
      C(M)=(0.,0.)                                                      
  140 CONTINUE                                                          
  150 NFOLD=NT/2+1                                                      
      T=REAL(NT)*DT                                                     
      DF=1./T                                                           
C                                                                       
C     FOURIER TRANSFORM                                                 
C                                                                       
      CALL FAST(NT,C,8192,-1)                                           
C                                                                       
C     FOURIER SPECTRUM                                                  
C                                                                       
      IF(IND.EQ.001) GO TO 180                                          
      DO 160 K=1,NFOLD                                                  
      F(K)=CABS(C(K))*DT                                                
  160 CONTINUE                                                          
      IF(IND.EQ.100) RETURN                                             
C                                                                       
C     POWER SPECTRUM                                                    
C                                                                       
      IF(IND.EQ.101) GO TO 180                                          
      G(1)=F(1)**2/T                                                    
      DO 170 K=2,NFOLD-1                                                
      G(K)=2.*F(K)**2/T                                                 
  170 CONTINUE                                                          
      G(NFOLD)=F(NFOLD)**2/T                                            
      IF(MOD(IND,10).EQ.0) RETURN                                       
C                                                                       
C     AUTOCORRELATION                                                   
C                                                                       
  180 DO 190 K=1,NT                                                     
      C(K)=C(K)*CONJG(C(K))                                             
  190 CONTINUE                                                          
      CALL FAST(NT,C,8192,+1)                                           
      R0=REAL(C(1))                                                     
      DO 200 J=1,NFOLD                                                  
      R(J)=REAL(C(J))/R0                                                
  200 CONTINUE                                                          
      RETURN                                                            
      END                                                               
