C  * * * * * * * * * * * * * * * * * * * * * *                          
C     SUBROUTINE FOR FAST FOURIER TRANSFORM                             
C  * * * * * * * * * * * * * * * * * * * * * *                          
C                                                                       
C                         CODED BY Y.OHSAKI                             
C                                                                       
C     PURPOSE                                                           
C       TO PERFORM FAST FOURIER OR INVERSE FOURIER TRANSFORM OF         
C       A SERIES OF EQUI-SPACED DATA                                    
C                                                                       
C     USAGE                                                             
C       CALL FAST(N,X,ND,IND)                                           
C                                                                       
C     DESCRIPTION OF ARGUMENSTS                                         
C       N     - TOTAL NUMBER OF COMPLEX DATA AND TRANSFORMED VALUES     
C       X(ND) - EQUI-SPACED COMPLEX DATA/TRANSFORMED VALUES AT          
C               CALL/RETURN                                             
C       ND    - DIMENSION OF X IN CALLING PROGRAM                       
C       IND   - IND=-1 FOR FOURIER TRANSFORM                            
C                   +1 FOR INVERSE FOURIER TRANSFORM                    
C                                                                       
C     REMARKS                                                           
C       (1) N MUST BE EQUAL TO POWER OF 2                               
C       (2) WHEN IND=-1, TRANSFORMED VALUES ARE MULTIPLIED BY N         
C       (3) EXAMPLE OF CALLING                                          
C               COMPLEX   A(1024)                                       
C               DIMENSION DATA(1024)                                    
C               DATA      NN/1024/                                      
C               DO 1 M=1,NN                                             
C               A(M)=CMPLX(DATA(M),0.0)                                 
C             1 CONTINUE                                                
C               CALL FAST(NN,A,1024,-1)                                 
C                                                                       
C     REFERENCE                                                         
C       N.M.BRENNER/ THREE FORTRAN PROGRAMS THAT PERFORM THE COOLEY-    
C       TUKEY FOURIER TRANSFORM / MIT JULY 1967                         
C                                                                       
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     
C       NONE                                                            
C                                                                       
      SUBROUTINE FAST(N,X,ND,IND)                                       
C                                                                       
      COMPLEX X(ND),TEMP,THETA                                          
C                                                                       
      J=1                                                               
      DO 140 I=1,N                                                      
      IF(I.GE.J) GO TO 110                                              
      TEMP=X(J)                                                         
      X(J)=X(I)                                                         
      X(I)=TEMP                                                         
  110 M=N/2                                                             
  120 IF(J.LE.M) GO TO 130                                              
      J=J-M                                                             
      M=M/2                                                             
      IF(M.GE.2) GO TO 120                                              
  130 J=J+M                                                             
  140 CONTINUE                                                          
      KMAX=1                                                            
  150 IF(KMAX.GE.N) RETURN                                              
      ISTEP=KMAX*2                                                      
      DO 170 K=1,KMAX                                                   
      THETA=CMPLX(0.0,3.141593*REAL(IND*(K-1))/REAL(KMAX))              
      DO 160 I=K,N,ISTEP                                                
      J=I+KMAX                                                          
      TEMP=X(J)*CEXP(THETA)                                             
      X(J)=X(I)-TEMP                                                    
      X(I)=X(I)+TEMP                                                    
  160 CONTINUE                                                          
  170 CONTINUE                                                          
      KMAX=ISTEP                                                        
      GO TO 150                                                         
      END                                                               
