#include <cmath>
#include <iostream>
#include <math.h> // for HUGE_VAL,M_PI
#include <rfftw.h>
#include <stdlib.h>

#include "XYGraphTool.H"

#ifdef USE_GTK
#include <gtk/gtkmain.h>
#include "GTKGUI.H"
#include "GTKGUIVirtualInput.H"
#include "GTKWindow.H"
#else
#include "GUI.H"
#include "GUIVirtualInput.H"
#endif

#include <time.h>





// changes: reformulate exact to follow the same pattern as the generated data, that is, store a 
// vector of "exact" following how the regular stuff is generated, except with an error tolerance
// of DBL_EPSILON or so.  Then we also need to change the formulae for the eqns we're solving to
// match the formulae given in the assignment.


enum HOW_STIFF{STIFF,NOT_STIFF};
HOW_STIFF how_stiff=NOT_STIFF;

enum ITERATIVE_CORRECTION_METHOD{FUNCTIONAL_ITERATION,EXACT_JACOBIAN,
  FINITE_DIFFERENCE_JACOBIAN,DIAGONAL_JACOBIAN};
ITERATIVE_CORRECTION_METHOD iterative_correction_method=
  FUNCTIONAL_ITERATION;

enum ERROR_FLAG{CONTROL_ASOLUTE_ERROR,CONTROL_RELATIVE_TO_CURRENT,
  CONTROL_RELATIVE_TO_MAX};
ERROR_FLAG error_flag=CONTROL_ASOLUTE_ERROR;
int ierror_flag=static_cast<int>(error_flag);

char *display_name=0;
double winsize=0.5;

double tstart=0.;
double tstop=20.;

double eps=sqrt(DBL_EPSILON);


int number_coarse_steps=20;

extern "C" {
  void F77NAME(epsode)(
    void (*diffun)(const int &n,const double &t,const double *y,
                   double *ydot),
    void (*pederv)(const int &n,const double &t,const double *y,
                   double *pd,const int &n0),
    const int &n,double &t0,double &h0,double *y0,const double &tout,
    const double &eps,const int &ierror,const int &mf,int &index);
}

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// variables are in order: y[0] = x, y[1] = z, y[2] = y, y[3] = q

void diffun(const int &n,const double &t,const double *y,double *ydot) { 

 double R = sqrt(y[0]*y[0]+y[2]*y[2]);

  ydot[0] = y[1];
  ydot[1] = -y[0]/(R*R*R);
  ydot[2] = y[3];
  ydot[3] = -y[2]/(R*R*R);
}
void pederv(const int &n,const double &t,const double *y,double *pd,
const int &n0) { 

  double R = sqrt(y[0]*y[0]+y[2]*y[2]);
  pd[0] = 0.;
  pd[1] = 1.;
  pd[2] = 0.;
  pd[3] = 0.;
  pd[4] = -1./(R*R*R) + 3. * y[0]*y[0]/(R*R*R*R*R);
  pd[5] = 0.;
  pd[6] = 3. * y[0]*y[2]/(R*R*R*R*R);
  pd[7] = 0.;
  pd[8] = 0.;
  pd[9] = 0.;
  pd[10] = 0.;
  pd[11] = 1.;
  pd[12] = 3. * y[0]*y[2]/(R*R*R*R*R);
  pd[13] = 0.;
  pd[14] = -1./(R*R*R) + 3. * y[2]*y[2]/(R*R*R*R*R);
  pd[15] = 0.;
}
// double exact(const double &t) { return y_init*exp(rate*t); }
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

