function dei ( x1 )

c*********************************************************************72
c
c  an exponential integral routine.
c
c  for x greater than 0, the exponential integral, ei, is defined by
c  ei(x) = integral ( exp ( t ) / t dt ), from t = -infinity to t = x
c  where the integral is to be interpreted as the cauchy principal
c  value.  for x less than 0, ei(x) = -e1(-x), where
c  e1(z) = integral ( exp ( -t ) / t dt ) from t = z to t = infinity.
c
c  modified:
c
c    04 october 2006
c
c  reference:
c
c    kathleen paciorek,
c    algorithm 385:
c    exponential integral ei(x),
c    communications of the acm,
c    volume 13, number 7, july 1970, pages 446-447.
c
      implicit none

      double precision a(6)
      double precision b(6)
      double precision c(8)
      double precision d(8)
      double precision dei
      double precision denm
      double precision e(8)
      double precision f(8)
      double precision frac
      integer i
      integer j
      integer l
      double precision p0(6)
      double precision p1(9)
      double precision p2(9)
      double precision p3(10)
      double precision p4(10)
      double precision px(10)
      double precision q0(6)
      double precision q1(9)
      double precision q2(8)
      double precision q3(9)
      double precision q4(9)
      double precision qx(10)
      double precision r
      double precision sump
      double precision sumq
      double precision t
      double precision w
      double precision x
      double precision x0
      double precision x1
      double precision xmx0
      double precision y

      save a
      save b
      save c
      save d
      save e
      save f
      save p0
      save p1
      save p2
      save p3
      save p4
      save q0
      save q1
      save q2
      save q3
      save q4
      save x0

      data a /
     &  -5.77215664901532863d-01,
     &   7.54164313663016620d-01,
     &   1.29849232927373234d-01,
     &   2.40681355683977413d-02,
     &   1.32084309209609371d-03,
     &   6.57739399753264501d-05 /
      data b /
     &  1.0d+00,
     &  4.25899193811589822d-01,
     &  7.9779471841022822d-02,
     &  8.30208476098771677d-03,
     &  4.86427138393016416d-04,
     &  1.30655195822848878d-05 /
      data c /
     &  8.67745954838443744d-08,
     &  9.99995519301390302d-01,
     &  1.18483105554945844d+01,
     &  4.55930644253389823d+01,
     &  6.99279451291003023d+01,
     &  4.25202034768840779d+01,
     &  8.83671808803843939d+00,
     &  4.01377664940664720d-01 /
      data d /
     &  1.0d+00,
     &  1.28481935379156650d+01,
     &  5.64433569561803199d+01,
     &  1.06645183769913883d+02,
     &  8.97311097125289802d+01,
     &  3.14971849170440750d+01,
     &  3.79559003762122243d+00,
     &  9.08804569188869219d-02 /
      data e /
     &  -9.99999999999973414d-01,
     &  -3.44061995006684895d+01,
     &  -4.27532671201988539d+02,
     &  -2.39601943247490540d+03,
     &  -6.16885210055476351d+03,
     &  -6.57609698748021179d+03,
     &  -2.10607737142633289d+03,
     &  -1.48990849972948169d+01 /
      data f /
     &  1.0d+00,
     &  3.64061995006459804d+01,
     &  4.94345070209903645d+02,
     &  3.19027237489543304d+03,
     &  1.03370753085840977d+04,
     &  1.63241453557783503d+04,
     &  1.11497752871096620d+04,
     &  2.37813899102160221d+03 /
      data p0 /
     &  1.0d+00,
     &  2.23069937666899751d+00,
     &  1.70277059606809295d+00,
     &  5.10499279623219400d-01,
     &  4.89089253789279154d-02,
     &  3.65462224132368429d-04 /
      data p1 /
     &  5.99569946892370010d+09,
     & -2.50389994886351362d+08,
     &  7.05921609590056747d+08,
     & -3.36899564201591901d+06,
     &  8.98683291643758313d+06,
     &  7.37147790184657443d+04,
     &  2.85446881813647015d+04,
     &  4.12626667248911939d+02,
     &  1.10639547241639580d+01 /
      data p2 /
     &  9.98957666516551704d-01,
     &  5.73116705744508018d+00,
     &  4.18102422562856622d+00,
     &  5.88658240753281111d+00,
     & -1.94132967514430702d+01,
     &  7.89472209294457221d+00,
     &  2.32730233839039141d+01,
     & -3.67783113478311458d+01,
     & -2.46940983448361265d+00 /
      data p3 /
     &  9.99993310616056874d-01,
     & -1.84508623239127867d+00,
     &  2.65257581845279982d+01,
     &  2.49548773040205944d+01,
     & -3.32361257934396228d+01,
     & -9.13483569999874255d-01,
     & -2.10574079954804045d+01,
     & -1.00064191398928483d+01,
     & -1.86009212172643758d+01,
     & -1.64772117246346314d+00 /
      data p4 /
     &  1.00000000000000486d+00,
     & -3.00000000320981266d+00,
     & -5.00006640413131002d+00,
     & -7.06810977895029359d+00,
     & -1.52856623636929637d+01,
     & -7.63147701620253631d+00,
     & -2.79798528624305389d+01,
     & -1.81949664929868906d+01,
     & -2.23127670777632410d+02,
     &  1.75338801265465972d+02 /
      data q0 /
     &  1.0d+00,
     &  2.73069937666899751d+00,
     &  2.73478695106925836d+00,
     &  1.21765962960151532d+00,
     &  2.28817933990526412d-01,
     &  1.31114151194977706d-02 /   
      data q1 /
     &  2.55926497607616350d+09,
     & -2.79673351122984591d+09,
     &  8.02827782946956507d+08,
     & -1.44980714393023883d+08,
     &  1.77158308010799884d+07,
     & -1.49575457202559218d+06,
     &  8.53771000180749097d+04,
     & -3.02523682238227410d+03,
     &  5.12578125d+01 /
      data q2 /
     &  1.14625253249016191d+00,
     & -1.99149600231235164d+02,
     &  3.41365212524375539d+02,
     &  5.23165568734558614d+01,
     &  3.17279489254369328d+02,
     & -8.38767084189640707d+00,
     &  9.65405217429280303d+02,
     &  2.63983007318024593d+00 /
      data q3 /
     &  1.00153385204534270d+00,
     & -1.09355619539109124d+01,
     &  1.99100447081774247d+02,
     &  1.19283242396860101d+03,
     &  4.42941317833792840d+01,
     &  2.53881931563070803d+02,
     &  5.99493232566740736d+01,
     &  6.40380040535241555d+01,
     &  9.79240359921729030d+01 /
      data q4 /
     &  1.99999999999048104d+00,
     & -2.99999894040324960d+00,
     & -7.99243595776339741d+00,
     & -1.20187763547154743d+01,
     &  7.04831847180424676d+01,
     &  1.17179220502086455d+02,
     &  1.37790390235747999d+02,
     &  3.97277109100414518d+00,
     &  3.97845977167414721d+04 /
      data x0 / 0.372507410781366634d+00 /

      x = x1
      if ( x .le. 0.0d+00 ) go to 100
      if ( x .ge. 12.0d+00 ) go to 60
      if ( x .ge. 6.0d+00 ) go to 40
