#include "fourier.h"
#include "winexception.h"
#include "utils.h"
#include "ressources.h"
#include "resstring.h"
#include <math.h>

#define SWAP(a, b) tempr=(a); (a)=(b); (b)=tempr

Fourier::Fourier(double *data, unsigned long *nn, int ndim, int isign):
    _data(data-1), _ndim(ndim), _isign(isign), 
    _nn(nn-1)
{
    Check();
}

Fourier::Fourier(const Fourier &rhs): _ndim(rhs._ndim),
    _isign(rhs._isign), _nn(rhs._nn), _data(rhs._data)
{}

Fourier::~Fourier()
{}

Fourier &Fourier::operator=(const Fourier &rhs)
{
    if(&rhs != this)
      {
      _ndim  = rhs._ndim;
      _isign = rhs._isign;
      _data  = rhs._data;
      _nn    = rhs._nn;
      }
    return (*this);
}

// Test if all dimensions are a power of two
void Fourier::Check()
{
    for(int i = 1; i <= _ndim; i++)
       {
       if(!IsPowerTwo(_nn[i]))
         {
         _data = NULL;
         _nn = NULL;
         _ndim = 0;
         _isign = 0;
         throw WinException(ResString(ERROR_FOURIER_SIZE));
         }
       }
}

/* Replaces data by its ndim-dimensional discrete Fourier transform, if isign is input as 1.
nn[1..ndim] is an integer array containing the lengths of each dimension (number of complex
values), which MUST all be powers of 2. data is a real array of length twice the product of
these lengths, in which the data are stored as in a multidimensional complex array: real and
imaginary parts of each element are in consecutive locations, and the rightmost index of the
array increases most rapidly as one proceeds along data. For a two-dimensional array, this is
equivalent to storing the array by rows. If isign is input as -1, data is replaced by its inverse
transform times the product of the lengths of all dimensions.

Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable
files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website
http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America) */

void Fourier::Fourn()
{
int idim;
unsigned long i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
unsigned long ibit, k1, k2, n, nprev, nrem, ntot;
double tempi, tempr;
double theta, wi, wpi, wpr, wr, wtemp;

for(ntot=1,idim=1;idim<=_ndim;idim++) //Compute total number of complex values
    ntot*= _nn[idim];

nprev = 1;

for(idim=_ndim; idim>=1; idim--)
   { //Main loop over the dimensions.
   n=_nn[idim];
   nrem=ntot / (n*nprev);

   ip1=nprev << 1;
   ip2=ip1*n;
   ip3=ip2*nrem;
   i2rev=1;

   for(i2=1;i2<=ip2;i2+=ip1)
      { //This is the bit-reversal section of the routine.
      if(i2 < i2rev)
        {
        for(i1=i2; i1<=i2+ip1-2; i1+=2)
           {
           for(i3=i1; i3<=ip3; i3+=ip2)
              {
              i3rev = i2rev+i3-i2;
              SWAP(_data[i3], _data[i3rev]);
              SWAP(_data[i3+1], _data[i3rev+1]);
              }
           }
       }

   ibit=ip2 >> 1;
   while(ibit >= ip1 && i2rev > ibit)
      {
      i2rev -= ibit;
      ibit >>= 1;
      }

   i2rev += ibit;
   }

ifp1=ip1; // Here begins the Danielson-Lanczos section of the routine.

while(ifp1 < ip2)
   {
   ifp2=ifp1 << 1;
   theta = _isign*6.28318530717959/(ifp2/ip1); // Initialize for the trig. recurrence.
   wtemp = sin(0.5*theta);

   wpr = -2.0*wtemp*wtemp;
   wpi = sin(theta);

   wr = 1.0;
   wi = 0.0;

   for(i3=1; i3<=ifp1; i3+=ip1)
      {
      for(i1=i3; i1<=i3+ip1-2; i1+=2)
         {
         for(i2=i1; i2<=ip3; i2+=ifp2)
            {
            k1 = i2; //Danielson-Lanczos formula:
            k2 = k1+ifp1;

            tempr = wr*_data[k2] - wi*_data[k2+1];
            tempi = wr*_data[k2+1] + wi*_data[k2];

            _data[k2] = _data[k1]-tempr;
            _data[k2+1] = _data[k1+1]-tempi;
            _data[k1] += tempr;
            _data[k1+1] += tempi;
            }
         }
      wr = (wtemp = wr)*wpr - wi*wpi + wr; // Trigonometric recurrence.
      wi = wi*wpr + wtemp*wpi + wi;
      }

   ifp1 = ifp2;
   }

nprev *= n;
  }
}

unsigned long* Fourier::Getnn()
{
    return _nn;
}

int Fourier::Getndim() const
{
    return _ndim;
}

int Fourier::Getisign() const
{
    return _isign;
}
