/**
 * 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
 */


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


#include <math.h>
#include <iostream>
#include <fstream>
#include <ctype.h>
#include <stdlib.h>
#include <stdio.h>
#include <vector>
#include <assert.h>
#include <time.h>
#include <string>
#include "mathutil.h"

namespace Mathutil {
 using namespace std;

 /** Utility functions */


  /************************************************************************
   * This random number generator originally appeared in "Toward a Universal
   *Random Number Generator" by George Marsaglia and Arif Zaman.
   * Florida State University Report: FSU-SCRI-87-50 (1987)
   *
   *It was later modified by F. James and published in "A Review of Pseudo-
   * random Number Generators"
   *
   *Converted from FORTRAN to C by Phil Linttell, James F. Hickling
   *Management Consultants Ltd, Aug. 14, 1989.
   *
   *THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
   *      (However, a newly discovered technique can yield
   *        a period of 10^600. But that is still in the development stage.)
   *
   *It passes ALL of the tests for random number generators and has a period
   *  of 2^144, is completely portable (gives bit identical results on all
   *  machines with at least 24-bit mantissas in the floating point
   *  representation).
   *
   *The algorithm is a combination of a Fibonacci sequence (with lags of 97
   *  and 33, and operation "subtraction plus one, modulo one") and an
   *  "arithmetic sequence" (using subtraction).
   *
   * On a Vax 11/780, this random number generator can produce a number in
   *   13 microseconds.
   *
   *
   * Ported to C++ by Sören Molander, 2005-05-07
   */

// using namespace std;


  float RandomNumber::random() {
    float vec[1];
    int len = 1;
    ranmar(vec,len);
    return vec[0];
  }  

  void RandomNumber::init() {
    u.resize(97);
    sow(seed1,seed2);
    rmarin(seed1,seed2);
  }
  
  /************************************************************************
   * This is the initialization routine for the random number generator RANMAR()
   * NOTE: The seed variables can have values between:    0 <= IJ <= 31328
   *                                                    0 <= KL <= 30081
   * The random number sequences created by these two seeds are of sufficient
   * length to complete an entire calculation with. For example, if several
   * different groups are working on different parts of the same calculation,
   * each group could be assigned its own IJ seed. This would leave each group
   * with 30000 choices for the second seed. That is to say, this random
   * number generator can create 900 million different subsequences -- with
   * each subsequence having a length of approximately 10^30.
   *
   * Use IJ = 1802 & KL = 9373 to test the random number generator. The
   * subroutine RANMAR should be used to generate 20000 random numbers.
   * Then display the next six random numbers generated multiplied by 4096*4096
   * If the random number generator is working properly, the random numbers
   * should be:
   *       6533892.0  14220222.0   7275067.0
   *        6172232.0   8354498.0  10633180.0
   */

  int RandomNumber::rmarin(int ij, int kl) {
    
    float s, t;
    int i, j, k, l, m;
    int ii, jj;
    
    /* Change FALSE to TRUE in the next statement to test the
       random routine.*/
    
    initialized = true;
    
    if ( ( ij < 0 || ij > 31328 ) ||
         ( kl < 0 || kl > 30081 ) )
      {
        printf ("RMARIN: The first random number seed must have a "
                "value between 0 and 31328\n");
        printf ("        The second random number seed must have a "
                "value between 0 and 30081");
        return 1;
      }
    
    i = (int)fmod(ij/177.0, 177.0) + 2;
    j = (int)fmod(ij      , 177.0) + 2;
    k = (int)fmod(kl/169.0, 178.0) + 1;
    l = (int)fmod(kl      , 169.0);
    
    for ( ii=0; ii<=96; ii++ )
      {
        s = (float)0.0;
        t = (float)0.5;
        for ( jj=0; jj<=23; jj++ )
          {
            m = (int)fmod( fmod(i*j,179.0)*k , 179.0 );
            i = j;
            j = k;
            k = m;
            l = (int)fmod( 53.0*l+1.0 , 169.0 );
            if ( fmod(l*m,64.0) >= 32)
              s = s + t;
            t = (float)(0.5 * t);
          }
        u[ii] = s;
      }
    
    c  = (float)(  362436.0 / 16777216.0);
    cd = (float)( 7654321.0 / 16777216.0);
    cm = (float)(16777213.0 / 16777216.0);
    
    i97 = 96;
    j97 = 32;
    
    initialized = true;
    
    return 0;
  }
  
