/** File fourier.c
 *  Compute Fourier transform
 *  See /usr/local/PV/GRACE_PRIVATE/grace-5.99.1dev5/src/fourier.c
 *  See /usr/local/PV/GRACE_PRIVATE/grace-5.99.1dev5/src/computils.c  1522 
 *  See compwin.c 411
 *  Parameters used to compute the transform are statically stored in "fou"
 *  HAVE_FFTW is #defined if FFTW2 is used, undef if not (e.g. with FFTW3)
 */

#include <stdlib.h>
#include <string.h>
#include <cmath.h>

#include "files.h"
#include "graphs.h"
#include "utils.h"
#include "noxprotos.h"
#include "ssdata.h"

#include "fourier.h"

#ifdef HAVE_FFTW
#  include <fftw.h>
#endif
#ifdef HAVE_FFTW3
#  include <fftw3.h>
#endif

extern void comment_and_legend_set (int gno ,int setno ,char *buf);

static int do_fourier_filter (double *v ,int len ,int filter ,double padding);

static FFTParams fou;

void fourier_set_defaults (void)
{
  fou.invflag 	   = 0;
  fou.xscale  	   = FFT_XSCALE_NU;  //FFT_XSCALE_OMEGA; // 
  fou.norm    	   = FFT_NORM_NONE;
  fou.insegmt 	   = FFT_INPUT_SEG_X;  //FFT_INPUT_SEG_INDEX;  //  
  fou.complexin    = FALSE;
  fou.dcdump       = FALSE;
  fou.oversampling = 1.0;
  fou.round2n      = FALSE;
  fou.filter       = FFT_WINDOW_NONE;
  fou.padding      = 1.0;
  fou.outsegmt     = FFT_OUTSEGMT_CENTERED; // FFT_OUTSEGMT_POSITIVE;  // 
  fou.output       = FFT_OUTPUT_MAGNITUDE;   //FFT_OUTPUT_REIM; // 
}

int fourier_params_get (FFTParams *p)
{
  if (p != NULL) {
    memcpy (p ,&fou ,sizeof(FFTParams));
    return RETURN_SUCCESS;
  } else {
    return RETURN_FAILURE;
  }
}

int fourier_params_set (FFTParams *p)
{
  if (p != NULL) {
    memcpy (&fou ,p ,sizeof(FFTParams));
    return RETURN_SUCCESS;
  } else {
    return RETURN_FAILURE;
  }
}

int fourier_nb_out_sets (FFT_ouput output)
{
  /* number of output sets */
  switch (output) {
  case FFT_OUTPUT_RE:
  case FFT_OUTPUT_IM:
  case FFT_OUTPUT_MAGNITUDE:
  case FFT_OUTPUT_PHASE:
    return 1;
    break;
  case FFT_OUTPUT_APHI:
  case FFT_OUTPUT_COEFF:
  case FFT_OUTPUT_REIM:
    return 2;
    break;
  }
  return -1;
}

/**
 *  the output sets setno2 and  setno22 must exists before the call
 *  adapted from do_fourier in grace-5.99
 *  See /usr/local/PV/GRACE_PRIVATE/grace-5.99.1dev5/src/computils.c  1522 
 *
 */
