function en(N,X)
C-------------------------------------------------------------------
C CALCULATES THE EXPONENTIAL INTEGRAL FOR ARBITRARY X AND N
C REF. ABRAMOWITZ AND STEGUN, IN ALL CASES EN(1,X) IS CALCULATED AND
C EN(N,X) IS OBTAINED FROM THE RECURRENCE RELATION 5.1.14
C FOR 0<X<1.0  RELATION 5.1.53 IS USED
C FOR 1.0<X< 15. RELATION 5.1.56 IS USED
C FOR X> 15. EN < 1.0E-7 AND IS SET TO ZERO
C-------------------------------------------------------------------
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION A(6),B(4),C(4)

       A(1)=-0.57721566D0
       A(2)=00.99999193D0
       A(3)=-0.24991055D0
       A(4)=-0.05519968D0
       A(5)=-0.00976004D0
       A(6)=0.00107857D0
       B(1)=8.5733287401D0
       B(2)=1.80590169730D1
       B(3)=8.6347608925D0
       B(4)=0.2677737343D0
       C(1)=9.5733223454D0
       C(2)=2.56329561486D1
       C(3)=2.10996530827D1
       C(4)=3.9584969228D0

       EN = 0.0D0
       IF(X .EQ. 0.0D0) THEN
	 IF(N .EQ. 1) RETURN
	 EN = 1.0D0/(N-1)
	 RETURN
       ELSE IF(X .LT. 1.0D-05) THEN
	  IF(N .EQ. 1) THEN
	     E1 = -DLOG(X) +A(1)
	  ELSE
	     EN = 1.0D0/(N-1)
	     RETURN
	  ENDIF
C
       ELSE IF((X.GT.1.0D-05).AND.(X.LE.1.0D0)) THEN
	  E1 = -DLOG(X)
	  DO 10 I = 1, 6
 10          E1 = E1 + A(I)*(X**(I-1))
       ELSE IF ((X.GT.1.0D0).AND.(X.LT.1.5D1)) THEN
	  TOP = X**4
	  BOT = X**5
	  DO 20 I = 1, 4
	     TOP = TOP + B(I)*(X**(4-I))
 20          BOT = BOT + C(I)*(X**(5-I))
	  E1 = DEXP(-X)*(TOP/BOT)
       ELSE IF(X .GE. 1.5D1) THEN
	  RETURN
       ENDIF
C
       en = E1
       IF(N .EQ. 1) RETURN
       en = DEXP(-X) - X*en
       IF(N .EQ. 2) RETURN
       en = 5.0D-1*(DEXP(-X) - X*en)
       RETURN
       END