  int RandomNumber::ranmar(float rvec[], int len) {
    float uni;
    int ivec;
    
    if ( !initialized )
      {
        printf ("RANMAR: Call the initialization routine (RMARIN) "
                "before calling RANMAR.\n");
        return 1;
      }
    
    for ( ivec=0; ivec < len; ivec++)
      {
        uni = u[i97] - u[j97];
        if ( uni < 0.0F )
          uni = uni + 1.0;
        u[i97] = uni;
        i97--;
        if ( i97 < 0 )
          i97 = 96;
        j97--;
        if ( j97 < 0 )
          j97 = 96;
        c = c - cd;
        if ( c < 0.0F )
          c = c + cm;
        uni = uni - c;
        if ( uni < 0.0F )
          uni = uni + 1.0;
        rvec[ivec] = uni;
      }
    return 0;
  }
  
  /* I use the following procedure in TC to generate seeds:
   *     
   * The sow() procedure calculates two seeds for use with the random number
   * generator from the system clock.  I decided how to do this myself, and
   * I am sure that there must be better ways to select seeds; hopefully,
   * however, this is good enough.  The first seed is calculated from the values
   * for second, minute, hour, and year-day; weighted with the second most
   * significant and year-day least significant.  The second seed weights the
   * values in reverse.
   */
  
  void RandomNumber::sow(int& seed1, int& seed2) {
    struct tm *tm_now;
    float s_sig, s_insig, maxs_sig, maxs_insig;
    long secs_now;
    int s, m, h, d, s1, s2;
    
    time(&secs_now);
    tm_now = localtime(&secs_now);
    
    s = tm_now->tm_sec + 1;
    m = tm_now->tm_min + 1;
    h = tm_now->tm_hour + 1;
    d = tm_now->tm_yday + 1;
    
    
    maxs_sig   = (float)(60.0 + 60.0/60.0 + 24.0/60.0/60.0 +
                         366.0/24.0/60.0/60.0);
    maxs_insig = (float)(60.0 + 60.0*60.0 + 24.0*60.0*60.0 +
                         366.0*24.0*60.0*60.0);
    
    
    s_sig      = (float)(s + m/60.0 + h/60.0/60.0 + d/24.0/60.0/60.0);
    s_insig    = (float)(s + m*60.0 + h*60.0*60.0 + d*24.0*60.0*60.0);
    
    
    s1 = (int)(s_sig   / maxs_sig   * 31328.0);
    s2 = (int)(s_insig / maxs_insig * 30081.0);
    
    seed1 = s1;
    seed2 = s2;
  }




/*********************************************************************
   Returns the Bessel functions All of them. Pick up our choice.
   C.A. Bertulani        May/16/2000
*********************************************************************/

/*********************************************************************
   Returns the Bessel function J_0(x) for any real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
  double bessj0(double x) {
    double ax,z;
    double xx,y,ans,ans1,ans2;		/* Acumulate polynomials in double */
    /* precision. */
    
