      function cndshe(dens,t,gamma)
*
* this function returns the thermal conductivity computed for a specific
* dens, temp, and gamma using the results if itoh et al (1984) for helium
* in the xtal regime.
*
      implicit double precision (a-z)
*
*  declare variables (find the values in the block data
*  piece of code).
*
      common/hedats/ a(5), b(3), e(5), f(3), i(5), j(3), p(5), q(3),
     1          alpha(4), beta(4), c, d, g, h, k, l, r, s
      data at, z, twopi, pisqr /4.,2., 6.283185308, 9.869604404/
 
      x13 = -1./3.
      x23 = -2./3.
      t8 = 10.**t / 1.e+08
      d6 = 10.**dens / 1.e+06 
      d6mu = d6 * z / at
      d6mu23 = d6mu**(2./3.)
*
* compute the interpolation coefficients v and w
*
      u = twopi * ( dens - 3. ) / 9.
      v = alpha(1) + alpha(2) * gamma**x13 + alpha(3) * gamma**x23
     1    + alpha(4) / gamma
      w = beta(1) + beta(2) * gamma**x13 + beta(3) * gamma**x23
     1    + beta(4) / gamma
*
* set up the variables (sines and cosines) for the loops
*
      su  = dsin(u)
      su2 = su * su 
      cu2 = 1. - su2
      cu  = dsqrt( cu2 )
      s2u = 2. * su * cu
      c2u = cu2 - su2
      s3u = su * ( 3. - 4. * su2 )
      c3u = cu * ( 4. * cu2 - 3. )
      c4u = -8. * cu2 * su2 + 1.
*
*  here's the loops, unwrapped
*
      asum = a(1) + a(2) * cu    + a(3) * c2u
     1            + a(4) * c3u + a(5) * c4u
      bsum = b(1) * su + b(2) * s2u + b(3) * s3u
     1       + c*u + d
      esum = e(1) + e(2) * cu    + e(3) * c2u
     1            + e(4) * c3u + e(5) * c4u
      fsum = f(1) * su + f(2) * s2u + f(3) * s3u
     1       + g*u + h
      isum = i(1) + i(2) * cu    + i(3) * c2u
     1            + i(4) * c3u + i(5) * c4u
      jsum = j(1) * su + j(2) * s2u + j(3) * s3u
     1       + k*u + l
      psum = p(1) + p(2) * cu    + p(3) * c2u
     1            + p(4) * c3u + p(5) * c4u
      qsum = q(1) * su + q(2) * s2u + q(3) * s3u
     1       + r*u + s
 
      isig171 = asum + bsum
      isig5k  = esum + fsum
      ikap171 = isum + jsum
      ikap5k  = psum + qsum
 
      isig    = (1.-v) * isig171 + v * isig5k
      ikap    = (1.-w) * ikap171 + w * ikap5k
 
      smgamma = .07832 * (z/at) * dsqrt(d6) / t8
      sg2     = smgamma * smgamma
 
      g0      = 13.00 / dsqrt( 1. + 0.0174 * sg2 ) 
      g2   = (sg2 / pisqr) * (1. + 0.0118 * sg2)**(-3./2.)
 
      fkap    = isig * g0 + ikap * g2
      cstd6   = 1.018 * d6mu23
      nukap   = 9.554e+16 * t8 * dsqrt( 1. + 1./cstd6 ) * fkap
 
      opacity = 4.146d+15 * d6mu / dsqrt( 1. + cstd6 )
     1          * t8 * 1.e+18 / nukap
 
      cndshe = dlog10( opacity )
 
      return
      end 

************************************************************************

