/**
 * Copyright Sören Molander, 2000-2004
 * The code maybe copied and modified at will  
 * 
 * Description: mathematical functions
 *
 * Author: C.A. Bertulani http://www.if.ufrj.br/persons/bertulani.html
 *         Sören Molander
 */


/**
* \file
* Mathematical functions and definitions.  
*
*
**/    

#ifndef INCLUDE_MATHUTIL_H
#define INCLUDE_MATHUTIL_H

#include <math.h>
#include <stdlib.h>
#include <vector>

/** Utilities for mathematical functions and definitions */
namespace Mathutil {


  /** A random number class (see definition .cpp file for credits) */
  class RandomNumber {
  private: 
    //   float* u;
    float  c, cd, cm;
    int i97, j97;
    bool initialized;
    std::vector<float> u;
    int seed1, seed2;
  public: 
    RandomNumber() {}
    void init();
    float random();
  private: 
    int rmarin (int ij, int kl);
    int ranmar(float rvec[], int len);
    void sow(int& seed1, int& seed2);
  };
  

  /** Constants */
  const double Pi             = 3.14159265358979323846;
  const double OneOverSqrt2Pi = 0.39894228040143270286;
  const double E              = 2.71828182845904523546;
  const double SqrtTwo        = 1.41421356237309504880;

  /** Utility functions */

  template<class T>  T Sign(T a, T b) {return (b >= 0.0 ? fabs(a) : -fabs(a));}
  template<class T>  int Sign(T x) {return x>=0 ? 1 :-1;} 
  template<class T>  T Min(T a,T b) {return a<b ? a : b;} 
  template<class T>  T Max(T a,T b) {return a>b ? a : b;} 
  template<class T>  T Min(T a,T b,T c) {return Min(a,Min(b,c));} 
  template<class T>  T Max(T a,T b,T c) {return Max(a,Max(b,c));} 
  template<class T>  T sqr (T a) {return T(a*a);} 
/*   template<class T>  T sqrt (T a) {return T(sqrt(float(a)));} */
  template<class T>  T abs(T a) {return a<0 ? -a : a;}

  inline float logB(double base, double x) {
    float out = float(fabs(x) > 0.00000000001 ? log(x)/log(base) : 0);
    return out;
  }
  
  inline double gaussian(double x) {return OneOverSqrt2Pi*exp(-0.5*x*x);}  

  inline double gaussian(double x, double s) {return 1.0/sqrt(s)*OneOverSqrt2Pi*exp(-0.5*x*x/s);}  

  // N-th order derivatives
  inline double dGaussianPol(double x, double s, int order) {
    switch (order) {
    case 0: return 1;
    case 1: return -x/s;
    case 2: return 1/s*(-1 + x*x/s);
    default: return -1 / s * (order-1)* dGaussianPol(x,s,order-2) + 
	       x * dGaussianPol(x,s,order-1);
    }
  }


  
  inline int fac(int n) {
    int f=1;
    n = (n <= 0 ? 1 : n);
    for (int k=1;k<=n;k++) 
      f *= k;
    return f;
  } 

  /** Random numbers */

  /** Uniformly distributed numbers in the range [Min..Max] by linear congruences */
  double UniformDistribution(double min, double max);

  /** Normally distributed variables with mean and stdv */
  double GaussianDistribution(double mean, double stdv);

   /** Bessel functions */

  double bessj0(double x);
  void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi);
  double bessi(int n, double x);
  double bessi0(double x);
  double bessi1(double x);
  double bessj(int n, double x);
  double bessj0(double x);
  double bessj1(double x);
  void bessjy(double x, double xnu, double *rj, double *ry, double *rjp, double *ryp);
  double bessk(int n, double x);
  double bessk0(double x);
  double bessk1(double x);
  double bessy(int n, double x);
  double bessy0(double x);
  double bessy1(double x);
  double chebev(double a, double b, double c[], int m, double x);
  void sphbes(int n, double x, double *sj, double *sy, double *sjp, double *syp);
  double scaleSpaceGaussian(double x,double t);
  double errorFunctionPhi(double x);
}


#endif