    if ((ax=fabs(x)) < 8.0) {		/* Direct rational function fit. */
      y=x*x;
      ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
					      +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
      ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
					    +y*(59272.64853+y*(267.8532712+y*1.0))));
      ans=ans1/ans2;
    } else {		/* Fitting function. */
      z=8.0/ax;
      y=z*z;
      xx=ax-0.785398164;
      ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
				      +y*(-0.2073370639e-5+y*0.2093887211e-6)));
      ans2 = -0.1562499995e-1+y*(0.1430488765e-3
				 +y*(-0.6911147651e-5+y*(0.7621095161e-6
							 -y*0.934935152e-7)));
      ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
    }
    return ans;
  }
  /*********************************************************************
 Reurns the Bessel function J_1(x) for any real x.
C.A. Bertulani        May/16/2000
  *********************************************************************/
  double bessj1(double x)  {
    double ax,z;
    double xx,y,ans,ans1,ans2;		/* Accumulate polynomials in double */
    /* precision. */
    
    
    if ((ax=fabs(x)) < 8.0) {		/* Rational function approximation. */
      y=x*x;
      ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
						+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
      ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
					     +y*(99447.43394+y*(376.9991397+y*1.0))));
      ans=ans1/ans2;
    } else {
      z=8.0/ax;
      y=z*z;
      xx=ax-2.356194491;
      ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
				 +y*(0.2457520174e-5+y*(-0.240337019e-6))));
      ans2=0.04687499995+y*(-0.2004690873e-3
			    +y*(0.8449199096e-5+y*(-0.88228987e-6
						   +y*0.105787412e-6)));
      ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
      if (x < 0.0) ans = -ans;
    }
    return ans;
  }
/*********************************************************************
   Returns the Bessel function J_n(x) for any real x and n>=2.
   C.A. Bertulani        May/16/2000
*********************************************************************/
// #define ACC 40.0		/* Make larger to increase accuracy. */
// #define BIGNO 1.0e10
// #define BIGNI 1.0e-10
  
  double bessj(int n, double x) {
    double ACC = 40.0;
    double BIGNO = 1.0e10;
    double BIGNI = 1.0e-10;
    
//     double bessj0(double x);
//     double bessj1(double x);
    int j,jsum,m;
    double ax,bj,bjm,bjp,sum,tox,ans;
    
    if (n < 2) std::cerr<< "Index n less than 2 in bessj\n";
    ax=fabs(x);
    if (ax == 0.0)
      return 0.0;
    else if (ax > (double) n) {		/* Upwards recurrence from J_0 and J_1 */
      tox=2.0/ax;
      bjm=bessj0(ax);
      bj=bessj1(ax);
      for (j=1;j<n;j++) {
	bjp=j*tox*bj-bjm;
	bjm=bj;
	bj=bjp;
      }
      ans=bj;
    } else {		/* Downwards recurrence from an even m here computed. */
      tox=2.0/ax;
      m=2*((n+(int) sqrt(ACC*n))/2);
      jsum=0;		/* jsum will alternate between 0 and 1; when it is 1, we */
      bjp=ans=sum=0.0;		/* accumulate in sum the even terms. */
      bj=1.0;
      for (j=m;j>0;j--) {		/* Downward recurrence. */
	bjm=j*tox*bj-bjp;
	bjp=bj;
	bj=bjm;
	if (fabs(bj) > BIGNO) {		/* Renormalize to prevent overflows. */
	  bj *= BIGNI;
	  bjp *= BIGNI;
	  ans *= BIGNI;
	  sum *= BIGNI;
	}
	if (jsum) sum += bj;		/* Accumulate the sum.  */
	jsum=!jsum;					/* Change 0 to 1 or vice versa. */
	if (j == n) ans=bjp;		/* Save the unnormalized answer. */
      }
      sum=2.0*sum-bj;
      ans /= sum;						/* Normalize the answer. */
    }
    return x < 0.0 && (n & 1) ? -ans : ans;
  }
/*********************************************************************
   Returns the Bessel function I_0(x) for any real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
// #undef ACC
// #undef BIGNO
// #undef BIGNI
  double bessi0(double x) {
    double ax,ans;
    double y;				/* Accumulate polynomials in double precision.  */
    
    if ((ax=fabs(x)) < 3.75) {	/* Polynomial fit. */	
      y=x/3.75;
      y*=y;
      ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
					   +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
    } else {
      y=3.75/ax;
      ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
					    +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
									       +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
														    +y*0.392377e-2))))))));
    }
    return ans;
}
/*********************************************************************
   Returns the Bessel function K_0(x) for any positive real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
  double bessk0(double x) {
    double bessi0(double x);
    double y,ans;			/* Accumulate polynomials in double precision.  */
    
    if (x <= 2.0) {	/* Polynomial fit. */	
      y=x*x/4.0;
      ans=(-log(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420
						  +y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
										    +y*(0.10750e-3+y*0.74e-5))))));
    } else {
      y=2.0/x;
      ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
					   +y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
										+y*(-0.251540e-2+y*0.53208e-3))))));
    }
    return ans;
  }
