subroutine cia_hhe(t,ah,ahe,ff,opac)
c     ====================================
c
c     CIA H-He opacity
c     data from Gustafsson M., Frommhold, L. 2001, ApJ 546, 1168
c
      IMPLICIT REAL*8(A-H,O-Z)
      parameter (nlines=43)
      dimension freq(nlines),temp(11),alpha(nlines,11)
      parameter (amagat=2.6867774d+19,fac=1./amagat**2)
      data temp / 1000.,  1500.,  2250., 3000.,  4000.,  5000.,
     *            6000.,  7000., 8000.,  9000., 10000./
      data ntemp /11/
      data ifirst /0/
      PARAMETER (CAS=2.997925D10)
c     input frequency in Hz but needed wave numbers in cm^-1
      f=ff/cas
c     read in CIA tables if this is the first call
      if (ifirst.eq.0) then
         write(*,'(a)') 'Reading in H-He CIA opacity tables...'
         open(10,file="./data/CIA_HHe.dat",status='old')
         do i=1,3
            read (10,*)
         enddo
         do i=1,nlines
            read (10,*) freq(i),(alpha(i,j),j=1,ntemp)
         enddo
         close(10)

c     take logarithm of tables prior to doing linear interpolations

         do i=1,nlines
            do j=1,ntemp
               alpha(i,j)=log(alpha(i,j))
            enddo
         enddo

         ifirst=1
      endif

c     locate position in temperature array
      call locate(temp,ntemp,t,j,ntemp)

      if (j.eq.0) then
         write(*,*)
         write(*,'(a,f6.0,a)')
     *   'Warning: requested temperature is below',temp(1),' K'
         write(*,'(a)') 'CIA H-He opacity set to 0'
         write(*,*)
         opac=0.
         return
      endif

c     locate position in frequency array
      call locate(freq,nlines,f,i,nlines)

c     linearly interpolate in frequency and temperature

      if (j.eq.ntemp) then
c     hold values constant if off high temperature end of table
         y1=alpha(i,j)
         y2=alpha(i+1,j)
         tt=(f-freq(i))/(freq(i+1)-freq(i))
         alp=(1.-tt)*y1 + tt*y2
      else if (i.eq.0 .or. i.eq.nlines) then
c     set values to a very small number if off frequency table
         alp=-50.
      else
c     interpolate linearly within table
         y1=alpha(i,j)
         y2=alpha(i+1,j)
         y3=alpha(i+1,j+1)
         y4=alpha(i,j+1)

         tt=(f-freq(i))/(freq(i+1)-freq(i))
         uu=(t-temp(j))/(temp(j+1)-temp(j))

         alp=(1.-tt)*(1.-uu)*y1 + tt*(1.-uu)*y2 + tt*uu*y3 +
     *       (1.-tt)*uu*y4
      endif

      alp=exp(alp)

c     final opacity

      opac=fac*ah*ahe*alp
c
      return
      end