### Program Code for Calculation of Exponential Integrals

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

```