c
c  x in (0,6).
c
      t = x + x
      t = t / 3.0d+00 - 2.0d+00
      px(10) = 0.0d+00
      qx(10) = 0.0d+00
      px(9) = p1(9)
      qx(9) = q1(9)
c
c  the rational function is expressed as a ratio of finite sums of
c  shifted chebyshev polynomials, and is evaluated by noting that
c  t*(x) = t(2*x-1) and using the clenshaw-rice algorithm found in
c  reference (4).
c
      do l = 2, 8
        i = 10 - l
        px(i) = t * px(i+1) - px(i+2) + p1(i)
        qx(i) = t * qx(i+1) - qx(i+2) + q1(i)
      end do

      r = ( 0.5d+00 * t * px(2) - px(3) + p1(1) ) 
     &  / ( 0.5d+00 * t * qx(2) - qx(3) + q1(1) )
c
c  ( x - x0 ) = ( x - x1 ) - x2, where x1 = 409576229586. / 2**40 and
c  x2 = -.7671772501993940d-12.
c
      xmx0 = ( x - 409576229586.0d+00 / 1099511627776.0d+00 ) 
     & - 0.7671772501993940d-12
      if ( dabs ( xmx0 ) .lt. 0.037d+00 ) go to 15
      dei = dlog ( x / x0 ) + xmx0 * r
      return