/*********************************************************************
   Returns the Bessel function I_1(x) for any real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
  double bessi1(double x) {
    double ax,ans;
    double y;			/* Accumulate polynomials in double precision.  */
    
    if ((ax=fabs(x)) < 3.75) {	/* Polynomial fit. */	
      y=x/3.75;
      y*=y;
      ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
						 +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
    } else {
      y=3.75/ax;
      ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
					   -y*0.420059e-2));
      ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
					 +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
      ans *= (exp(ax)/sqrt(ax));
    }
    return x < 0.0 ? -ans : ans;
  }
/*********************************************************************
   Returns the Bessel function K_1(x) for any positive real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
  double bessk1(double x) {
    double bessi1(double x);
    double y,ans;			/* Accumulate polynomials in double precision.  */
    
    if (x <= 2.0) {	/* Polynomial fit. */	
      y=x*x/4.0;
      ans=(log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144
						 +y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1
										   +y*(-0.110404e-2+y*(-0.4686e-4)))))));
    } else {
      y=2.0/x;
      ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619
					   +y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2
										+y*(0.325614e-2+y*(-0.68245e-3)))))));
    }
    return ans;
  }
  /*********************************************************************
  Returns the Bessel function K_n(x) for any positive real x and n>=2.
  C.A. Bertulani        May/16/2000
  *********************************************************************/
  double bessk(int n, double x) {
    double bessk0(double x);
    double bessk1(double x);
    int j;
    double bk,bkm,bkp,tox;
    
    if (n < 2) std::cerr << "Index n less than 2 in bessk";
    tox=2.0/x;
    bkm=bessk0(x);		/* Upward recurrence for all x ...   */
    bk=bessk1(x);
    for (j=1;j<n;j++) {		/* ... and there it is.   */
      bkp=bkm+j*tox*bk;
      bkm=bk;
      bk=bkp;
    }
    return bk;
  }
/*********************************************************************
   Returns the Bessel function I_n(x) for any real x and n>=2.
   C.A. Bertulani        May/16/2000
*********************************************************************/
//#define ACC 40.0		/* Make larger to increase accuracy.  */
//#define BIGNO 1.0e10
//#define BIGNI 1.0e-10
  
  double bessi(int n, double x) {
    double ACC = 40.0;
    double BIGNO = 1.0E10;
    double BIGNI = 1.0e-10;

    double bessi0(double x);
    int j;
    double bi,bim,bip,tox,ans;
    
    if (n < 2) std::cerr << "Index n less than 2 in bessi\n";
    if (x == 0.0)
      return 0.0;
    else {
      tox=2.0/fabs(x);
      bip=ans=0.0;
      bi=1.0;
      for (j=2*(n+(int) sqrt(ACC*n));j>0;j--) {  /* Download recurrence from even m.  */
	bim=bip+j*tox*bi;
	bip=bi;
	bi=bim;
	if (fabs(bi) > BIGNO) {		/* Renormalize to prevent overflows. */
	  ans *= BIGNI;
	  bi *= BIGNI;
	  bip *= BIGNI;
	}
	if (j == n) ans=bip;
      }
      ans *= bessi0(x)/bi;		/* Normalize with bessi0.  */
//       std::cerr << """bessi("<<n<<","<<x<<") = "<<ans<<endl;
      return x < 0.0 && (n & 1) ? -ans : ans;
    }
  }
  //#undef ACC
  //#undef BIGNO
  //#undef BIGNI
