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 real *8 (a-h,o-z)
      real *8 dexpi,dt,dum,dtemp,de(6)
      include "consts.inc"
      include "es.inc"
c
c     new statement 9/4/93
      include "depolt.inc"
c     end of new statement

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 
      return
      end