Forward FFT of N-dimensional complex array

The supplied data
values in X and Y are replaced by the coefficients of the Fourier transform of the supplied data. The
coefficients are normalised so that a subsequent call to PDA_NFFTB to perform a backward FFT
would restore the original data values.

The multi-dimensional FFT is implemented using one-dimensional FFTPACK routines. First each row (i.e. a line of pixels parallel to the first axis) in the supplied array is transformed, the Fourier coefficients replacing the supplied data. Then each column (i.e. a line of pixels parallel to the second axis) is transformed. Then each line of pixels parallel to the third axis is transformed, etc. Each dimension is transformed in this way. Most of the complications in the code come from needing to work in an unknown number of dimensions. Two addressing systems are used for each pixel; 1) the vector (i.e. one-dimensional) index into the supplied arrays, and 2) the corresponding Cartesian pixel indices.

CALL PDA_NFFTF( NDIM,
DIM, X, Y, WORK, ISTAT )

The number of
dimensions. This should be no more than 20.

The size of
each dimension.

Supplied holding the real parts of
the complex data values. Returned holding the real parts of the Fourier coefficients. The
array should have the number of elements implied by NDIM ande DIM.

Supplied holding the imaginary parts of the complex data values.
Returned holding the imaginary parts of the Fourier coefficients. The array should have
the number of elements implied by NDIM ande DIM.

A work array. This should have at least ( 6*DimMax + 15 ) elements where
DimMax is the maximum of the values supplied in DIM.

If the value of NDIM is greater than 20 or less than 1, then ISTAT is returned equal to 1,
and the values in X and Y are left unchanged. Otherwise, ISTAT is returned equal to 0.

A double precision version PDA_DNFFTF of the routine exists.