#ifndef __FFTW_H
#define __FFTW_H

#include <boost/multi_array.hpp>
#include <complex>
#include <functional>

namespace mimas {

/** @defgroup fourierTransforms Fourier transformations using boost::multi_array
    Several functions offered by the fftw3 library are wrapped using
    \c boost::multi_array.
    There are several types of fourier-transforms. The current implementation
    only allows fourier-transforms on double precision numbers.
    \li \c fft is the complex-to-complex fourier transform.
    \li \c invfft is the inverse complex-to-complex fourier transform.
    \li \c rfft is the real-to-complex fourier transform (last dimension
    of input-array must be even).
    \li \c invrfft is the inverse real-to-complex fourier transform.

    The real-to-complex fourier transform will result in an array of about
    half the size, because the fourier transform of a real function always is
    symmetric:
    \f$\mathrm{mm\_rfft}:\mathbf{R}^{n_1\times n_2\times\ldots n_d}\rightarrow
    \mathbf{C}^{n_1\times n_2\times\ldots (n_d/2+1)}\f$

    \include fftw/main.cc

    Have a look at the <A HREF="info:/fftw">fftw info-pages</A>.
    Also see <A HREF="http://www.fftw.org">fftw docs</A> and
    <A HREF="http://www.cs.rug.nl/~cosmin/tip/">TiP-0.0.1</A>.
    @author Bala Amavasai (bpa@amavasai.org)
    @author Jan Wedekind (jan@wedesoft.de)
    @date Wed Aug 04 22:48:27 UTC 2004
    @{ */
/** Discrete fourier transform using fftw.
    @see fft */
template< typename Real, std::size_t Dim, typename TPtr >
struct fft_t:
public std::unary_function< boost::const_multi_array_ref< std::complex< Real >, Dim, TPtr >,
                              boost::multi_array< std::complex< Real >, Dim > >
{
  ///
  boost::multi_array< std::complex< Real >, Dim > operator()
    ( const boost::const_multi_array_ref< std::complex< Real >, Dim, TPtr > &field );
};

/** Discrete fourier transform of boost::multi_array.
    @see fft_t */
template< typename Real, std::size_t Dim, typename TPtr >
boost::multi_array< std::complex< Real >, Dim >
fft( const boost::const_multi_array_ref< std::complex< Real >, Dim, TPtr > &field )
{
  return fft_t< Real, Dim, TPtr >()( field );
}

/** Inverse discrete fourier transform using fftw.
    @see invfft */
template< typename Real, std::size_t Dim, typename TPtr >
struct invfft_t:
public std::unary_function< boost::const_multi_array_ref< std::complex< Real >, Dim, TPtr >,
                              boost::multi_array< std::complex< Real >, Dim > >
{
  boost::multi_array< std::complex< Real >, Dim > operator()
    ( const boost::const_multi_array_ref< std::complex< Real >, Dim, TPtr > &field );
};

/** Inverse discrete fourier transform of boost::multi_array.
    @see invfft_t */
  template< typename Real, std::size_t Dim, typename TPtr >
boost::multi_array< std::complex< Real >, Dim >
  invfft( const boost::const_multi_array_ref< std::complex< Real >, Dim, TPtr > &field )
{
  return invfft_t< Real, Dim, TPtr >()( field );
}

/** Discrete fourier transform of real values using fftw.
    The fourier transform of a real (multidimensional) function is
    point-symmetric. Thus it is sufficient to compute only half of the
    values of the corresponding discrete fourier transform.
    
    If the input array is an element of \f$\mathbf{R}^{n_1\times n_2\times
    \ldots\times n_d}\f$, then the result will be element of
    \f$\mathbf{C}^{n_1\times n_2\times\ldots\times (n_d/2+1)}\f$.

    In this implementation \f$n_d\f$ always must be even, so that no additional
    information for a later inverse transform is required!
    @see rfft */
template< typename Real, std::size_t Dim, typename TPtr >
struct rfft_t:
public std::unary_function< boost::const_multi_array_ref< Real, Dim, TPtr >,
                              boost::multi_array< std::complex< Real >, Dim > >
{
  ///
  boost::multi_array< std::complex< Real >, Dim > operator()
    ( const boost::const_multi_array_ref< Real, Dim, TPtr > &field );
};

/** Discrete fourier transform of real-valued boost::multi_array.
    @see rfft_t */
template< typename Real, std::size_t Dim, typename TPtr >
boost::multi_array< std::complex< Real >, Dim >
  rfft( const boost::const_multi_array_ref< Real, Dim, TPtr > &field )
{
  return rfft_t< Real, Dim, TPtr >()( field );
}


/** Inverse discrete fourier transform of symmetric complex array using fftw.
    The fourier transform of a real (multidimensional) function is
    point-symmetric. Thus it is sufficient to compute only half of the
    values of the corresponding discrete fourier transform.
    
    If the input array is an element of \f$\mathbf{C}^{n_1\times n_2\times
    \ldots\times n_d}\f$, then the result will be element of
    \f$\mathbf{R}^{n_1\times n_2\times\ldots\times (n_d-1)*2)}\f$.

    In this implementation the nth dimension of the resulting array is assumed
    to be even!
    @see invrfft */
template< typename Real, std::size_t Dim, typename TPtr >
struct invrfft_t:
    public std::unary_function< boost::const_multi_array_ref< std::complex< Real >, Dim, TPtr >,
                              boost::multi_array< Real, Dim > >
{
  boost::multi_array< Real, Dim > operator()
    ( const boost::const_multi_array_ref< std::complex< Real >, Dim, TPtr > &field );
};

/** Inverse discrete fourier resulting in real-valued boost::multi_array.
    @see invrfft_t */
template< typename Real, std::size_t Dim, typename TPtr >
boost::multi_array< Real, Dim >
invrfft( const boost::const_multi_array_ref< std::complex< Real >, Dim, TPtr > &field )
{
  return invrfft_t< Real, Dim, TPtr >()( field );
}

///@}

};

#include "fourier.tcc"

#endif