C
C PURPOSE: COMPUTES THE N'TH EXPONENTIAL INTEGRAL OF X
C
C INPUTS:
C     X  INDEPENDENT VARIABLE (-100. .LE. X .LE. +100.)
C     N  ORDER OF DESIRED EXPONENTIAL INTEGRAL (1 .LE. N .LE. 8)
C
C OUTPUTS:
C     EXPINT,  THE DESIRED RESULT
C COMMON:
C
C COMMENTS: OCTOBER 6, 1999, P. JUDGE
C 
      FUNCTION EXPINT(N,X,EX)
C
C  OUTPUT - EXPINT,  THE DESIRED RESULT
C           EX,  EXPF(-X)
C  NOTE   RETURNS WITH E1(0)=0, (NOT INFINITY).
C  BASED ON THE SHARE ROUTINE NUEXPI, WRITTEN BY J. W. COOLEY,
C  COURANT INSTITUTE OF MATHEMATICAL SCIENCES, NEW YORK UNIVERSITY
C  OBTAINED FROM RUDOLF LOESER
C-----GENERAL COMPILATION OF 1 AUGUST 1967.
C
C  MODIFIED FOR STANDARD F77 BY P. JUDGE 19 MAR 2001 (OBSOLETE FEATURES REMOVED)
C
      INCLUDE 'PREC'
      DIMENSION TAB(20),XINT(7)
      DATA XINT/1.,2.,3.,4.,5.,6.,7./
      DATA TAB /.2707662555,.2131473101,.1746297218,.1477309984,
     1.1280843565,.1131470205,.1014028126,.0919145454,.0840790292,
     1.0774922515,.0718735405,.0670215610,.0627878642,.0590604044,
     1.0557529077,.0527977953,.0501413386,.0477402600,.0455592945,
     1.0435694088/
      DATA XSAVE /0./
C
      U=X
      IF(U .EQ. 0.)THEN 
         EX=1.
         IF((N-1) .LE. 0) THEN 
            EXPINT=0.
         ELSE 
            EXPINT=1./XINT(N-1)
         ENDIF
         RETURN
      ENDIF
      IF((U-XSAVE) .NE. 0.) THEN 
         XSAVE=U
         XM=-U
         EMX=EXP(XM)
C     
C  SELECT METHOD FOR COMPUTING EI(XM)
C
         IF((XM-24.5) .GE. 0.) GOTO 400
 501     IF((XM-5.) .GE. 0.) GOTO 300
 502     IF((XM+1.) .LT. 0.) THEN 
            GOTO 100
         ELSE 
            GOTO 200
         ENDIF
      ENDIF
 503  EISAVE=-ARG
      EXSAVE=EMX
C
C  NOW RECURSE TO HIGHER ORDERS
C
      IF(N .GT. 1) THEN 
 505     DO 506 I=2,N
            EISAVE=(U*EISAVE-EXSAVE)/(-XINT(I-1))
 506     CONTINUE
      ENDIF
 507  EXPINT=EISAVE
      EX=EXSAVE
      RETURN
C     
C     EI(XM) FOR XM .LT. -1.0
C     HASTINGS POLYNOMIAL APPROXIMATION
C     
  100 ARG=((((((U+8.573328740 )*U+18.05901697  )*U+8.634760893 )*U
     *+.2677737343)/XM)*EMX)/((((U+9.573322345 )*U+25.63295615  )*U
     *+21.09965308  )*U+3.958496923 )
      GOTO 503
C     EI(XM) FOR -1. .LE. XM .LT. 5.0
C     POWER SERIES EXPANSION ABOUT ZERO
  200 ARG=LOG(ABS(XM))
      ARG=((((((((((((((((.41159050E-14*XM+.71745406E-13)*XM+.76404637E-
     *12)*XM+.11395905E-10)*XM+.17540077E-9)*XM+.23002666E-8)*XM+.275360
     *18E-7)*XM+.30588626E-6)*XM+.31003842E-5)*XM+.28346991E-4)*XM+.2314
     *8057E-3)*XM+.0016666574)*XM+.010416668)*XM+.055555572)*XM+.25)*XM+
     *.99999999)*XM+.57721566)+ARG
      GOTO 503
C     
C     EI(XM) FOR 5.0 .LE. XM .LT. 24.5
C     TABLE LOOK-UP AND INTERPOLATION
C     
 300  I=XM+.5
      XZERO=I
      DELTA=XZERO-XM
      ARG=TAB(I-4)
      IF(DELTA .NE. 0.) THEN
         Y=ARG
         DELTAX=DELTA/XZERO
         POWER=1./DELTAX
         DO 304 I=1,7
            POWER=POWER*DELTAX
            Y=((Y-POWER/XZERO)*DELTA)/XINT(I)
            ARG=ARG+Y
            IF((ABS(Y/ARG)-1.E-8) .LT. 0) GOTO 305
 304     CONTINUE
      ENDIF
 305  ARG=EMX*ARG
      GOTO 503
C     EI(XM) FOR 24.5 .LE. XM
C     TRUNCATED CONTINUED FRACTION
 400  ARG=((((XM-15.)*XM+58.)*XM-50.)*EMX)/((((XM-16.)*XM+72.)*XM-96.)
     *     *XM+24.)
      GOTO 503
      END
C
C***************************************************************
C