#include "include/func.hpp"
#include "include/error.hpp"
#include <gsl/gsl_deriv.h>
#include <cmath>
#include <limits>

namespace kwantix{

//The static variables...
double realfunc::eps = 0;
double realfunc::sqrteps = 0;
double realfunc::root3eps = 0;
bool realfunc::initialised = false;

//Static variables
point gsl_function_wrapper::x;
size_t gsl_function_wrapper::index = 1;
realfunc gsl_function_wrapper::myfunc(0);
gsl_function* gsl_function_wrapper::f = 0;

//******************* Wrapper functions *****************************
gsl_function_wrapper::gsl_function_wrapper( const realfunc &thefunc, 
                                            point p, size_t idx){
  myfunc = thefunc;
  x = p;
  index = idx;
  f -> function = &takemyaddress;
  f -> params = 0;
   
}
  
void gsl_function_wrapper::set_params(const realfunc &thefunc,
                                      point p, size_t idx){
  myfunc = thefunc;
  x = p;
  index = idx;
  f -> function = &takemyaddress;
  f -> params = 0;
}

gsl_function* gsl_function_wrapper::get_gsl_function() const{
  return f;
}

double gsl_function_wrapper::takemyaddress(double xi, void* nothing){
  x(index) = xi;
  nothing = 0;
  return myfunc(x);
}

// **************** realfunc functions *********************************
realfunc::realfunc() : myfunc(0){
  if(!initialised){
    eps = std::numeric_limits<double>::epsilon();
    sqrteps = sqrt(eps);
    root3eps = pow(eps,1/3.0);
    initialised = true;
  }
}

realfunc::realfunc( double(*f)(const point&)) : myfunc(f) {
  if(!initialised){
    eps = std::numeric_limits<double>::epsilon();
    sqrteps = sqrt(eps);
    root3eps = pow(eps,1/3.0);
    initialised = true;
  }
}

void realfunc::set_function_ptr(double (f_in)(const point &p)){
  myfunc = f_in;
}

double realfunc::operator()(const point& p) const{
  return at(p);
}
double realfunc::at(const point& p) const{
  if(myfunc == 0){
    throw no_init(__LINE__,__FILE__);
  }

  return myfunc(p);
}

double realfunc::d(const point& p, size_t k) const{
  gsl_function_wrapper gfw(*this,p,k);
  double result,  abserror;
  double x = p(1); 
  double typx = (1 > log(x) ? 1 : log(x));
  double h = sqrteps*( fabs(x) > typx ? fabs(x) : typx);

  gsl_deriv_central(gfw.get_gsl_function(), x, h, &result, &abserror);
  return result;    
}

double realfunc::d2(const point& p, size_t k1, size_t k2) const{
  //FIXME
  //Figure this out later.
  k1 = k2;
  p.size();
  return 0;
}

kwantix::badArgument
realfunc::no_init(int line, string file) const{
  kwantix::badArgument exc;
  exc.line = line;
  exc.file = file;
  exc.reason = "Did not assign a function pointer to a realfunc object.";
  return exc;
}
}