int do_fourier (int gno1 ,int setno1 ,int gno2 ,int setno2 ,int setno22)
{
  int i ,j ,inlen ,buflen ,outlen ,ns2;
  double *in_x ,*in_re ,*in_im ,*buf_re ,*buf_im ,*out_x ,*out_x1 ,*out_y ,*out_y1;
  double xspace ,amp_correction ,avg ,dummy;
  char buf[256] ,buf1[256];
  double xmin ,xmax ,a ,b ,c ,arg ,dX;
  int imin ,imax;
  double txc ,*out_xc ,*out_yc ,*out_y1c ,Xk ,cAXk ,tre;

  /* get input and perform safety check on it */
  inlen = getsetlength (gno1 ,setno1);
  if (inlen < 2) {
    errmsg("Set length < 2");
    return RETURN_FAILURE;
  }

  in_re = getcol (gno1 ,setno1 ,DATA_Y);
  if (!in_re) {
    /* should never happen */
    return RETURN_FAILURE;
  }
  in_im = NULL;
  if (fou.complexin) in_im =  getcol (gno1 ,setno1 ,DATA_Y1);
  
  in_x = getcol (gno1 ,setno1 ,DATA_X);
  if (monospaced (in_x ,inlen ,&xspace) != TRUE) {
    errmsg("do_fourier error: abscissas not monospaced.\n\t (you can interpolate the set before doing FFT if you want)");
    return RETURN_FAILURE;
  } else {
    if (xspace == 0.0) {
      errmsg("do_fourier: The set spacing is 0, can't continue");
      return RETURN_FAILURE;
    }
  }

  /* Definition constants of the Fourier transform used with FFT_INPUT_SEG_X */
  switch (fou.xscale) {
  case FFT_XSCALE_INDEX:
    a = 1.0;
    c = 1.0;
    break;
  case FFT_XSCALE_NU:
    a = 1.0;
    c = 2.0 * M_PI;
    break;
  case FFT_XSCALE_OMEGA:
    a = 1.0 / sqrt (2.0 * M_PI);
    c = 1.0;
    break;
  default:
    a = 1.0;
    c = 1.0;
    break;
  }
  if (fou.invflag) c = -c;
    

  /* copy input data to buffer storage which will be used then to hold out of fftw */    
  buf_re = copy_data_column (in_re ,inlen);
  if (in_im) {
    buf_im = copy_data_column (in_im ,inlen);
  } else {
    buf_im = xcalloc(inlen, SIZEOF_DOUBLE);
  }
  if (!buf_re || !buf_im) {
    xfree (buf_re);
    xfree (buf_im);
    return RETURN_FAILURE;
  }
    
  /* dump the DC component */
  if (fou.dcdump) {
    stasum (buf_re, inlen, &avg, &dummy);
    for (i = 0; i < inlen; i++) {
      buf_re[i] -= avg;
    }
  }
    
  /* apply data filtering */
  do_fourier_filter (buf_re, inlen, fou.filter, fou.padding);
  if (in_im) {
    do_fourier_filter (buf_im, inlen, fou.filter, fou.padding);
  }
    
  /* a safety measure */
  fou.oversampling = MAX2 (1.0 ,fou.oversampling);
    
  /* round to the closest 2^N ,but NOT throw away any data */
  buflen = (int) rint(inlen*fou.oversampling);
  if (fou.round2n) {
    int i2 = (int) rint (ilog2 ((double) buflen));
    buflen = (int) pow (2.0 ,(double)i2);
    if (buflen < inlen) {
      buflen *= 2;
    }
//    } else {
//      /* it is a nonsense to compute the FFT of an odd number of points,
//       *  thus we add a point to make it even and add half-sum (see below)
//       */
//      if (buflen%2) buflen++;
  }
  if (buflen != inlen) {
    buf_re = xrealloc (buf_re ,SIZEOF_DOUBLE*buflen);
    buf_im = xrealloc (buf_im ,SIZEOF_DOUBLE*buflen);
    if (!buf_re || !buf_im) {
      xfree(buf_re);
      xfree(buf_im);
      return RETURN_FAILURE;
    }
    if (buflen == inlen+1) {
      /* Add the average between first and last values  is better
	 then zero to lower Gibbs oscillations */
      buf_re[inlen] = 0.5 * (buf_re[0] + buf_re[inlen-1]);
      buf_im[inlen] = 0.5 * (buf_im[0] + buf_im[inlen-1]);
    } else {
      /* stuff the added data with zeros */
      for (i = inlen; i < buflen; i++) {
	buf_re[i] = 0.0;
	buf_im[i] = 0.0;
      }
    }
  }
    
  /* EFFECTIVE COMPUTATION of the Fourier transform */
  if (fourier (buf_re ,buf_im ,buflen ,fou.invflag) == RETURN_FAILURE) {
    xfree(buf_re);
    xfree(buf_im);
    return RETURN_FAILURE;
  }

  /* amplitude correction due to the zero padding etc */
  if (fou.insegmt == FFT_INPUT_SEG_X) {
    amp_correction = 1.0;
  } else  {
    amp_correction = (double) buflen/inlen;
    if ((fou.invflag  && fou.norm == FFT_NORM_BACKWARD) ||
	(!fou.invflag && fou.norm == FFT_NORM_FORWARD)) {
      amp_correction /= buflen;
    } else if (fou.norm == FFT_NORM_SYMMETRIC) {
      amp_correction /= sqrt((double) buflen);
    }
  }


  /******* output preparation ********/

  /* number of sets produced */
  ns2 = fourier_nb_out_sets (fou.output);

  /* Output segment is positive half length */
  if (fou.outsegmt == FFT_OUTSEGMT_HALF_LENGTH) {
    outlen = buflen/2 + 1;
    if (fou.complexin) {
      errmsg ("do_fourier error: Output segment positive half length required with complex input");
      return RETURN_FAILURE;
    } else if (fou.output == FFT_OUTPUT_MAGNITUDE) {
      /* DC component */
      buf_re[0] = buf_re[0];
      buf_im[0] = 0.0;
      for (i = 1; i < outlen; i++) {
	/* carefully get amplitude of complex form: 
	   use abs(a[i]) + abs(a[-i]) except for zero term */
	buf_re[i] += buf_re[buflen - i];
	buf_im[i] -= buf_im[buflen - i];
      }
    }
  } else {
    outlen = buflen;
  }

  /* Reallocate set(s) */
  setlength (gno2 ,setno2 ,outlen);
  out_x = getx (gno2 ,setno2);
  out_y = gety (gno2 ,setno2);
  if (ns2 == 2) {
    setlength (gno2 ,setno22 ,outlen);
    out_x1 = getx (gno2 ,setno22);
    out_y1 = gety (gno2 ,setno22);
  } else {
    out_x1 = out_y1 = NULL;
  }

  /* Compute abscissas */
  for (i = 0; i < outlen; i++) {
    switch (fou.xscale) {
    case FFT_XSCALE_NU:
      out_x[i] = (double) i/(xspace*buflen);
      break;
    case FFT_XSCALE_OMEGA:
      out_x[i] = 2*M_PI*i/(xspace*buflen);
      break;
    default:
      out_x[i] = (double) i;
      break;
    }
    if (ns2 == 2) out_x1[i] = out_x[i];
  }

  /* Taking account input translation */
  if (fou.insegmt == FFT_INPUT_SEG_X) {
    minmax (in_x ,inlen ,&xmin ,&xmax ,&imin ,&imax);
    b = xspace;
    arg = c * xmin;
    dX = 2.0 * M_PI / (b * c * buflen);
    for (i = 0; i < outlen; i++) {
      out_x[i] = (double)i * dX;
      Xk = out_x[i];
      if (i >  outlen/2) Xk = (double)(i-outlen) * dX;
      cAXk = c * xmin * Xk;
      switch (fou.output) {
      case FFT_OUTPUT_RE:
	buf_re[i] = a * b * (cos(cAXk)*buf_re[i] + sin(cAXk)*buf_im[i]);
	break;
      case FFT_OUTPUT_IM:
	buf_im[i] = a * b * (cos(cAXk)*buf_im[i] - sin(cAXk)*buf_re[i]);
	break;
      case FFT_OUTPUT_MAGNITUDE:
      case FFT_OUTPUT_PHASE:
      case FFT_OUTPUT_APHI:
      case FFT_OUTPUT_COEFF:
      case FFT_OUTPUT_REIM:
	tre =  cos(cAXk)*buf_re[i] + sin(cAXk)*buf_im[i];
	buf_im[i] = a * b * (cos(cAXk)*buf_im[i] - sin(cAXk)*buf_re[i]);
	buf_re[i] = a * b * tre;
	break;
      }
    }
  }
      
  /* Compute final values */
  for (i = 0; i < outlen; i++) {
    switch (fou.output) {
    case FFT_OUTPUT_MAGNITUDE:
      out_y[i]  = amp_correction * hypot (buf_re[i], buf_im[i]);
      break;
    case FFT_OUTPUT_RE:
      out_y[i]  = amp_correction * buf_re[i];
      break;
    case FFT_OUTPUT_IM:
      out_y[i]  = amp_correction * buf_im[i];
      break;
    case FFT_OUTPUT_PHASE:
      out_y[i]  = atan2 (buf_im[i], buf_re[i]);
      break;
    case FFT_OUTPUT_REIM:
      out_y[i]  = amp_correction * buf_re[i];
      out_y1[i] = amp_correction * buf_im[i];
      break;
    case FFT_OUTPUT_APHI:
      out_y[i]  = amp_correction * hypot (buf_re[i], buf_im[i]);
      out_y1[i] = atan2 (buf_im[i], buf_re[i]);
      break;
    case FFT_OUTPUT_COEFF:        // A FAIRE
      printf ("do_fourier: FFT_OUTPUT_COEFF A FAIRE\n");
      break;
    }
  }

  /* Write appropriate comment */
  switch (fou.output) {
  case FFT_OUTPUT_MAGNITUDE:
    sprintf (buf, "| FFT(G%d.s%d) |", gno1 ,setno1);
    break;
  case FFT_OUTPUT_RE:
    sprintf (buf, "Re FFT(G%d.s%d)", gno1 ,setno1);
    break;
  case FFT_OUTPUT_IM:
    sprintf (buf, "Im FFT(G%d.s%d)", gno1 ,setno1);
    break;
  case FFT_OUTPUT_PHASE:
    sprintf (buf, "Phase FFT(G%d.s%d)", gno1 ,setno1);
    break;
  case FFT_OUTPUT_REIM:
    sprintf (buf  ,"Re FFT(G%d.s%d)" ,gno1 ,setno1);
    sprintf (buf1, "Im FFT(G%d.s%d)", gno1 ,setno1);
    comment_and_legend_set (gno2 ,setno22 ,buf1);
    break;
  case FFT_OUTPUT_APHI:
    sprintf (buf  ,"(A, Phi) FFT(G%d.s%d)" ,gno1 ,setno1);
    break;
  case FFT_OUTPUT_COEFF:        // A FAIRE
    printf ("do_fourier: FFT_OUTPUT_COEFF A FAIRE\n");
    break;
  }
  comment_and_legend_set (gno2 ,setno2 ,buf);
      
  xfree (buf_re);
  xfree (buf_im);

  /* Center the output segment */
  if (fou.outsegmt == FFT_OUTSEGMT_CENTERED) {
    txc =  outlen * out_x[1];
    out_xc     = xmalloc (outlen * SIZEOF_DOUBLE);
    out_yc     = xmalloc (outlen * SIZEOF_DOUBLE);
    if (ns2 == 2) {
      out_y1c = xmalloc (outlen * SIZEOF_DOUBLE);
    } else {
      out_y1c = NULL;
    }
    if (out_xc != NULL && out_yc != NULL) {
      j = 0;
      for (i = outlen/2; i < outlen; i++) {
	out_xc[j] 		 = out_x[i] - txc;
	out_yc[j] 		 = out_y[i];
	if (ns2 == 2) out_y1c[j] = out_y1[i];
	j++;
      }
      for (i = 0; i < outlen/2; i++) {
	out_xc[j]                = out_x[i];
	out_yc[j]                = out_y[i];
	if (ns2 == 2) out_y1c[j] = out_y1[i];
	j++;
      }
      for (i = 0; i < outlen; i++) {
	out_x[i] = out_xc[i];
	out_y[i] = out_yc[i];
	if (ns2 == 2) {
	  out_x1[i] = out_x[i];
	  out_y1[i] = out_y1c[i];
	}
      }
      xfree (out_xc);
      xfree (out_yc);
      if (ns2 == 2) xfree (out_y1c);
    } else {
      return RETURN_FAILURE;
    }
  }
  return RETURN_SUCCESS;
}