/*********************************************************************
   Returns the Bessel function Y_0(x) for any positive real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
  double bessy0(double x)  {
    double bessj0(double x);
    double z;
    double xx,y,ans,ans1,ans2;	/* Accumulate polynomials in double precision.  */
    
    if (x < 8.0) {		/* Direct rational function fit. */
      y=x*x;
      ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
					      +y*(10879881.29+y*(-86327.92757+y*228.4622733))));
      ans2=40076544269.0+y*(745249964.8+y*(7189466.438
					   +y*(47447.26470+y*(226.1030244+y*1.0))));
      ans=(ans1/ans2)+0.636619772*bessj0(x)*log(x);
    } else {		/* Fitting function. */
      z=8.0/x;
      y=z*z;
      xx=x-0.785398164;
      ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
				      +y*(-0.2073370639e-5+y*0.2093887211e-6)));
      ans2 = -0.1562499995e-1+y*(0.1430488765e-3
				 +y*(-0.6911147651e-5+y*(0.7621095161e-6
							 +y*(-0.934945152e-7))));
      ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
    }
	return ans;
  }
/*********************************************************************
   Returns the Bessel function Y_1(x) for any positive real x.
   C.A. Bertulani        May/16/2000
*********************************************************************/
  double bessy1(double x) {
    double bessj1(double x);
    double z;
    double xx,y,ans,ans1,ans2;	/* Accumulate polynomials in double precision.  */
    
	if (x < 8.0) {		/* Rational function approximation.  */
	  y=x*x;
	  ans1=x*(-0.4900604943e13+y*(0.1275274390e13
				      +y*(-0.5153438139e11+y*(0.7349264551e9
							      +y*(-0.4237922726e7+y*0.8511937935e4)))));
	  ans2=0.2499580570e14+y*(0.4244419664e12
				  +y*(0.3733650367e10+y*(0.2245904002e8
							 +y*(0.1020426050e6+y*(0.3549632885e3+y)))));
	  ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x);
	} else {			/* Fitting function. */
	  z=8.0/x;
	  y=z*z;
	  xx=x-2.356194491;
	  ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
				     +y*(0.2457520174e-5+y*(-0.240337019e-6))));
	  ans2=0.04687499995+y*(-0.2004690873e-3
				+y*(0.8449199096e-5+y*(-0.88228987e-6
						       +y*0.105787412e-6)));
	  ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
	}
	return ans;
  }
/*********************************************************************
  Returns the Bessel function Y_n(x) for any positive real x and n>=2.
  C.A. Bertulani        May/16/2000
*********************************************************************/
  double bessy(int n, double x) {
    double bessy0(double x);
    double bessy1(double x);
    int j;
    double by,bym,byp,tox;
    
    if (n < 2) std::cerr << "Index n less than 2 in bessy";
    tox=2.0/x;
    by=bessy1(x);			/* Starting values for recurrence. */
    bym=bessy0(x);
    for (j=1;j<n;j++) {		/* Recurrence. */
      byp=j*tox*by-bym;
      bym=by;
      by=byp;
    }
	return by;
  }
/*********************************************************************
   Returns the spherical Bessel function j_n(x), yn(x), and their
   respective derivatives, for any positive real x and n>=0.
   For details, see Numerical Recipes, Chap. 6.	
  C.A. Bertulani        May/16/2000
*********************************************************************/
  
// #define RTPIO2 1.2533141
  
  void sphbes(int n, double x, double *sj, double *sy, double *sjp, double *syp) {
   //  void bessjy(double x, double xnu, double *rj, double *ry, double *rjp,
// 		double *ryp);
    const double RTPIO2 = 1.2533141; 
    double factor,order,rj,rjp,ry,ryp;
    
    if (n < 0 || x <= 0.0) std::cerr << "bad arguments in sphbes\n";
    order=n+0.5;
    bessjy(x,order,&rj,&ry,&rjp,&ryp);
    factor=RTPIO2/sqrt(x);
    *sj=factor*rj;
    *sy=factor*ry;
    *sjp=factor*rjp-(*sj)/(2.0*x);
    *syp=factor*ryp-(*sy)/(2.0*x);
  }
  // #undef RTPIO2