15    y = xmx0 / x0
c
c  a rational approximation to log ( x / x0 ) * log ( 1 + y ), 
c  where y = ( x - x0 ) / x0, and dabs ( y ) is less than 0.1,
c  that is for dabs ( x - x0 ) less than 0.037.
c
      sump = (((( p0(6) 
     &      * y + p0(5) ) 
     &      * y + p0(4) ) 
     &      * y + p0(3) ) 
     &      * y + p0(2) ) 
     &      * y + p0(1)

      sumq = (((( q0(6) 
     &      * y + q0(5) ) 
     &      * y + q0(4) ) 
     &      * y + q0(3) ) 
     &      * y + q0(2) ) 
     &      * y + q0(1)

      dei = ( sump / ( sumq * x0 ) + r ) * xmx0
      return
c
c  x in (6,12).
c
40    denm = p2(9) + x
      frac = q2(8) + x
c
c  the rational function is expressed as a j-fraction.
c
      do j = 2, 8
        i = 9 - j
        denm = p2(i+1) + x + frac
        frac = q2(i) / denm
      end do

      dei = dexp ( x ) * ( ( p2(1) + frac ) / x )
      return

60    if ( x .ge. 24.0d+00 ) go to 80
c
c  x in (12,24).
c
      denm = p3(10) + x
      frac = q3(9) / denm
c
c  the rational function is expressed as a j-fraction.
c
      do j = 2, 9
        i = 9 - j
        denm = p3(i+1) + x + frac
        frac = q3(i) / denm
      end do

      dei = dexp ( x ) * ( ( p3(1) + frac ) / x )
      return
c
c  x greater than 24.
c
80    if ( x .le. 174.673d+00 ) go to 90
c
c  x is greater than 174.673 and dei is set to infinity.
c
      dei = 7.2d+75
      return
90    y = 1.0d+00 / x
      denm = p4(10) + x
      frac = q4(9) / denm
c
c  the rational function is expressed as a j-fraction.
c
      do j = 2, 9
        i = 10 - j
        denm = p4(i+1) + x + frac
        frac = q4(i) / denm
      end do

      dei = dexp ( x ) * ( y + y * y * ( p4(1) + frac ) )
      return

100   if ( x .ne. 0.0d+00 ) go to 101
c
c  x = 0 and dei is set to -infinity.
c
      dei = -7.2d+75
      write ( *, 500 )
500   format ( 
     & 'dei called with a zero argument, result set to -infinity' )
      return
101   y = -x
      w = 1.0d+00 / y
      if ( y .gt. 4.0d+00 ) go to 300
      if ( y .gt. 1.0d+00 ) go to 200
c
c  x in (-1,0).
c
      dei = dlog ( y ) - (((((
     &        a(6)
     &  * y + a(5) )
     &  * y + a(4) )
     &  * y + a(4) )
     &  * y + a(2) )
     &  * y + a(1) ) / (((((
     &        b(6)
     &  * y + b(5) )
     &  * y + b(4) )
     &  * y + b(4) )
     &  * y + b(2) )
     &  * y + b(1) )
      return
c
c  x in (-4,-1).
c
200   dei = -dexp ( -y ) * ((((((((
     &        c(8)
     &  * w + c(7) )
     &  * w + c(6) )
     &  * w + c(5) )
     &  * w + c(4) )
     &  * w + c(3) )
     &  * w + c(2) )
     &  * w + c(1) ) / (((((((
     &        d(8)
     &  * w + d(7) )
     &  * w + d(6) )
     &  * w + d(5) )
     &  * w + d(4) )
     &  * w + d(3) )
     &  * w + d(2) )
     &  * w + d(1) ) )
      return
c
c  x less than -4.
c
300   dei = -dexp ( -y ) * ( w * ( 1.0d+00 + w * (((((((
     &        e(8)
     &  * w + e(7) )
     &  * w + e(6) )
     &  * w + e(5) )
     &  * w + e(4) )
     &  * w + e(3) )
     &  * w + e(2) )
     &  * w + e(1) ) / (((((((
     &        f(8)
     &  * w + f(7) )
     &  * w + f(6) )
     &  * w + f(5) )
     &  * w + f(4) )
     &  * w + f(3) )
     &  * w + f(2) )
     &  * w + f(1) ) ) )

      return
      end