/**
 *  do_fourier_filter is the function "apply_window" in grace-5.99
 */
static int do_fourier_filter (double *v ,int len ,int filter ,double padding)
{
  int i;

  if (len < 2 || !v)             return RETURN_FAILURE;
  if (filter == FFT_WINDOW_NONE) return RETURN_SUCCESS;
    
  for (i = 0; i < len; i++) {
    double c ,w ,tmp;
    w = 2 * M_PI * i / (len - 1);
    switch (filter) {
    case FFT_WINDOW_TRIANGULAR:
      c = 1.0 - fabs((i - 0.5*(len - 1))/(0.5*(len - 1)));
      break;
    case FFT_WINDOW_PARZEN:
      c = 1.0 - fabs((i - 0.5*(len - 1))/(0.5*(len + 1)));
      break;
    case FFT_WINDOW_WELCH:
      tmp = (i - 0.5*(len - 1))/(0.5*(len + 1));
      c = 1.0 - tmp*tmp;
      break;
    case FFT_WINDOW_HANNING:
      c = 0.5 - 0.5*cos(w);
      break;
    case FFT_WINDOW_HAMMING:
      c = 0.54 - 0.46*cos(w);
      break;
    case FFT_WINDOW_BLACKMAN:
      c = 0.42 - 0.5*cos(w) + 0.08*cos(2*w);
      break;
    case FFT_WINDOW_FLATTOP:
      c = 0.2810639 - 0.5208972*cos(w) + 0.1980399*cos(2*w);
      break;
#ifdef HAVE_GSL
    case FFT_WINDOW_KAISER:
      tmp = (double) (2*i - len + 1)/(len - 1);
      c = gsl_sf_bessel_I0(padding*sqrt(1.0 - tmp*tmp))/
	gsl_sf_bessel_I0(padding);
      break;
#endif
    default:	/* should never happen */
      c = 0.0;
      return RETURN_FAILURE;
      break;
    }
    
    v[i] *= c;
  }
    
  return RETURN_SUCCESS;
}

