CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C									      C
C Subroutine FASTF							      C
C									      C
C Radix 4 complex discrete fast fourier transform                             C
C									      C
C XREAL = Array which on input contains REAL part of data for transformation  C
C and on output gives REAL part of the result, type REAL, dimension           C
C IABS(ISIZE) or greater.                   				      C
C									      C
C XIMAG = Array which on input contains IMAGINARY part of data for            C
C transformation and on output gives imaginary part of the result, type real, C
C dimension iabs(isize) or greater.               			      C
C									      C
C ISIZE = INTEGER variable or constant specifying length and type of          C
C transform. The length is IABS(ISIZE) which must be a power of 2. Minimum    C
C length 4, maximum 2**20. If ISIZE is positive the forward transform is      C
C calculated, if negative the inverse transform is found.		      C
C									      C
C This subroutine was written by D.Monro. Imperial College, London.   1971    C
C									      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      SUBROUTINE FASTF(XREAL,XIMAG,ISIZE)
C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XREAL(32384),XIMAG(32384)
      INTEGER L(20)                                                      
      EQUIVALENCE(L1,L(1)),(L2,L(2)),(L3,L(3)),(L4,L(4)),                
     1  (L5,L(5)),(L6,L(6)),(L7,L(7)),(L8,L(8)),(L9,L(9)),               
     2  (L10,L(10)),(L11,L(11)),(L12,L(12)),(L13,L(13)),                 
     3  (L14,L(14)),(L15,L(15)),(L16,L(16)),(L17,L(17)),                 
     4  (L18,L(18)),(L19,L(19)),(L20,L(20))                              
      PIE=4.*ATAN(1.)                                                    
      N=IABS(ISIZE)                                                      
      IF(N-4)24,1,1                                                      
C SET UP INITIAL VALUES OF TRANSFORM SPLIT                               
    1 IFACA=N/4                                                          
      IF(ISIZE)2,4,4                                                     
C IF THIS IS TO BE AN INVERSE TRANSFORM, CONJUGATE THE DATA              
    2 DO 3 K=1,N                                                         
    3   XIMAG(K)=-XIMAG(K)                                               
    4 ITIME=0                                                            
    5 IFCAB=IFACA*4                                                      
      ITIME=ITIME+2                                                      
C DO THE TRANSFORMS REQUIRED BY THIS STAGE                               
      Z=PIE/FLOAT(IFCAB)                                                 
      BCOS=-2.*(SIN(Z)**2)                                               
      BSIN=SIN(2.*Z)                                                     
      CW1=1.                                                             
      SW1=0.                                                             
      DO 10 LITLA=1,IFACA                                                
        DO 8 I0=LITLA,N,IFCAB                                            
C THIS IS THE MAIN CALCULATION OF RADIX 4 TRANSFORMS                     
          I1=I0+IFACA                                                    
          I2=I1+IFACA                                                    
          I3=I2+IFACA                                                    
          XS0=XREAL(I0)+XREAL(I2)                                        
          XS1=XREAL(I0)-XREAL(I2)                                        
          YS0=XIMAG(I0)+XIMAG(I2)                                        
          YS1=XIMAG(I0)-XIMAG(I2)                                        
          XS2=XREAL(I1)+XREAL(I3)                                        
          XS3=XREAL(I1)-XREAL(I3)                                        
          YS2=XIMAG(I1)+XIMAG(I3)                                        
          YS3=XIMAG(I1)-XIMAG(I3)                                        
          XREAL(I0)=XS0+XS2                                              
          XIMAG(I0)=YS0+YS2                                              
          X1=XS1+YS3                                                     
          Y1=YS1-XS3                                                     
          X2=XS0-XS2                                                     
          Y2=YS0-YS2                                                     
          X3=XS1-YS3                                                     
          Y3=YS1+XS3                                                     
          IF(LITLA-1)24,6,7                                              
    6     XREAL(I2)=X1                                                   
          XIMAG(I2)=Y1                                                   
          XREAL(I1)=X2                                                   
          XIMAG(I1)=Y2                                                   
          XREAL(I3)=X3                                                   
          XIMAG(I3)=Y3                                                   
          GO TO 8                                                        
