C    -------------------------------------------------------------------------
C    |  File name:   fftd.for              Sean O'Connor             9-22-87 |
C    |                                                               7-01-17 |
C    ------------------------------------------------------------------------- 

C    -------------------------------------------------------------------------
C    |  VAX 11/785       VMS 4.5                   VAX EXTENDED FORTRAN 77   |
C    |  MacBook Pro      Ubuntu LTS 16.04 Linux    gfortran                  |
C    -------------------------------------------------------------------------


C  ****************************************************************************
C  *                                                                          *
C  *                            FFTD                                          *
C  *                                                                          *
C  ****************************************************************************
C
C
C  DESCRIPTION
C 
C      This program tests the subroutine FFT in the file FFT.FOR  It does a 
C      forward discrete Fourier transform (DFT).  It then does an inverse
C      transform of the forward transform to check the roundoff error.
C
C
C  INPUT 
C
C      The number of complex points in the transform.
C      Whether to do a forward (1) or inverse (-1) transform.
C      The data points (complex numbers).
C
C
C  OUTPUT
C
C      The DFT of the input.
C      The inverse DFT of the first DFT.
C
C
C  EXAMPLE 
C
C      $ ./FFTD
C
C      Enter number of complex points in transform: 4
C
C      Enter 1 (forward) or -1 (inverse) transform:  1
C
C      Enter point no.   0 (Real, Imag) >  (1,1)
C      Enter point no.   1 (Real, Imag) >  (2,2)
C      Enter point no.   2 (Real, Imag) >  (3,3)
C      Enter point no.   3 (Real, Imag) >  (4,4)
C
C      num_points =          4
C      direction  =          1
C
C      Point No.   0 = (   1.0000000000000       1.0000000000000     ) 
C      Point No.   1 = (   2.0000000000000       2.0000000000000     ) 
C      Point No.   2 = (   3.0000000000000       3.0000000000000     ) 
C      Point No.   3 = (   4.0000000000000       4.0000000000000     ) 
C
C      The transform of the input is . . .
C 
C      Point No.   0= (  5.0000000000000       5.0000000000000     ) 
C      Point No.   1= ( -2.0000000000000      0.00000000000000E+000) 
C      Point No.   2= ( -1.0000000000000      -1.0000000000000     ) 
C      Point No.   3= (-0.59604644775391E-007 -2.0000000000000     ) 
C
C      The inverse transform of the transform above is . . .
C  
C      Point No.   0= (  1.0000000000000       1.0000000000000     ) 
C      Point No.   1= (  2.0000000000000       2.0000000000000     ) 
C      Point No.   2= (  3.0000000000000       3.0000000000000     ) 
C      Point No.   3= (  4.0000000000000       4.0000000000000     ) 
C
C
C  INSTALLATION
C
C      To compile, type
C
C          $ gfortran FFTD.FOR FORTRAN FFT.FOR -o FFT
C
C      To run, type
C
C          $ ./FFTD
C       
C 
C  ****************************************************************************
C  *                                                                          *
C  ****************************************************************************

      PROGRAM fftd

C  -------------------------------------------------------------------------
C  |                            MACRO DEFINITIONS                          |
C  ------------------------------------------------------------------------- 

      PARAMETER ( MAX_NUM_POINTS = 255 )   !  Maximum number of points in FFT.


C -------------------------------------------------------------------------
C |                          LOCAL VARIABLES                              |
C -------------------------------------------------------------------------

        INTEGER i,             !  Loop counter
     $          num_points,    !  Number of points in the transform
     $          direction      !  Direction of transform 1 for forward,
                               ! -1 for inverse

        COMPLEX X( 0 : MAX_NUM_POINTS )  !  Array of complex data points
                                         !  to be transformed.


C -------------------------------------------------------------------------
C |                          PROGRAM CODE                                 |
C -------------------------------------------------------------------------

C
C    Read the number of points in the FFT
C
        WRITE( 6, 100 )
100     FORMAT ( ' Enter number of complex points in transform:  ' )

        READ ( 5 , * ) num_points      ! Free format input



C
C       Decide whether to do the forward or inverse FFT
C

        WRITE( 6, 200 )
200     FORMAT ( //, ' Enter 1 for a forward transform ', /,
     $               ' or -1 for an inverse transform: ' )

        READ ( 5, * ) direction     ! Free format input


C
C       Check to see that the number of points in the FFT
C       is valid.  If not, halt.
C

        IF ( (num_points .LE. 0) .OR.
     $       (num_points .GT. MAX_NUM_POINTS) ) THEN

                WRITE( 6, 300 )
300             FORMAT( ' The number of points is out of range ')
                STOP

        END IF



C
C       Read in the complex data points, one by one 
C
        DO i = 0, num_points - 1

            WRITE( 6, 400 ) i
400         FORMAT ( ' Enter point no. ', I3,' (Real, Imag) > ')
            READ ( 5 , * ) X( i )      ! Free format input

        END DO


C
C       Print out the points just read in.
C

        WRITE( 6, 500 ) num_points
500     FORMAT ( //, ' num_points = ', I10, / )

        WRITE( 6, 600 ) direction
600     FORMAT ( ' direction =  ', I10, // )

        DO i = 0, num_points - 1

            WRITE( 6 , 700 ) i , X ( i )
700         FORMAT ( ' Point no. ', I3, ' = ( ', 2G22.14E3, ') ' )

        END DO


C
C       Call the FFT subroutine.
C
        CALL FFT( num_points, X, direction )


C
C       Print the transform results.
C
        WRITE( 6, 800 )
800     FORMAT ( //, ' The transform of the input data is . . . ', // )


        DO i = 0, num_points - 1

            WRITE( 6, 900 ) i , X( i )
900         FORMAT (' Point No. ', I3, '= (', 2G22.14E3, ') ' )

        END DO


C
C       Call the FFT subroutine with -direction to do the inverse transform.
C
        CALL FFT( num_points, X, -direction )



C
C       Print the inverse transform coefficients.
C
        WRITE ( 6, 1000 )
1000    FORMAT ( //,' The inverse transform ',
     $               'of the transform above is . . .', // )


        DO i = 0, num_points - 1

            WRITE ( 6, 1100 ) i , X( i )
1100        FORMAT (' Point No. ', I3, '= (', 2G22.14E3, ') ' )

        END DO



C
C  Test the error checking in the FFT routine.
C

C  Number of points is too large.  Currently, MAX_TABLE_SIZE in
C  subroutine FFT is 20, so we can transform up to 2 ** 21 points.

        PRINT *, ' Calling FFT with too many points.'

        CALL FFT( 2 ** 22, X, 1 )


C  Number of points is less than 2.

        PRINT *, ' Calling FFT with 1 point.'

        CALL FFT( 1, X, 1 )


C  Calling FFT with 2^ 15 + 1 points (not a power of 2)

        PRINT *, ' Calling FFT with 2 ^ 15 + 1 points'

        CALL FFT( 2 * 15 + 1, X, 1 )


C  Direction is not 1 or -1.

        PRINT *, ' Calling FFT with direction = 2'

        CALL FFT( 16, X, 2)


        STOP
        END   ! of program fftd