/**
 *  do_fourier_command : interface called in  pars.yacc
 *  Adapted from do_fourier_command/do_fourier in grace-5.1.22
 */
void do_fourier_command (int gno1 ,int setno1 ,int invflag ,FFT_filters filter)
{
  fou.invflag   = invflag;
  fou.filter    = filter;
  int setno22;
  int setno2  = nextset (gno1);
  do_copyset (gno1 ,setno1 ,gno1 ,setno2);
  if (fourier_nb_out_sets (fou.output) == 2) {
    setno22 =  nextset (gno1);
    do_copyset (gno1 ,setno1 ,gno1 ,setno22);
  } else {
    setno22 = setno2;
  }
  do_fourier (gno1 ,setno1 ,gno1 ,setno2 ,setno22);
}


void do_fourier_command_1 (int gno1 ,int setno1 			 
			   ,int invflag
			   ,FFT_xscale xscale
			   ,int complexin
			   ,FFT_filters filter
			   ,FFT_ouput output
			   ,int insegmt ,int outsegmt)
{
  fou.output    = output;
  fou.invflag   = invflag;
  fou.xscale    = xscale;
  fou.complexin = complexin;
  fou.filter    = filter;
  if (insegmt  >= 0) fou.insegmt  = insegmt;
  if (outsegmt >= 0) fou.outsegmt = outsegmt;
  int setno22;
  int setno2  = nextset (gno1);
  do_copyset (gno1 ,setno1 ,gno1 ,setno2);
  if (fourier_nb_out_sets (fou.output) == 2) {
    setno22 =  nextset (gno1);
    do_copyset (gno1 ,setno1 ,gno1 ,setno22);
  } else {
    setno22 = setno2;
  }
  do_fourier (gno1 ,setno1 ,gno1 ,setno2 ,setno22);
#ifdef WITH_GTK
  gg_fourier_update ();
#endif
}

