Re: [ROOT] exponential integrals

From: Damir Buskulic (buskulic@lapp.in2p3.fr)
Date: Thu Feb 20 2003 - 21:22:35 MET


Since Numerical Recipes is not free software, in the sense that you have
to pay a minimal fee, I don't think that it is possible to include this
code in ROOT. Am I wrong ?

Cheers

Damir

=====================================================================
| Damir Buskulic                  | Universite de Savoie/LAPP       |
|                                 | Chemin de Bellevue, B.P. 110    |
| Tel : +33 (0)450091600          | F-74941 Annecy-le-Vieux Cedex   |
| e-mail: buskulic@lapp.in2p3.fr  | FRANCE                          |
=====================================================================
mailto:buskulic@lapp.in2p3.fr

On Fri, 21 Feb 2003, Allister Levi Sanchez wrote:

> Hi Rooters,
> 
> I was working on some integrals lately and I stumbled on "exponential 
> integrals".  I was surprised that they were not implemented in TMath. 
>  Anyway, I decided to look them up from the online Numerical Recipes in 
> C.  I simply copied the codes for E_{n}(x) and Ei(x), where x>0.  By the 
> way,
> 
> E_{n}(x) = \int_{1}^{infty} { exp(-xt)/t^{n} } dt      and
> Ei(x) = -\int_{-x}^{infty} { exp(-t)/t } dt
> 
> However, since I needed to input negative values for Ei(x), I used the 
> fact that
> E_{1}(x) = -Ei(-x)
> and placed it into the code.
> 
> Anyway, here's the code, just in case someone else will need them like I 
> did.
> 
> //************** EXPONENTIAL INTEGRAL Ei ******
> // define: ei(x) = -\int_{-x}^{\infty}{exp(-t)/t}dt,  for x>0
> // power series: ei(x) = eulerconst + ln(x) + x/(1*1!) + x^2/(2*2!) + ...
> double ei(double x)
> { // taken from Numerical Recipes in C
>   const double euler = 0.57721566; // Euler's constant, gamma
>   const int maxit = 100;           // max. no. of iterations allowed
>   const double fpmin = 1.0e-40;    // close to smallest floating-point 
> number
>   const double eps = 1.0e-30;       // relative error, or absolute error 
> near
>                                    // the zero of Ei at x=0.3725
>    //  I actually changed fpmin and eps into smaller values than in NR
> 
>   int k;
>   double fact, prev, sum, term;
> 
>   // special case
>   if(x < 0) return -expint(1,-x);
> 
>   if(x == 0.0) { cout << "Bad argument for ei(x)" << endl; return -1; }
>   if(x < fpmin) return log(x)+euler;
>   if(x <= -log(eps)) {
>     sum = 0;
>     fact = 1;
>     for(k=1; k<=maxit; k++) {
>       fact *= x/k;
>       term = fact/k;
>       sum += term;
>       if(term < eps*sum) break;
>     }
>     if(k>maxit) { cout << "Series failed in ei(x)" << endl; return -1; }
>     return sum+log(x)+euler;
>   } else {
>     sum = 0;
>     term = 1;
>     for(k=1; k<=maxit; k++) {
>       prev = term;
>       term *= k/x;
>       if(term<eps) break;
>       if(term<prev) sum+=term;
>       else {
>     sum -= prev;
>     break;
>       }
>     }
>     return exp(x)*(1.0+sum)/x;
>   }
> }
> //*********************************************
> 
> //************** EXPONENTIAL INTEGRALS En *****
> // define: E_n(x) = \int_1^infty{exp(-xt)/t^n}dt, x>0, n=0,1,...
> double expint(int n, double x) {
>   // based on Numerical Recipes in C
>   const double euler = 0.57721566; // Euler's constant, gamma
>   const int maxit = 100;           // max. no. of iterations allowed
>   const double fpmin = 1.0e-30;    // close to smallest floating-point 
> number
>   const double eps = 6.0e-8;       // relative error, or absolute error near
>                                    // the zero of Ei at x=0.3725
> 
>   int i, ii, nm1;
>   double a,b,c,d,del,fact,h,psi,ans;
> 
>   nm1=n-1;
>   if(n<0 || x<0 || (x==0 && (n==0 || n==1))) {
>     cout << "Bad argument for expint(n,x)" << endl; return -1;
>   }
>   else {
>     if(n==0) ans=exp(-x)/x;
>     else {
>       if(x==0) ans=1.0/nm1;
>       else {
>     if(x>1) {
>       b=x+n;
>       c=1.0/fpmin;
>       d=1.0/b;
>       h=d;
>       for(i=1; i<maxit; i++) {
>         a = -i*(nm1+i);
>         b += 2.0;
>         d=1.0/(a*d+b);
>         c=b+a/c;
>         del=c*d;
>         h *= del;
>         if(fabs(del-1.0)<eps) {
>           ans=h*exp(-x);
>           return ans;
>         }
>       }
>       cout << "***continued fraction failed in expint(n,x)!!!" << endl;
>       return -1;
>     } else {
>       ans = (nm1!=0 ? 1.0/nm1 : -log(x)-euler);
>       fact=1;
>       for(i=1; i<=maxit; i++) {
>         fact *= -x/i;
>         if(i!=nm1) del = -fact/(i-nm1);
>         else {
>           psi = -euler;
>           for(ii=1; ii<=nm1; ii++) psi += 1.0/ii;
>           del = fact*(-log(x)+psi);
>         }
>         ans += del;
>         if(fabs(del)<fabs(ans)*eps) return ans;
>       }
>       cout << "***series failed in expint!!!" << endl;
>       return -1;
>     }
>       }
>     }
>   }
> 
>   return ans;
> }
> //*********************************************
> 
> 
> 
> 



This archive was generated by hypermail 2b29 : Thu Jan 01 2004 - 17:50:09 MET