subroutine bremsolo(dens,t,gamma,qbremso)
*
* This routine returns the bremsstrahlung neutrino rates computed for 
* Solid Oxygen using the results of 
* Itoh, Kohyama, Matsumoto, & Seki 1984, ApJ, 285, 304 
* (+Errata 1987, ApJ, 322, 584)
* Includes lattice and phonon (lattice vibrations) contributions.
*
      implicit double precision (a-z)
*
* declare variables (find the values in the block data piece of code).
*
      common/soldato/ 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
      common/phondato/ ap(5), bp(3), ip(5), jp(3),
     1          alphap(4), betap(4), cp, dp, kp, lp

      at = 16.
      z = 8.
      twopi = 6.283185308
      pisqr = 9.869604404
      en = 2.
      s2theta = 0.217
      cv = 0.5 + 2.*s2theta
      ca = 0.5
      cvp = 1. - cv
      cap = 1. - ca
      x13 = -1./3.
      x23 = -2./3.
      t8 = 10.**t / 1.e+08
*
* compute the interpolation coefficients v and w
*
      u = twopi * ( dens - 3. ) / 9.
*
* lattice variables If gamma>5000., set to 5000. for lattice results.
*
      if(gamma.gt.5000.)then
        gamma2 = 5000.
       else
        gamma2 = gamma
      endif
      v = alpha(1) + alpha(2) * gamma2**x13 + alpha(3) * gamma2**x23
     1    + alpha(4) / gamma2
      w = beta(1) + beta(2) * gamma2**x13 + beta(3) * gamma2**x23
     1    + beta(4) / gamma2
*
* phonon variables
*
      vp = alphap(1) + alphap(2)*gamma**x13 + alphap(3)*gamma**x23
     1    + alphap(4) / gamma
      wp = betap(1) + betap(2)*gamma**x13 + betap(3)*gamma**x23
     1    + betap(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.
*
*  Lattice contributions
*
      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
*
*  Phonon contributions
*
      apsum = ap(1) + ap(2)*cu + ap(3)*c2u
     1       + ap(4)*c3u + ap(5)*c4u
      bpsum = bp(1)*su + bp(2)*s2u + bp(3)*s3u
     1       + cp*u + dp
      ipsum = ip(1) + ip(2)*cu + ip(3)*c2u
     1       + ip(4)*c3u + ip(5)*c4u
      jpsum = jp(1)*su + jp(2)*s2u + jp(3)*s3u
     1       + kp*u + lp
*
* sum up Lattice contributions
*
      flat171 = asum + bsum
      flat5k  = esum + fsum
      glat171 = isum + jsum
      glat5k  = psum + qsum
      flat    = (1.-v) * flat171 + v * flat5k
      glat    = (1.-w) * glat171 + w * glat5k
*
* sum up Phonon contributions
* We just extrapolate this is gamma > 5000.
*
      fplat171 = vp*(apsum + bpsum)
      gplat171 = wp*(ipsum + jpsum)
*
* sum up for total neutrino loss contributions in solid phase
*
      ftot = flat + fplat171
      gtot = glat + gplat171
*
* compute neutrino loss rate
*
      first = 0.5*(cv*cv + ca*ca + en*cvp*cvp + en*cap*cap)*ftot
      sec = 0.5*(cv*cv - ca*ca + en*cvp*cvp - en*cap*cap)*gtot
      qbremso = 0.5738*(z*z/at)*(t8**6)*(first - sec)

      return
      end 

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