/**
 *  used by gg_fext_routine to compute frequency and period
 */
void do_fourier_featext (int gno1 ,int setno1 ,int *setno2)
{
  /* Save current params values */
  FFTParams p = fou;
  /* set values as in grace-5.1.22 */
  fou.invflag      = 0;
  fou.xscale       = FFT_XSCALE_NU;
  fou.norm         = FFT_NORM_NONE;
  fou.insegmt      = FFT_INPUT_SEG_INDEX;
  fou.complexin    = 0;
  fou.dcdump       = 0;
  fou.oversampling = 1.0;
  fou.round2n      = 0;
  fou.filter       = FFT_WINDOW_NONE;
  fou.padding      = 0.0;
  fou.outsegmt     = FFT_OUTSEGMT_HALF_LENGTH;
  fou.output       = FFT_OUTPUT_MAGNITUDE;
  /* create working set and compute Fourier transform */
  *setno2  = nextset (gno1);
  do_copyset (gno1 ,setno1 ,gno1 ,*setno2);
  do_fourier (gno1 ,setno1 ,gno1 ,*setno2 ,*setno2);
  /* restore Fourier params defaults */
  fou = p;
  /* working set will be destroyed in gg_fext_routine */
}


#if defined(HAVE_FFTW) || defined(HAVE_FFTW3)