/*********************************************************************
   Needed to calculate spherical Bessel functions. Returns the Bessel
   fnctions rj=J_nu, ry = Y_nu and their derivatives rjp and ryp, for
   positive x and for xnu (nu) >=0. 
   For details, see Numerical Recipes in C, Chap. 6.		
   C.A. Bertulani        May/16/2000
*********************************************************************/
//#include "butil.h"
// #define EPS 1.0e-16
// #define FPMIN 1.0e-30
// #define MAXIT 10000
// #define XMIN 2.0
// #define PI 3.141592653589793
  
  void bessjy(double x, double xnu, double *rj, double *ry, double *rjp, double *ryp) {
    const double EPS = 1.0e-16;
    const double FPMIN = 1.0e-30;
    const int MAXIT = 10000;
    const double XMIN = 2.0;
    const double PI = 3.141592653589793;
//     void beschb(double x, double *gam1, double *gam2, double *gampl,
// 		double *gammi);
    int i,isign,l,nl;
    double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2,
      fact3,ff,gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,rjl,
      rjl1,rjmu,rjp1,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sum1,
      temp,w,x2,xi,xi2,xmu,xmu2;
    
    if (x <= 0.0 || xnu < 0.0) std::cerr << "bad arguments in bessjy\n";
    nl=(x < XMIN ? (int)(xnu+0.5) : Max(0,(int)(xnu-x+1.5)));
    /* n is the number of downward recurrences of the J's and upward recurrences */
    /* of Y's. xmu lies between -1/2 and 1/2 for x < XMIN, while it is chosen so */
    /* that x is greater than the turning point for x>= XMIN. */
    xmu=xnu-nl;
    xmu2=xmu*xmu;
    xi=1.0/x;
    xi2=2.0*xi;
    w=xi2/PI;			/* The Wronskian.  */
    isign=1;			/* Lentz's method */
    h=xnu*xi;			/* isign keeps track of sign changes in the denominator. */
    if (h < FPMIN) h=FPMIN;
    b=xi2*xnu;
    d=0.0;
	c=h;
	for (i=1;i<=MAXIT;i++) {
	  b += xi2;
	  d=b-d;
	  if (fabs(d) < FPMIN) d=FPMIN;
	  c=b-1.0/c;
	  if (fabs(c) < FPMIN) c=FPMIN;
	  d=1.0/d;
	  del=c*d;
	  h=del*h;
	  if (d < 0.0) isign = -isign;
	  if (fabs(del-1.0) < EPS) break;
	}
	if (i > MAXIT)std::cerr << "x too large in bessjy; try asymptotic expansion\n";
	rjl=isign*FPMIN;	/* Initialize J_nu and J'_nu for downward recurrence. */
	rjpl=h*rjl;
	rjl1=rjl;			/* Store values for later r`escaling. */
	rjp1=rjpl;
	fact=xnu*xi;
	for (l=nl;l>=1;l--) {
	  rjtemp=fact*rjl+rjpl;
	  fact -= xi;
	  rjpl=fact*rjtemp-rjl;
	  rjl=rjtemp;
	}
	if (rjl == 0.0) rjl=EPS;
	f=rjpl/rjl;
	if (x < XMIN) {	/* Now have unnormalized J_nu and J'_nu. */
	  x2=0.5*x;	/* Use series. */
	  pimu=PI*xmu;
	  fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu));
	  d = -log(x2);
	  e=xmu*d;
	  fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e);
	  beschb(xmu,&gam1,&gam2,&gampl,&gammi);
	  ff=2.0/PI*fact*(gam1*cosh(e)+gam2*fact2*d);
	  e=exp(e);
	  p=e/(gampl*PI);
	  q=1.0/(e*PI*gammi);
	  pimu2=0.5*pimu;
	  fact3 = (fabs(pimu2) < EPS ? 1.0 : sin(pimu2)/pimu2);
	  r=PI*pimu2*fact3*fact3;
	  c=1.0;
	  d = -x2*x2;
	  sum=ff+r*q;
	  sum1=p;
	  for (i=1;i<=MAXIT;i++) {
	    ff=(i*ff+p+q)/(i*i-xmu2);
	    c *= (d/i);
	    p /= (i-xmu);
	    q /= (i+xmu);
	    del=c*(ff+r*q);
	    sum += del;
	    del1=c*p-i*del;
	    sum1 += del1;
	    if (fabs(del) < (1.0+fabs(sum))*EPS) break;
	  }
 	  if (i > MAXIT)std::cerr << "bessy series failed to converge\n";
	  rymu = -sum;
	  ry1 = -sum1*xi2;
	  rymup=xmu*xi*rymu-ry1;
	  rjmu=w/(rymup-f*rymu);
	} else {
	  a=0.25-xmu2;
	  p = -0.5*xi;
	  q=1.0;
	  br=2.0*x;
	  bi=2.0;
	  fact=a*xi/(p*p+q*q);
	  cr=br+q*fact;
	  ci=bi+p*fact;
	  den=br*br+bi*bi;
	  dr=br/den;
	  di = -bi/den;
	  dlr=cr*dr-ci*di;
	  dli=cr*di+ci*dr;
	  temp=p*dlr-q*dli;
	  q=p*dli+q*dlr;
	  p=temp;
	  for (i=2;i<=MAXIT;i++) {
	    a += 2*(i-1);
	    bi += 2.0;
	    dr=a*dr+br;
	    di=a*di+bi;
	    if (fabs(dr)+fabs(di) < FPMIN) dr=FPMIN;
	    fact=a/(cr*cr+ci*ci);
	    cr=br+cr*fact;
	    ci=bi-ci*fact;
	    if (fabs(cr)+fabs(ci) < FPMIN) cr=FPMIN;
	    den=dr*dr+di*di;
	    dr /= den;
	    di /= -den;
	    dlr=cr*dr-ci*di;
	    dli=cr*di+ci*dr;
	    temp=p*dlr-q*dli;
	    q=p*dli+q*dlr;
	    p=temp;
	    if (fabs(dlr-1.0)+fabs(dli) < EPS) break;
	  }
 	  if (i > MAXIT)std::cerr << "cf2 failed in bessjy\n";
	  gam=(p-f)/q;
	  rjmu=sqrt(w/((p-f)*gam+q));
	  rjmu=Sign(rjmu,rjl);
	  rymu=rjmu*gam;
	  rymup=rymu*(p+q/gam);
	  ry1=xmu*xi*rymu-rymup;
	}
	fact=rjmu/rjl;
	*rj=rjl1*fact;		/* Scale original J_nu and Y_nu.  */
	*rjp=rjp1*fact;
	for (i=1;i<=nl;i++) {	/* Upward recurrence of Y_nu. */
	  rytemp=(xmu+i)*xi2*ry1-rymu;
	  rymu=ry1;
	  ry1=rytemp;
	}
	*ry=rymu;
	*ryp=xnu*xi*rymu-ry1;
  }
