#include <cmath>
#include <rfftw.h>
// #ifdef USE_GTK
// #include <gtkmm/main.h>
// #endif
#include <iostream>
#include <math.h> // for HUGE_VAL,M_PI
#include <stdlib.h>

using namespace std;
//#define GAUSSIAN
#define NPTS 1024
// #define INDEF

// #include "/home/faculty/johnt/math224/classes/graphics/GTKColormap.H"
// #include "/home/faculty/johnt/math224/classes/graphics/XYGraphTool.H"

#include "GTKColormap.H"
#include "XYGraphTool.H"


/*double funcn(double x) { 
	return x >= 0.;
}*/

int main(int argc,char** argv) {
    double xmin=-1.*M_PI;
    double xmax=M_PI;

    double ymin=-1.;
    double ymax=1.;

    int nmin=2;
    int nmax=8;
    int power2=__cmath_power(2,nmax);

    Palette pal;
    XYGraphTool::WINDOW_TYPE::COLOR_MAP_TYPE cmap(&pal);
    XYGraphTool gt("interpolation","x","y",xmin,xmax,ymin,ymax,&cmap,0,0.5);
    gt.setbgColor("white");

    fftw_complex *f=OPERATOR_NEW_BRACKET(fftw_complex,power2);
    fftw_complex *fhat=OPERATOR_NEW_BRACKET(fftw_complex,power2);
    fftw_complex *F=OPERATOR_NEW_BRACKET(fftw_complex,NPTS);
    fftw_complex *FHAT=OPERATOR_NEW_BRACKET(fftw_complex,NPTS);

    double *errors=OPERATOR_NEW_BRACKET(double,nmax+1);

    double errors_min=DBL_MAX;
    double errors_max=-DBL_MAX;
    double log10=log(10.);
    for (int i=0;i<=nmax;i++) {
      errors[i]=HUGE_VAL;
    }

    fftw_plan plan_backward=
      fftw_create_plan(NPTS,FFTW_BACKWARD,FFTW_ESTIMATE);

    for (int n=4,nm=nmin;nm<=nmax;n*=2,nm++) {

#ifdef INDEF
    for (int i=0;i<power2;i++) {
       f[i].re=f[i].im=HUGE_VAL;
    }

	for (int i=0;i<NPTS;i++) {
		FHAT[i].re=FHAT[i].im=HUGE_VAL;
	}
#endif // INDEF

	
	double dx=(xmax-xmin)/static_cast<double>(n);
      double x=xmin;
      for (int j=0;j<n;j++,x+=dx) {
// 		f[j].re=funcn(x);
		f[j].re=(x>=0);
        	f[j].im=0.;
	}

	fftw_plan plan_forward=
        fftw_create_plan(n,FFTW_FORWARD,FFTW_ESTIMATE);
      fftw_one(plan_forward,f,fhat);

	// Now fhat contains the frequency amplitudes computed from f

	fftw_real factor=1./static_cast<fftw_real>(n);

// 	for (int j=0;j<=n/2;j++) {
//         FHAT[j].re=fhat[j].re*factor;
//         FHAT[j].im=fhat[j].im*factor;
//       }
//       for (int j=n/2+1;j<=NPTS-n/2;j++) FHAT[j].re=FHAT[j].im=0.;
//       for (int j=1;j<n/2;j++) {
//         FHAT[NPTS-j].re=fhat[n-j].re*factor;
//         FHAT[NPTS-j].im=fhat[n-j].im*factor;
//       }

	for (int j=0;j<=n/4;j++) {
		FHAT[j].re=fhat[j].re*factor;
        	FHAT[j].im=fhat[j].im*factor;
      }
      for (int j=n/4+1;j<=NPTS-n/4;j++) FHAT[j].re=FHAT[j].im=0.;
      for (int j=1;j<n/4;j++) {
		FHAT[NPTS-j].re=fhat[n-j].re*factor;
        	FHAT[NPTS-j].im=fhat[n-j].im*factor;
      }

	fftw_one(plan_backward,FHAT,F);

	// now F contains the function represented by frequencies from FHAT.

	dx=(xmax-xmin)/static_cast<double>(NPTS);
//    if (n==4) {
        x=xmin;
        gt.newPage();
        gt.setfgColor("black");
        gt.drawAxes();
        gt.setfgColor("red");
        for (int j=0;j<NPTS;j++,x+=dx) {
//           double y=funcn(x);
          double y=(x>=0);
          if (j==0) gt.movePen(x,y);
          else gt.drawLine(x,y);
        }
//    }
      double maxerr=0.;
//    gt.setfgColor(n-nmin+1,nmax-nmin);
      gt.setfgColor("blue");
      x=xmin;
      for (int j=0;j<NPTS;j++,x+=dx) {
	fftw_real y=F[j].re;
	if (j==0) gt.movePen(x,y);
        else gt.drawLine(x,y);
//         maxerr=max(maxerr,abs(y-funcn(x)));
        maxerr=max(maxerr,abs(y-(x>=0)));
	
//      cout << t << " " << z << endl;
      }
      maxerr=max(maxerr,DBL_EPSILON);
      gt.setfgColor("green");
      x=xmin;
      dx=(xmax-xmin)/static_cast<double>(n);
      for (int j=0;j<n;j++,x+=dx) {
//         double y=funcn(x);
	double y=(x>=0);
        gt.drawPlus(x,y,dx*0.25);
      }
      errors[nm]=log(maxerr)/log10;
      errors_min=min(errors_min,errors[nm]);
      errors_max=max(errors_max,errors[nm]);
      gt.flush();
      cout << "errors[" << n << "] = " << maxerr << endl;
      XYGraphTool::WINDOW_TYPE::QuitButton qb;

    	fftw_destroy_plan(plan_forward);
	
	}



    XYGraphTool gt2("errors","log(numPoints)","log(error)",static_cast<double>(nmin)*log(2.)/log10,
      static_cast<double>(nmax)*log(2.)/log10,errors_min,errors_max,
      &cmap, 0,0.5);
    gt2.setbgColor("white");
    gt2.setfgColor("black");
    gt2.drawAxes();
    gt2.setfgColor("blue");
    double t=static_cast<double>(nmin)*log(2.)/log10;
    gt2.movePen(t,errors[nmin]);
    for (int n=nmin+1;n<nmax;n++) {
      t=static_cast<double>(n)*log(2.)/log10;
      gt2.drawLine(t,errors[n]);
    }
    gt2.flush();
    XYGraphTool::WINDOW_TYPE::QuitButton qb;



	delete [] f;
    delete [] fhat;
    delete [] F;
    delete [] FHAT;
    fftw_destroy_plan(plan_backward);
	delete [] errors;

	return 0;
}