/* FFTW2  interface was originally written by Marcus H. Mendenhall */

static char *wisdom_file = NULL;
static char *initial_wisdom = NULL;
static int  using_wisdom = FALSE;

static void save_wisdom(void){
    char *final_wisdom;
    final_wisdom = fftw_export_wisdom_to_string();
        if (!initial_wisdom ||
        compare_strings (initial_wisdom, final_wisdom) != TRUE) {
        FILE *wf;
	wf=fopen(wisdom_file,"w");
        if (wf) {
            fftw_export_wisdom_to_file(wf);
	    fclose(wf);
        }
    } 
        fftw_free(final_wisdom);
    if (initial_wisdom) {
        fftw_free(initial_wisdom);
    }
}


static void init_wisdom(void)
{
  static int wisdom_inited = FALSE;
  if (!wisdom_inited)  {
    char *ram_cache_wisdom;
    wisdom_inited = TRUE;
    wisdom_file      = getenv("GRACE_FFTW_WISDOM_FILE");
    ram_cache_wisdom = getenv("GRACE_FFTW_RAM_WISDOM");
    if (ram_cache_wisdom) {
      sscanf(ram_cache_wisdom, "%d", &using_wisdom);
    }
    /* turn on wisdom if it is requested even without persistent storage */
    if (wisdom_file && wisdom_file[0] ) {
      /* if a file was specified in GRACE_FFTW_WISDOM_FILE, try to read it */
      FILE *wf;
#ifdef HAVE_FFTW
      fftw_status fstat;
#else                     /* HAVE_FFTW3 */
      int fstat;
#endif            
      wf = grace_openr (wisdom_file, SOURCE_DISK);
      if (wf) {
	fstat = fftw_import_wisdom_from_file (wf);
	grace_close (wf);
	initial_wisdom = fftw_export_wisdom_to_string ();
      } else {
	initial_wisdom = NULL;
      }
      atexit (save_wisdom);
      /* if a file is specified, always use wisdom */
      using_wisdom = TRUE;
    }
  }
}