// #undef EPS
// #undef FPMIN
// #undef MAXIT
// #undef XMIN
// #undef PI
/*********************************************************************
   Needed to calculate spherical Bessel function j_n(x).
   For details, see Numerical Recipes in C, Chap. 6.	
   C.A. Bertulani        May/16/2000
*********************************************************************/
// #define NUSE1 7
// #define NUSE2 8
  
  void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi) {
//     double chebev(double a, double b, double c[], int m, double x);
    const int NUSE1 = 7;
    const int NUSE2 = 8;
    double xx;
    static double c1[] = {
      -1.142022680371172e0,6.516511267076e-3,
	3.08709017308e-4,-3.470626964e-6,6.943764e-9,
	3.6780e-11,-1.36e-13};
    static double c2[] = {
      1.843740587300906e0,-0.076852840844786e0,
	1.271927136655e-3,-4.971736704e-6,-3.3126120e-8,
	2.42310e-10,-1.70e-13,-1.0e-15};
    
    xx=8.0*x*x-1.0;
    *gam1=chebev(-1.0,1.0,c1,NUSE1,xx);
    *gam2=chebev(-1.0,1.0,c2,NUSE2,xx);
    *gampl= *gam2-x*(*gam1);
    *gammi= *gam2+x*(*gam1);
  }