int main(int argc, char* argv[]){

#ifdef USE_GTK
   GTKWindow::gtkInit(argc,argv);
#endif

   double emax=-HUGE_VAL;
   double emin=HUGE_VAL;
//   double ymax=max(exact(tstart),exact(tstop));
//   double ymin=min(exact(tstart),exact(tstop));
  double log10=log(10.);

   Palette pal;
  XYGraphTool::WINDOW_TYPE::COLOR_MAP_TYPE cmap(&pal);
   XYGraphTool gt("xGuess versus yGuess","x","y",-3.,3.,-2.,2.,&cmap,0,
     winsize);
  clock_t start,end;

//    XYGraphTool gt2("xExact versus yExact","x","y",-3.,3.,-2.,2.,&cmap,0,
//      winsize);
//    gt2.newPage();
//    gt2.setbgColor("white");
//    gt2.setfgColor("black");
//    gt2.drawAxes();

   double *errors=OPERATOR_NEW_BRACKET(double,9);
   double *times=OPERATOR_NEW_BRACKET(double,9);

  int n=4;

  double *yExact=OPERATOR_NEW_BRACKET(double,n);
   double *yGuess=OPERATOR_NEW_BRACKET(double,n);

	// initialization
	// variables are in order: x, z, y, q

double tol;
int counter = 0;

  double dt=(tstop-tstart)/static_cast<double>(number_coarse_steps);
  double h0=dt;
  double h1=dt;
  double t0=tstart;
  double t1=tstart;
  double tout=tstart;

  int mf=10*static_cast<int>(how_stiff)
        +static_cast<int>(iterative_correction_method);


  yExact[0]= .1;
  yExact[1] = 0;
  yExact[2] = 0;
  yExact[3] = sqrt(19);
  
  int index=1;

	for (int i=1;i<number_coarse_steps;i++) {

	  tout=tstart+i*dt;
	  F77NAME(epsode)(diffun,pederv,n,t0,h0,yExact,tout,eps,ierror_flag,mf,index);
	}


for(tol = .01; tol > .0000000001; tol /= 10.){


   gt.newPage();
   gt.setbgColor("white");
   gt.setfgColor("black");
   gt.drawAxes();

  dt=(tstop-tstart)/static_cast<double>(number_coarse_steps);
  h0=dt;
  h1=dt;
  t0=tstart;
  t1=tstart;
  tout=tstart;

   yGuess[0]= .1;
   yGuess[1] = 0;
   yGuess[2] = 0;
   yGuess[3] = sqrt(19);


  index = 1;


  gt.movePen(yGuess[0],yGuess[2]);
//   gt2.movePen(yExact[0],yExact[2]);

  start = clock();

	for (int i=1;i<number_coarse_steps;i++) {

	  tout=tstart+i*dt;
// 	  F77NAME(epsode)(diffun,pederv,n,t0,h0,yExact,tout,tol,ierror_flag,mf,index);
//  	  tout=tstart+i*dt;
 	  F77NAME(epsode)(diffun,pederv,n,t1,h1,yGuess,tout,tol,ierror_flag,mf,index);
	  ASSERT(index==0);
// 	cout << "\txGuess(" << t1 << ") = " << yGuess[0] << "\tyGuess(" << t1 << ") = " << yGuess[2] << endl;
// 	cout << "\txExact(" << t0 << ") = " << yExact[0] << "\tyExact(" << t0 << ") = " << yExact[2] << endl;
	//     times[i]=t0;
	//     errors[i]=log(abs(y0[0]-exact(t0)))/log10;
	//     emax=max(emax,errors[i]);
	//     emin=min(emin,errors[i]);
	     
		gt.drawLine(yGuess[0],yGuess[2]);
// 		gt2.drawLine(yExact[0],yExact[2]);
		
		
	}
	index=2;
	tout=tstop;
// 	F77NAME(epsode)(diffun,pederv,n,t0,h0,yExact,tout,tol,ierror_flag,mf,index);
	tout=tstop;
	F77NAME(epsode)(diffun,pederv,n,t1,h1,yGuess,tout,tol,ierror_flag,mf,index);

  	end = clock();

	cout << "tol: " << tol << "\txGuess(" << t1 << ") = " << yGuess[0] << "\tyGuess(" << t1 << ") = " << yGuess[2] << endl;
	errors[counter] = yGuess[0];
	times[counter] = (100*(end-start))/(double)(CLOCKS_PER_SEC);

	cout << endl << "error: " << log(abs(errors[counter]-yExact[0]))/log10 << "\ttime: " << times[counter] << endl;

	counter++;
// 	cout << "\txExact(" << t0 << ") = " << yExact[0] << "\tyExact(" << t0 << ") = " << yExact[2] << endl;

	

	ASSERT(index==0);

	gt.drawLine(yGuess[0],yGuess[2]);
// 	gt2.drawLine(yExact[0],yExact[2]);
	  gt.flush();
// 	gt2.flush();
	  {
	  XYGraphTool::WINDOW_TYPE::QuitButton qb;
	  }


}

	for(int j = 0; j < 9; j++){
		errors[j] = log(abs(errors[j]-yExact[0]))/log10;
	}


//   times[number_coarse_steps-1]=t0;
//   errors[number_coarse_steps-1]=log(abs(y0[0]-exact(t0)))/log10;
   emax=max(emax,errors[0]);
   emin=min(emin,errors[8]);
//   gt.drawLine(t0,y0[0]);
//   gt.flush();
//   {
//   XYGraphTool::WINDOW_TYPE::QuitButton qb;
//   }
// 
  XYGraphTool gte("error versus computer time","computer time","log_10(error)",0.,tstop,
    emin,emax,&cmap,0,winsize);
  gte.setbgColor("white");
  gte.setfgColor("black");
  gte.drawAxes();
  gte.movePen(times[0],errors[0]);
  for (int i=1;i<number_coarse_steps;i++) {
    gte.drawLine(times[i],errors[i]);
  }
  gte.flush();
  {
  XYGraphTool::WINDOW_TYPE::QuitButton qb;
  }

  delete [] yExact;
   delete [] yGuess;
//   delete [] times;
//   delete [] errors;

  if (display_name) delete display_name; display_name=0;
}