/**
 * fourier : compute Fourier transform with FFT hypothesis
 *           iflag=1 backward (exp+), iflag=0 forward (exp-)
 */ 
int fourier(double *jr, double *ji, int n, int iflag)
{
    int i;
    int plan_flags;
    fftw_plan plan;
    fftw_complex *cbuf;
    init_wisdom();
    cbuf = xcalloc(n, sizeof(fftw_complex));
    if (!cbuf) return RETURN_FAILURE;

#ifdef HAVE_FFTW
    for (i = 0; i < n; i++) {
        cbuf[i].re = jr[i];
        cbuf[i].im = ji[i];
    }
    plan_flags = using_wisdom ? (FFTW_USE_WISDOM | FFTW_MEASURE):FFTW_ESTIMATE;
    plan_flags |= FFTW_IN_PLACE;
    plan = fftw_create_plan (n, iflag ? FFTW_BACKWARD:FFTW_FORWARD, plan_flags);
    fftw_one (plan, cbuf, NULL);
    for (i = 0; i < n; i++) {
        jr[i] = cbuf[i].re;
        ji[i] = cbuf[i].im;
    }
#endif

#ifdef HAVE_FFTW3
    for (i = 0; i < n; i++) {
        cbuf[i][0] = jr[i];
        cbuf[i][1] = ji[i];
    }
    plan_flags = using_wisdom ? FFTW_MEASURE : FFTW_ESTIMATE;
    plan = fftw_plan_dft_1d (n ,cbuf ,cbuf ,iflag ? FFTW_BACKWARD:FFTW_FORWARD, plan_flags);
    fftw_execute (plan);
    for (i = 0; i < n; i++) {
        jr[i] = cbuf[i][0];
        ji[i] = cbuf[i][1];
    }
#endif

    xfree(cbuf);
    fftw_destroy_plan (plan);
    
    return RETURN_SUCCESS;
}

#else   /* HAVE_FFTW */

/* Legacy code without FFTW library */

/*
* Bit swapping routine in which the bit pattern of the integer i is reordered.
* See Brigham's book for details
*/
static int bit_swap(int i, int nu)
{
    int ib, i1, i2;

    ib = 0;

    for (i1 = 0; i1 != nu; i1++) {
	i2 = i/2;
	ib = ib*2 + i - 2*i2;
	i = i2;
    }
    return (ib);
}

/*
 * log base 2
 */
//  static int ilog2(int n)
//  {
//      int i = 0;
//      int n1 = n;
//  
//      while (n1 >>= 1) {
//  	i++;
//      }
//      if (1 << i != n) {
//  	return -1;
//      } else {
//  	return i;
//      }
//  }


/*  DFT by definition (modified to work as FFTW) */
int dft(double *jr, double *ji, int n, int iflag)
{
    int i, j, sgn;
    double sumr, sumi, tpi, w, *xr, *xi, on = 1.0 / (double)n;
    double *cov, *siv, co, si;
    int iwrap;

    sgn = iflag ? -1 : 1;
    tpi = 2.0 * M_PI;
    xr = xcalloc(n, SIZEOF_DOUBLE);
    xi = xcalloc(n, SIZEOF_DOUBLE);
    cov = xcalloc(n, SIZEOF_DOUBLE);
    siv = xcalloc(n, SIZEOF_DOUBLE);
    if (xr == NULL || xi == NULL || cov == NULL || siv == NULL) {
	XCFREE(xr);
	XCFREE(xi);
	XCFREE(cov);
	XCFREE(siv);
      return RETURN_FAILURE;
    }
    for (i = 0; i < n; i++) {
	w = tpi * i * on;
	cov[i] = cos(w);
	siv[i] = sin(w)*sgn;
	xr[i] = jr[i];
	xi[i] = ji[i];
    }
    for (j = 0; j < n; j++) {
	sumr = 0.0;
	sumi = 0.0;
	for (i = 0, iwrap=0; i < n; i++, iwrap += j) {
	    if(iwrap >= n) iwrap -= n;
	    co = cov[iwrap];
	    si = siv[iwrap];
	    sumr = sumr + xr[i] * co + sgn * xi[i] * si;
	    sumi = sumi + xi[i] * co - sgn * xr[i] * si;
	}
	jr[j] = sumr;
	ji[j] = sumi;
    }
    XCFREE(xr);
    XCFREE(xi);
    XCFREE(cov);
    XCFREE(siv);
  return RETURN_SUCCESS;
}

