c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/dexpk.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine dexpk(a,b)
c
c***********************************************************************
c
cccc
c     subroutine dexpk
c
c     purpose-
c
c        calculate exponential integrals of orders 1 to 6 and the
c        elements of k1 matrix  (dave's eqn. 4-11)
c
c     method-
c         recursion relation for the exponential integrals  is  used
c            see handbook of math functions-abramowitz and stegun (p229)
c
c     calling sequence-  call dexpk(a,b)
c
c        variables     type   i/o  description
c        ---------     ----   ---  -----------
c        a             r*8     i   abs(a-b) is the argument of the
c        b             r*8     i   exponential  integrals
c
c     result of the subroutine is returned in common block  'es'
c        e(6)          r*8     o   exponential  integrals
c        eek(6)        r*8     o   elements of the k-matrix
c
c     external references -   dexpi
c
c     modifications-  by  k.klenk  8/19/77
c
c        functions k1(1) and k1(2) are also calculated and are returned
c        in eek(4) and eek(5) respectively (dave's eqn. 4.7)
c
c     last modified 9/4/93....zia ahmad
c     purpose: to account for molecular depolarization
c
c     last modified 03/14/95...dave flittner
c     purpose: set pressure scale height used in gravity correction
c     to rayleigh scattering od.  Create new variable pscaleforg and
c     pass in common block consts.
cccc
c***********************************************************************
c
      implicit none
      real*8 a,b
      include 'es.cmn'
      include 'consts.cmn'
      include 'depolt.cmn'
      include 'des.cmn'
c local
      integer i,j
      real *8 dexpi,dt,dum,de(6)
      external dexpi
c

c
c*****calculate exponential integrals
c
      dt = dabs(a-b)
      de(1) = - dexpi(-dt)
      dum = dexp(-dt)
      do 70 j = 2, 6
   70 de(j) = odan(j)*(dum-dt*de(j-1))
      do 80 i=1,6
   80 e(i) = de(i)
c
c*****calculate the k matrix and functions k1(1) and k(2)
c     new statements added here 9/4/93
      if(ipol .eq. 0)then
         eek(1) = 0.375*(e(1) + e(5))
         eek(2) = c38sq2 * (de(3) - de(5))
         eek(3) = 0.75 * (de(1)-2.d+00 * de(3) + de(5))
         eek(4)=0.375*(e(1)+e(3)-2.0*e(5))
         eek(5)=0.1875*(e(1)+2.0*e(3)+e(5))
      else
         eek(1)=q12s*((sdp**2+q**2)*e(1)+2.0d0*q*e(3)+e(5))
         eek(2)=q12s*delp*(q*(e(1)-e(3))+e(3)-e(5))
         eek(3)=q12s*delp**2*(e(1)-2.0d0*e(3)+e(5))
         eek(4)=0.375d0*q2*(e(1)+e(3)-2.0d0*e(5))
         eek(5)=0.1875d0*q2*(e(1)+2.0d0*e(3)+e(5))
      endif 
c
c now for the derivative of the exponential integral with respect
c to the arguement
c
      if(dt.gt.0.0)then
        dedx(1)=-(dum/dt)
      else
        dedx(1)=0.0
      endif
      do j=2,6
         dedx(j)=-de(j-1)
      enddo
c
      if(ipol .eq. 0)then
         deekdx(1) = 0.375*(dedx(1) + dedx(5))
         deekdx(2) = c38sq2 * (dedx(3) - dedx(5))
         deekdx(3) = 0.75 * (dedx(1)-2.d+00 * dedx(3) + dedx(5))
         deekdx(4)=0.375*(dedx(1)+dedx(3)-2.0*dedx(5))
         deekdx(5)=0.1875*(dedx(1)+2.0*e(3)+dedx(5))
      else
         deekdx(1)=q12s*((sdp**2+q**2)*dedx(1)+2.0d0*q*dedx(3)+dedx(5))
         deekdx(2)=q12s*delp*(q*(dedx(1)-dedx(3))+dedx(3)-dedx(5))
         deekdx(3)=q12s*delp**2*(dedx(1)-2.0d0*dedx(3)+dedx(5))
         deekdx(4)=0.375d0*q2*(dedx(1)+dedx(3)-2.0d0*dedx(5))
         deekdx(5)=0.1875d0*q2*(dedx(1)+2.0d0*dedx(3)+dedx(5))
      endif 
      return
      end