#ifndef SHRIMPS_Tools_Special_Functions_H
#define SHRIMPS_Tools_Special_Functions_H

#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Math/MathTools.H"
#include "ATOOLS/Math/Gauss_Integrator.H"

using namespace ATOOLS;

namespace SHRIMPS {
  /*!
    \file
    This file contains an assortment of special functions
    neccessary for the SHRIMPS::Form_Factor.  
    
    \todo Further expand and move to ATOOLS::MathTools
  */
  class Special_Functions {
  public:
    double LnGamma(const double & x) const
    {
      static const double coeff[6] = 
	{ 76.18009172947146, -86.50532032941677, 24.01409824083091, 
	  -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5 };
      double y(x), tmp(x+5.5);
      tmp -= (x+0.5)*log(tmp);
      double arg(1.000000000190015);
      for (short int j=0;j<6;j++) arg += coeff[j]/++y;
      return -tmp+log(2.5066282746310005*arg/x);
    }

    double IncompleteGamma(const double & a, const double & x) const
    {
      double incgam(0.);
      if (x<0. || a<0.) {
	msg_Error()<<"Error in "<<METHOD<<"("<<a<<", "<<x<<"):\n"
		   <<"   Out of bounds. "
		   <<"Will return 0 and hope for the best.\n";
	return incgam;
      }
      if (a==0.) {
	incgam = -0.577215664902-log(x);
	double pref(1.);
	for (int i=1;i<21;i++) {
	  incgam += pref*pow(x,i);
	  pref *= -double(i)/sqr(double(i+1));
	}
	return incgam;
      }
      double lngamma(LnGamma(a));
      if (x<a+1.) {
	if (x==0.) return incgam;
	else {
	  double aprime(a), del(1./a), sum(1./a);
	  for (short int i=0;i<100;i++) {
	    ++aprime;
	    del *= x/aprime;
	    sum += del;
	    if (dabs(del)<dabs(sum)*1.e-12) {
	      incgam = 1.-sum*exp(-x+a*log(x)-lngamma);
	      break;
	    }
	  }
	}
      }
      else {
	msg_Error()<<"Error in "<<METHOD<<"("<<a<<", "<<x<<") :\n"
		   <<"   Not implemented yet. "
		   <<"Will return 0 and hope for the best.\n";
      }
      return incgam;
    }

    double Jn(const int order,const double & arg) const 
    {
      double jn(0.),darg(dabs(arg));
      if (order==0) {
	if (darg<=1.e-12) return 1.;
	if (darg<8.) {
	  double x2(darg*darg);
	  double num = 
	    ((((-184.9052456*x2+77392.33017)*x2-
	       11214424.18)*x2+651619640.7)*x2-13362590354.)*x2+57568490574.;
	  double den =
	    ((((x2+267.8532712)*x2+59272.64853)*
	      x2+9494680.718)*x2+1029532985)*x2+57568490411.;
	  jn = num/den;
	}
	else {
	  double x2(64.0/sqr(darg)), xx(darg-0.785398164);
	  double pref1 = 
	    (((0.2993887211e-6*x2-0.2073370639e-5)*x2+
	      0.2734510407e-4)*x2-0.1098628627e-2)*x2+1.;
	  double pref2 =
	    (((-0.934945152e-7*x2+0.7621095161e-6)*x2-
	      0.6911147651e-5)*x2+0.1430488765e-3)*x2-0.1562499995e-1;
	  jn = sqrt(0.636619772/darg)*(cos(xx)*pref1-8./darg*sin(xx)*pref2);
	}
      }
      else {
	msg_Error()<<"Error in "<<METHOD<<"("<<order<<", "<<arg<<"):\n"
		   <<"   Not implemented yet.  Exit the run.\n";
	exit(1);
      }
      return jn;
    }
    