/*
   real_data ... ptr. to real part of data to be transformed
   imag_data ... ptr. to imag  "   "   "   "  "      "
   inv ..... Switch to flag normal or inverse transform
   n_pts ... Number of real data points
   nu ...... logarithm in base 2 of n_pts e.g. nu = 5 if n_pts = 32.
*/
// THE RESULTS OBTAINED WITH THIS FUNCTION ARE NOT CORRECT
int fft(double *real_data, double *imag_data, int n_pts, int nu, int inv)
{
    int n2, i, ib, mm, k;
    int sgn, tstep;
    double tr, ti, arg; /* intermediate values in calcs. */
    double c, s;        /* cosine & sine components of Fourier trans. */
    static double *sintab = NULL;
    static int last_n = 0;

    n2 = n_pts / 2;
    
    if (n_pts != last_n) { /* allocate new sin table */
        arg = 2*M_PI/n_pts;
        last_n = 0;
        sintab = xrealloc(sintab, n_pts*SIZEOF_DOUBLE);
        if (sintab == NULL) {
            return RETURN_FAILURE;
        }
        for (i = 0; i < n_pts; i++) {
            sintab[i] = sin(arg*i);
        }
        last_n = n_pts;
    }
    /* sign change for inverse transform  */
    sgn = inv ? 1:-1;
    /* do bit reversal of data in advance */
    for (k = 0; k != n_pts; k++) {
	ib = bit_swap(k, nu);
	if (ib > k) {
	    fswap((real_data + k), (real_data + ib));
	    fswap((imag_data + k), (imag_data + ib));
	}
    }
/* Calculate the components of the Fourier series of the function */
    tstep = n2;
    for (mm = 1; mm < n_pts; mm *= 2) {
        int sinidx = 0, cosidx = n_pts/4;
        for (k=0; k<mm; k++) {
	    c = sintab[cosidx];
	    s = sgn*sintab[sinidx];
	    sinidx += tstep;
            cosidx += tstep;
	    if (sinidx >= n_pts) {
              sinidx -= n_pts;
            }
	    if (cosidx >= n_pts) {
                cosidx -= n_pts;
            }
	    for (i = k; i < n_pts; i += mm*2) {
	        double re1, re2, im1, im2;  
	        re1 = real_data[i];
	        im1 = imag_data[i];
                re2 = real_data[i + mm];
                im2 = imag_data[i + mm];

	        tr = re2*c + im2*s;
	        ti = im2*c - re2*s;
	        real_data[i+mm] = re1 - tr;
	        imag_data[i+mm] = im1 - ti;
	        real_data[i] = re1 + tr;
	        imag_data[i] = im1 + ti;
	    }
        }
        tstep /= 2;
    }
    return RETURN_SUCCESS;
}

// COMMENTED DUE TO PB WITH FFT LEGACY CODE:
//  FOURIER WITHOUT FFTW CALLS ALWAYS dft()

int fourier (double *jr, double *ji, int n, int iflag)
{
//  int i2;
//    if ((i2 = ilog2(n)) > 0) {
//      return fft(jr, ji, n, i2, iflag);
//    } else {
  return dft(jr, ji, n, iflag);
//    }
}

#endif  /* HAVE_FFTW */

