/*
 * Purpose : Perform a one dimensional FFT in place.
 *
 * Parameters : input   [char *] array with input values
 *    npoints [int]  length of array.
 *    dir   [int]  if dir =  1 then forward FFT,
 *     if dir = -1 then inverse FFT.
 *
 *
 *
Purpose : Perform a one dimensional FFT in place.

Parameters:
 input [char *] array with input values
 npoints [int]  length of array.
 dir  [int]  if dir =  1 then forward FFT,
     if dir = -1 then inverse FFT.


*/
#include <math.h>
#include "complex.h"


#define PI 3.14159265358979

int fft(double *real,double *imag,int npoints,int dir)

{
   double temp;
   register int i,j,k;
   int index,swapindex;
   register double tr,ti,angle,wr,wi;

  /* Swap the input elements for the decimation in time algorithm. */

   for (index = 1, swapindex = 0; index < npoints; index++)
   {
      k = npoints;
      do  k /= 2;  while ((swapindex + k) >= npoints);
      swapindex = (swapindex % k) + k;
      if (swapindex <= index) continue;
      temp = real[index];
      real[index] = real[swapindex];
      real[swapindex] = temp;
      temp = imag[index];
      imag[index] = imag[swapindex];
      imag[swapindex] = temp;
   }

  /*
   * Do the butterfly computations.
   *
   *   stage index:   k = 1, 2, 4, 8 . . .
   *     For npoints = 8, for example, there will
   *     be three stages of butterfly computations.
   *   b'fly indices:  i and j will be separated by a distance
   *     dependent on the current stage.
   *     k is used as the separation constant.
   */

   for (k = 1 ; k < npoints ; k *= 2)
   {
      for (index = 0 ; index < k ; index++)
      {
  angle = PI * ((double) (index * -dir)) / ((double) k);
  wr = cos(angle);
  wi = sin(angle);
  for (i = index ; i < npoints ; i += 2*k)
  {
     j = i + k;
     tr = (wr * real[j]) - (wi * imag[j]);
     ti = (wr * imag[j]) + (wi * real[j]);
     
     real[j] = real[i] - tr;
     imag[j] = imag[i] - ti;
     real[i] += tr;
     imag[i] += ti;
  }
      }
   }

   wr = sqrt(1./npoints);

/* scale the output by sqrt(1./N) */

   for (i = 0 ; i < npoints ; i++)
   {
      real[i] *= wr;
      imag[i] *= wr;
   }

   return (0);
}
/*+
 reorders elements of an block for a real fft

 programmer:   L v.d. Wal
 rewritten in C by:  Frits v.d. Putten
 date:    sept 1980
     April 1988 (conversion to C)

 ****************************************************************

 input variables :
  c   complex data
  n   number of elements in c (power of 2)

 output variables :
  c   complex data
-*/

shufft(COMPLEX *c,int n,int dir)

{
   register int i, n2;
   register double scale;
   double d1;
   double d2;
   COMPLEX cda1;
   COMPLEX cda2;
   COMPLEX w1;
   COMPLEX w2;

   d1 = c[0].real;
   d2 = c[0].imag;
   c[0].real = d1+d2;
   c[0].imag = d1-d2;
   w1.real = w2.real = cos(PI/n);
   w1.imag = w2.imag = -dir*sin(PI/n);
   n2 = n/2;
   if (dir == 1)
      scale = 0.5;
   else
      scale = 1.0;

   for (i=1; i<=n2; i++)
   {
     cda1.real  = scale * (c[i].real + c[n-i].real);
     cda1.imag  = scale * (c[i].imag - c[n-i].imag);
     cda2.real  = scale * (c[n-i].imag + c[i].imag);
     cda2.imag  = scale * (c[n-i].real - c[i].real);
     c[n-i].real = cda1.real - dir*(cda2.real*w1.real - cda2.imag*w1.imag);
     c[n-i].imag = dir*(cda2.real*w1.imag + cda2.imag*w1.real) - cda1.imag;
     c[i].real  = cda1.real + dir*(cda2.real*w1.real - cda2.imag*w1.imag);
     c[i].imag  = cda1.imag + dir*(cda2.real*w1.imag + cda2.imag*w1.real);
     d1   = w1.real * w2.real - w1.imag * w2.imag;
     d2   = w1.real * w2.imag + w1.imag * w2.real;
     w1.real  = d1;
     w1.imag  = d2;
   }

   return;
}
#undef PI