C MULTIPLY BY TWIDDLE FACTORS IF REQUIRED                                
    7     XREAL(I2)=X1*CW1+Y1*SW1                                        
          XIMAG(I2)=Y1*CW1-X1*SW1                                        
          XREAL(I1)=X2*CW2+Y2*SW2                                        
          XIMAG(I1)=Y2*CW2-X2*SW2                                        
          XREAL(I3)=X3*CW3+Y3*SW3                                        
          XIMAG(I3)=Y3*CW3-X3*SW3                                        
    8     CONTINUE                                                       
        IF(LITLA-IFACA)9,10,24                                           
C CALCULATE A NEW SET OF TWIDDLE FACTORS                                 
    9   Z=CW1*BCOS-SW1*BSIN+CW1                                          
        SW1=BCOS*SW1+BSIN*CW1+SW1                                        
        TEMPR=1.5-0.5*(Z*Z+SW1*SW1)                                      
        CW1=Z*TEMPR                                                      
        SW1=SW1*TEMPR                                                    
        CW2=CW1*CW1-SW1*SW1                                              
        SW2=2.*CW1*SW1                                                   
        CW3=CW1*CW2-SW1*SW2                                              
        SW3=CW1*SW2+CW2*SW1                                              
   10   CONTINUE                                                         
      IF(IFACA-1)14,14,11                                                
C SET UP THE TRANSFORM SPLIT FOR THE NEXT STAGE                          
   11 IFACA=IFACA/4                                                      
      IF(IFACA)24,12,5                                                   
C THIS IS THE CALCULATION OF A RADIX TWO STAGE                           
   12 DO 13 K=1,N,2                                                      
        TEMPR=XREAL(K)+XREAL(K+1)                                        
        XREAL(K+1)=XREAL(K)-XREAL(K+1)                                   
        XREAL(K)=TEMPR                                                   
        TEMPR=XIMAG(K)+XIMAG(K+1)                                        
        XIMAG(K+1)=XIMAG(K)-XIMAG(K+1)                                   
   13   XIMAG(K)=TEMPR                                                   
      ITIME=ITIME+1                                                      
   14 IF(ISIZE)15,19,17                                                  
C IF THIS WAS AN INVERSE TRANSFORM, CONJUGATE THE RESULT                 
   15 DO 16 K=1,N                                                        
   16   XIMAG(K)=-XIMAG(K)                                               
      GO TO 19                                                           
C IF THIS WAS A FORWARD TRANSFORM, SCALE THE RESULT                      
   17 Z=1./FLOAT(N)                                                      
      DO 18 K=1,N                                                        
        XREAL(K)=XREAL(K)*Z                                              
   18   XIMAG(K)=XIMAG(K)*Z                                              
C UNSCRAMBLE THE RESULT                                                  
   19 I1=20-ITIME                                                        
      DO 20 K=1,I1                                                       
   20   L(K)=1                                                           
      II=1                                                               
      I1=I1+1                                                            
      DO 21 K=I1,20                                                      
        II=II*2                                                          
   21   L(K)=II                                                          
      II=1                                                               
      DO 23 J1=1,L1                                                      
       DO 23 J2=J1,L2,L1                                                 
        DO 23 J3=J2,L3,L2                                                
         DO 23 J4=J3,L4,L3                                               
          DO 23 J5=J4,L5,L4                                              
           DO 23 J6=J5,L6,L5                                             
            DO 23 J7=J6,L7,L6                                            
             DO 23 J8=J7,L8,L7                                           
              DO 23 J9=J8,L9,L8                                          
               DO 23 J10=J9,L10,L9                                       
                DO 23 J11=J10,L11,L10                                    
                 DO 23 J12=J11,L12,L11                                   
                  DO 23 J13=J12,L13,L12                                  
                   DO 23 J14=J13,L14,L13                                 
                    DO 23 J15=J14,L15,L14                                
                     DO 23 J16=J15,L16,L15                               
                      DO 23 J17=J16,L17,L16                              
                       DO 23 J18=J17,L18,L17                             
                        DO 23 J19=J18,L19,L18                            
                         DO 23 J20=J19,L20,L19                           
                          IF(II-J20)22,23,23                             
   22                     TEMPR=XREAL(II)                                
                          XREAL(II)=XREAL(J20)                           
                          XREAL(J20)=TEMPR                               
                          TEMPR=XIMAG(II)                                
                          XIMAG(II)=XIMAG(J20)                           
                          XIMAG(J20)=TEMPR                               
   23                     II=II+1                                        
   24 RETURN                                                             
      END