// #undef NUSE1
// #undef NUSE2
/*********************************************************************
   Needed to calculate spherical Bessel function j_n(x).
   For details, see Numerical Recipes in C, Chap. 6.	
   C.A. Bertulani        May/16/2000
*********************************************************************/
  double chebev(double a, double b, double c[], int m, double x) {
    double d=0.0,dd=0.0,sv,y,y2;
    int j;
    
     if ((x-a)*(x-b) > 0.0) std::cerr << "x not in range in routine chebev";
    y2=2.0*(y=(2.0*x-a-b)/(b-a));
    for (j=m-1;j>=1;j--) {
      sv=d;
      d=y2*d-dd+c[j];
      dd=sv;
    }
    return y*d-dd+0.5*c[0];
  }
/*********************************************************************/


/********************************************************************
                   Discrete version of Gaussian kernel                           
********************************************************************/
  double scaleSpaceGaussian(double x,double t) {
    int n = int(fabs(x));
    if (n == 0)
      return exp(-t)*bessi0(t);
    else if (n == 1)
      return exp(-t)*bessi1(t);
    else
      return exp(-t)*bessi(n,t);
  }
  
  /** Uniformly distributed numbers in the range [Min..Max] by linear congruences. Seed from timestamp */
  double UniformDistribution(double low, double high) {
    static bool initialized = false;
    static RandomNumber randnum;
    if (!initialized) {
      randnum.init();
      initialized = true;
    }
    double ran = double(randnum.random()*(high-low) + low);
    return ran;

//     assert(high > low);
//     static bool seed_done = false;
//     static int seed;
//     if (!seed_done) {
//       Imlib::timeS t;
//       time(&t);
//       seed = int(t);
//       seed_done = true;
//     }
//     srand(seed);
//     int intrand = rand();
//     double ran = double(intrand / (double)RAND_MAX) * (high - low) + low;
//     seed = intrand;
//     return ran;
  }
  
  /** Normally distributed variables with mean and stdv. Use the law of large numbers (where 12 is "large")
   *  y = Sum(xi), i=1..12, 
   * where E[xi] = 1/2, V[xi] = 1/12.
   * E[y] - 6 = 0
   * V[y] = 1
   * Transform to scaled random variable 
   * z = stdv * y + mean
   **/
  double GaussianDistribution(double mean, double stdv) {
    double rsum=0.0;
    for (int i=0;i<12;i++)
      rsum += UniformDistribution(0.0,1.0);
    return stdv*(rsum - 6) + mean;

  }


  /** An function approximation of the normalized cumulative distribution function Phi(x)
   *          
   *                      x 
   *            1         / 
   * Phi(x) = -----      | exp(-t^2/2) dt = prob( y < x)
   *         sqrt(2Pi)   /
   *                   -00
   *
   *
   * This is related to the error function integral Erf(x)
   *
   *                        x
   *            2          / 
   * Erf(x) = -----       | exp(-t^2) dt
   *         sqrt(Pi)     /
   *                     0
   * 
   * and
   *
   * Phi(x) = 1/2 (1+Erf(x/sqrt(2)))
   *
   *
   *
   * from Bagby (http://mathworld.wolfram.com/NormalDistributionFunction.html)
   *
   * Phi2(x) ~ Phi(x) = 1/2(1 + Sqrt[1-1/30*(7*exp(-x^2/2) + 16*exp(-x^2(2-Sqrt[2])))
   *                                                   + (7+Pi*x^2/4)*exp(-x^2))]
   * The function has a maximal error of 3 10^-5 over the interval {0,Infinity}
   *
   *
   *
   *
   */
  double errorFunctionPhi(double x) {
    return 0.5 + (0.5*sqrt(1.0-1/30.0*(7*exp(-x*x/2) + 16*exp(-x*x*(2.0-sqrt(2.0)) +
                                                      (7+Pi*x*x/4)*exp(-x*x)))));

  }





} // namaspace mathutil