/********************************************************************************
*  This PL/1 subroutine calculates the Exponential Integrals E1-E4, where
*        En(x) = 1∫∞dt e-x.t/tn ,
*  based on the series approximation and the recursive relationships given in 
*  Abramowitz and Stegun, Handbook of Mathematical Functions):
*
*  Input parameter is:
*    X: argument of the Exponential Integral
*
*  Output parameters are:
*   E1,E2,E3,E4: exponential integrals
*
*
*  This program has worked well and correct for me in the past (using IBM PL/1 compilers),  
*  but I can't guarantee this for all compilers, especially if the code is translated
*  into a different language.
*
*  If there are any problems, please let me know and I mention the issue here.
*  (for contact details see  http://www.plasmaphysics.org.uk/feedback.htm ).
*
*                                               Thomas Smid, May 2007
*****************************************************************************/

 EXPIN: PROCEDURE(X,E1,E2,E3,E4) ;					
 DEFAULT RANGE(A:Z) FLOAT DECIMAL VALUE (DEC FLOAT(16)) ;
 DCL A(14) INIT(-0.57721566,0.99999193,-0.24991055,0.05519968,
 -0.00976004,0.00107857,8.5733287401,18.0590169730,8.6347608925,
 0.2677737343,9.5733223454,25.6329561486,21.0996530827,3.9584969228)
 REAL DEC FLOAT(16) ;
 ON UNDERFLOW ; 

 IF X<=0 THEN DO ;
 E2=1; E3=0.5; E4=1/3 ; 
 GO TO ENDE ;
 END ;

 EXPX=EXP(-X) ; 

 IF X<=1 THEN E1=A(1)+A(2)*X+A(3)*X**2+A(4)*X**3+A(5)*X**4
 +A(6)*X**5 -LOG(X) ;
 ELSE E1=EXPX*(X**4+A(7)*X**3+A(8)*X**2+A(9)*X+A(10))/
 (X**5+A(11)*X**4+A(12)*X**3+A(13)*X**2+A(14)*X) ;

 E2=EXPX-X*E1 ; 
 E3=(EXPX-X*E2)/2 ;
 E4=(EXPX-X*E3)/3 ;

 ENDE: END EXPIN ;