//
// written by Yue Zhao 
//
#include <fftw3.h>
#include <stdio.h>
#include <math.h>

#define PI 3.1415926535897932384626433
#define N 24

double maxval(double *x, double *y, int n);

int main(int argc, char** argv)
{
  // declare variables
  double *in, *out, *x, *u, *uprime, *fmx;
  double xlen, dx, t1;
  fftw_plan plan_r2hc, plan_hc2r;
  int i;

  // allocate memory spaces
  x = (double*) fftw_malloc(sizeof(double) * N);
  u = (double*) fftw_malloc(sizeof(double) * N);
  uprime = (double*) fftw_malloc(sizeof(double) * N);
  fmx = (double*) fftw_malloc(sizeof(double) * N);
  in = (double*) fftw_malloc(sizeof(double) * N);
  out = (double*) fftw_malloc(sizeof(double) * N);

  // make fftw plans
  printf("planning FFTW_R2HC with FFTW_PATIENT\n");
  plan_r2hc = fftw_plan_r2r_1d(N, in, out, FFTW_R2HC, FFTW_PATIENT);
  printf("planning FFTW_HC2R with FFTW_PATIENT\n");
  plan_hc2r = fftw_plan_r2r_1d(N, out, in, FFTW_HC2R, FFTW_PATIENT);
  
  // initialize data
  xlen = 2.0*PI;
  dx = xlen /(double)N;
  for (i=0; i<N; i++) {
    x[i] = i * dx;
    u[i] = exp(sin(x[i]));
    uprime[i] = cos(x[i])*u[i];
    in[i] = u[i]; 
  }
  fmx[0] = 0.0;
  for (i=1; i<N/2; i++) {
    fmx[i]   = -2.0*PI*i/xlen;
    fmx[N-i] = 2.0*PI*i/xlen;
  }
  fmx[N/2] = 0.0;

  // do calculations
  printf("r2hc transform\n");
  fftw_execute_r2r(plan_r2hc, in, out);

  printf("d/dx via fourier multiples\n");
  in[0] = fmx[0]*out[0];
  for (i=1; i<N; i++) {
    in[i] = fmx[i]*out[N-i];
  }

  printf("hc2r tranform");
  fftw_execute_r2r(plan_hc2r, in, out);
  t1 = 1.0/N;
  for (i=0; i<N; i++) {
    out[i] = t1*out[i];
  }

  // output results
  printf(" \n");
  printf("N = %d, error = %E\n", N, maxval(out, uprime, N));
  printf(" \n");

  printf("cleaning up\n");
  fftw_free(in);
  fftw_free(out);
  fftw_free(x);
  fftw_free(u);
  fftw_free(uprime);
  fftw_free(fmx);

  fftw_destroy_plan(plan_r2hc);
  fftw_destroy_plan(plan_hc2r);
  
  return 0;
}


double maxval(double *x, double *y, int n) {
  int i;
  double tmp, res = fabs(x[0] - y[0]);

  for (i=1; i<n; i++) {
    tmp = fabs(x[i] - y[i]);
    if (res < tmp) res = tmp;
  }

  return tmp;
}