    double In(const int order,const double & arg) const 
    {
      double in(0.),darg(dabs(arg)), y;
      switch (order) {
      case 1:
	if (darg<3.75) {
	  y  = sqr(arg/3.75);
	  in = darg*((((((0.32411e-3*y+0.301532e-2)*y+0.2658733e-1)*y+
			0.15084934)*y+0.51498869)*y+0.87890594)*y+0.5);
	}
	else {
	  y  = 3.75/darg;
	  in = ((((((((-0.420059e-2*y+0.1787654e-1)*y-0.2895312e-1)*y+
		     0.2282967)*y-0.1031555e-1)*y+0.163801e-2)*y-
		  0.362018e-2)*y-0.3988024e-1)*y+0.39894228) * 
	    exp(darg)/sqrt(darg);
	}
	break;
      case 0:
	if (darg<3.75) {
	  y  = sqr(arg/3.75);
	  in = (((((0.45813e-2*y+0.360768e-1)*y+0.2659723)*y+1.2067492)*y
		 +3.0899424)*y+3.5156229)*y+1.0;
	}
	else {
	  y  = 3.75/darg;
	  in = ((((((((0.392377e-2*y-0.1647633e-1)*y+0.2635537e-1)*y-
		     0.2057706)*y+0.916281e-2)*y-0.157565e-2)*y+
		  0.225319e-2)*y+0.1328592e-1)*y+0.39894228) * 
	    exp(darg)/sqrt(darg);
	}    
	break;
      default:
	msg_Error()<<"Error in "<<METHOD<<"("<<order<<", "<<arg<<"):\n"
		   <<"   Not implemented yet.  Exit the run.\n";
	exit(1);
      }
      return in;
    }


    double I0(const double & arg)  const { return In(0,arg); }

    double Kn(const int order,const double & arg) const
    {
      if (arg<=0.) return 0.;
      double kn(0.),darg(dabs(arg)), y;
      if (order==1) {
	if (darg<=2.0) {
	  y   = sqr(darg)/4.;
	  kn  = log(darg/2.)*In(1,darg);
	  kn += ((((((-0.4686e-4*y-0.110404e-2)*y-0.1919402e-1)*y-0.18156897)*y
		   -0.67278579)*y+0.15443144)*y+1.)/darg;
	}
	else {
	  y   = 2./darg;
	  kn  = exp(-darg)/sqrt(darg);
	  kn *= ((((((-0.68245e-3*y+0.325614e-2)*y-0.780353e-2)*y+
		    0.1504268e-1)*y-0.3655620e-1)*y+0.23498619)*y+1.25331414);
	}
      }
      else {
	msg_Error()<<"Error in "<<METHOD<<"("<<order<<", "<<arg<<")\n:"
		   <<"   Not implemented yet.  Exit the run.\n";
	exit(1);
      }
      return kn;
    }

    void TestBessel(const std::string & dirname) const
    {
      double arg(0.);
      std::ofstream was;
      std::string filename = dirname+std::string("/BesselJ0I0I1K1.dat");
      was.open(filename.c_str());
      long int cnt(0);
      while (arg<10.) {
	was<<arg<<"   "<<Jn(0,arg)<<"   "
	   <<In(0,arg)<<"   "<<In(1,arg)<<"   "<<Kn(1,arg)<<std::endl;
	//if (cnt==0 || !(cnt%500)) 
	//  msg_Out()<<"  J_0("<<arg<<") = "<<Jn(0,arg)
	//	   <<"  K_1("<<arg<<") = "<<Kn(1,arg)<<std::endl;
	arg += 0.001;
	cnt++;
      }
      was.close();
    }
  };

  extern Special_Functions SF;


  class ExpInt : public ATOOLS::Function_Base {
  public:
    virtual double GetValue(double T) { return -exp(-T)/T; }
    virtual double operator()(double T) { return GetValue(T); }
    virtual double GetExpInt(double X){
	if(X>0.) exit(1);
  	ATOOLS::Gauss_Integrator integrator(this);
  	return integrator.Integrate(-X,100000.,1.e-2,1);
        }
  };

}
#endif