#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

// 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",0.,50.,0.,150.,&cmap,0,
//      winsize);
//    gt.newPage();
//    gt.setbgColor("white");
//    gt.setfgColor("black");
//    gt.drawAxes();

   XYGraphTool gt2("xExact versus yExact","x","y",0.,50.,0.,150.,&cmap,0,
     winsize);
   gt2.newPage();
   gt2.setbgColor("white");
   gt2.setfgColor("black");
   gt2.drawAxes();

//   double *errors=OPERATOR_NEW_BRACKET(double,number_coarse_steps);
//   double *times=OPERATOR_NEW_BRACKET(double,number_coarse_steps);

  int n=4;
  double dt=(tstop-tstart)/static_cast<double>(number_coarse_steps);
  double h0=dt;
  double t0=tstart;
  double tout=tstart;
  double *yExact=OPERATOR_NEW_BRACKET(double,4);
//   double *yGuess=OPERATOR_NEW_BRACKET(double,4);

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

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

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


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


  double tol = .5;

//for(tol = .1; tol > .000000001; tol /= 10.){

	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);
// 	  tout=tstart+i*dt;
// 	  F77NAME(epsode)(diffun,pederv,n,t0,h0,yGuess,tout,tol,ierror_flag,mf,index);
	  ASSERT(index==0);
/*	cout << "\txGuess(" << tout << ") = " << yGuess[0] << "\tyGuess(" << tout << ") = " << yGuess[2] << endl;*/
	cout << "\txExact(" << tout << ") = " << yExact[0] << "\tyExact(" << tout << ") = " << 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]);
	     if (i==0) {
// 		gt.movePen(yGuess[0],yGuess[2]);
		gt2.movePen(yExact[0],yExact[2]);
		}
	     else {
// 		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,eps,ierror_flag,mf,index);
// 	tout=tstop;
// 	F77NAME(epsode)(diffun,pederv,n,t0,h0,yGuess,tout,tol,ierror_flag,mf,index);
// 	cout << "\txGuess(" << tout << ") = " << yGuess[0] << "\tyGuess(" << tout << ") = " << yGuess[2] << endl;
	cout << "\txExact(" << tout << ") = " << yExact[0] << "\tyExact(" << tout << ") = " << 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;
	  }


//}



//   times[number_coarse_steps-1]=t0;
//   errors[number_coarse_steps-1]=log(abs(y0[0]-exact(t0)))/log10;
//   emax=max(emax,errors[number_coarse_steps-1]);
//   emin=min(emin,errors[number_coarse_steps-1]);
//   gt.drawLine(t0,y0[0]);
//   gt.flush();
//   {
//   XYGraphTool::WINDOW_TYPE::QuitButton qb;
//   }
// 
//   XYGraphTool gte("error versus t","t","log_10(error)",tstart,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;
}