subroutine photon(flam,t2,rho,mue,cv,cvp,ca,cap,qphot)
*
* Photoneutrinos
* a(i) terms are different in the newest paper. Pain to crank out too.
* Coefficients valid for 10^7 to 10^8K, then 10^8 to 10^9K, and
* finally above 10^9K. 
*
       implicit double precision (a-h,o-z)

       dimension a(0:2),b(3),c(3,7),d(3,5)
       double precision mue
       xi=(((rho/mue)/1.d9)**(1./3.))/flam
*
* b(i) terms same as before. Load 'em up.
*
       b(1)=6.290d-3
       b(2)=7.483d-3
       b(3)=3.061d-4
*
* Calculate a(i) numerical factors.
* and tabulated coefficients for sines and cosines.
* tabulated stuff is in form c(i,j), starting with 
* i=1, j=1 instead of 0,0
*
      pi = 3.1415926535898
      half = 1.d0 / 2.d0
      facpi = 5.d0 * pi / 3.d0 
*
* 10^7 to 10^8 K
*
      if(t2.ge.7.0 .and. t2.lt.8.0)then
        tau = t2 - 7.0
        cee = 0.5654 + tau
        c(1,1) = 1.008d+11 * half
        c(1,2) = 0.0
        c(1,3) = 0.0
        c(1,4) = 0.0
        c(1,5) = 0.0
        c(1,6) = 0.0
        c(1,7) = 0.0
        c(2,1) =  8.156d+10 * half
        c(2,2) =  9.728d+8
        c(2,3) = -3.806d+9
        c(2,4) = -4.384d+9
        c(2,5) = -5.774d+9
        c(2,6) = -5.249d+9
        c(2,7) = -5.153d+9 * half
        c(3,1) =  1.067d+11 * half
        c(3,2) = -9.782d+9
        c(3,3) = -7.193d+9
        c(3,4) = -6.936d+9
        c(3,5) = -6.893d+9
        c(3,6) = -7.041d+9
        c(3,7) = -7.193d+9 * half
        d(1,1) = 0.0
        d(1,2) = 0.0
        d(1,3) = 0.0
        d(1,4) = 0.0
        d(1,5) = 0.0
        d(2,1) = -1.879d+10
        d(2,2) = -9.667d+9
        d(2,3) = -5.602d+9
        d(2,4) = -3.370d+9
        d(2,5) = -1.825d+9
        d(3,1) = -2.919d+10
        d(3,2) = -1.185d+10
        d(3,3) = -7.270d+9
        d(3,4) = -4.222d+9
        d(3,5) = -1.560d+9
*
* 10^8 to 10^9 K
*
      elseif(t2.ge.8.0 .and. t2.lt.9.0)then
        tau = t2 - 8.0
        cee = 1.5654
        c(1,1) =  9.889d+10 * half
        c(1,2) = -4.524d+8
        c(1,3) = -6.088d+6
        c(1,4) =  4.269d+7
        c(1,5) =  5.172d+7
        c(1,6) =  4.910d+7
        c(1,7) =  4.388d+7 * half
        c(2,1) =  1.813d+11 * half
        c(2,2) = -7.556d+9
        c(2,3) = -3.304d+9
        c(2,4) = -1.031d+9
        c(2,5) = -1.764d+9
        c(2,6) = -1.851d+9
        c(2,7) = -1.928d+9 * half
        c(3,1) =  9.750d+10 * half
        c(3,2) =  3.484d+9
        c(3,3) =  5.199d+9
        c(3,4) = -1.695d+9
        c(3,5) = -2.865d+9
        c(3,6) = -3.395d+9
        c(3,7) = -3.418d+9 * half
        d(1,1) = -1.135d+8
        d(1,2) =  1.256d+8
        d(1,3) =  5.149d+7
        d(1,4) =  3.436d+7
        d(1,5) =  1.005d+7
        d(2,1) =  1.652d+9
        d(2,2) = -3.119d+9
        d(2,3) = -1.839d+9
        d(2,4) = -1.458d+9
        d(2,5) = -8.956d+8
        d(3,1) = -1.549d+10
        d(3,2) = -9.338d+9
        d(3,3) = -5.899d+9
        d(3,4) = -3.035d+9
        d(3,5) = -1.598d+9
*
* Temperature is greater than 10^9 K
*
      elseif(t2.ge.9.0)then
        tau = t2 - 9.0
        cee = 1.5654
        c(1,1) =  9.581d+10 * half
        c(1,2) =  4.107d+8
        c(1,3) =  2.305d+8
        c(1,4) =  2.236d+8
        c(1,5) =  1.580d+8
        c(1,6) =  2.165d+8
        c(1,7) =  1.721d+8 * half
        c(2,1) =  1.459d+12 * half
        c(2,2) =  1.314d+11
        c(2,3) = -1.169d+11
        c(2,4) = -1.765d+11
        c(2,5) = -1.867d+11
        c(2,6) = -1.983d+11
        c(2,7) = -1.896d+11 * half
        c(3,1) =  2.424d+11 * half
        c(3,2) = -3.669d+9
        c(3,3) = -8.691d+9
        c(3,4) = -7.967d+9
        c(3,5) = -7.932d+9
        c(3,6) = -7.987d+9
        c(3,7) = -8.333d+9 * half
        d(1,1) =  4.724d+8
        d(1,2) =  2.976d+8
        d(1,3) =  2.242d+8
        d(1,4) =  7.937d+7
        d(1,5) =  4.859d+7
        d(2,1) = -7.094d+11
        d(2,2) = -3.697d+11
        d(2,3) = -2.189d+11
        d(2,4) = -1.273d+11
        d(2,5) = -5.705d+10
        d(3,1) = -2.254d+10
        d(3,2) = -1.551d+10
        d(3,3) = -7.793d+9
        d(3,4) = -4.489d+9
        d(3,5) = -2.185d+9
      endif
*
* set up the variables (sines and cosines) for the loops
*
      su1 = dsin(facpi*1.d0*tau)
      su2 = dsin(facpi*2.d0*tau)
      su3 = dsin(facpi*3.d0*tau)
      su4 = dsin(facpi*4.d0*tau)
      su5 = dsin(facpi*5.d0*tau)
      cu1 = dcos(facpi*1.d0*tau)
      cu2 = dcos(facpi*2.d0*tau)
      cu3 = dcos(facpi*3.d0*tau)
      cu4 = dcos(facpi*4.d0*tau)
      cu5 = dcos(facpi*5.d0*tau)
      cu  = dcos(facpi*10.d0*tau)
*
*  Sum up the c's and d's to make a's
*
      csum1 = c(1,1) + c(1,2) * cu1 + c(1,3) *cu2 + 
     1        c(1,4) * cu3 + c(1,5) * cu4 + c(1,6) * cu5
      dsum1 = d(1,1) * su1 + d(1,2) * su2 + d(1,3) *su3 + 
     1        d(1,4) * su4 + d(1,5) * su5 
      a(0)  = csum1 + dsum1 + c(1,7) * cu
      csum2 = c(2,1) + c(2,2) * cu1 + c(2,3) *cu2 + 
     1        c(2,4) * cu3 + c(2,5) * cu4 + c(2,6) * cu5
      dsum2 = d(2,1) * su1 + d(2,2) * su2 + d(2,3) *su3 + 
     1        d(2,4) * su4 + d(2,5) * su5 
      a(1)  = csum2 + dsum2 + c(2,7) * cu
      csum3 = c(3,1) + c(3,2) * cu1 + c(3,3) *cu2 + 
     1        c(3,4) * cu3 + c(3,5) * cu4 + c(3,6) * cu5
      dsum3 = d(3,1) * su1 + d(3,2) * su2 + d(3,3) *su3 + 
     1        d(3,4) * su4 + d(3,5) * su5 
      a(2)  = csum3 + dsum3 + c(3,7) * cu
*
* There, now we have the a(i) coefficients! Lot's o'work.
*
       call qcalc(a,b,cee,xi,flam,f)
       n=2
       temp=1.875d8*flam+1.653d8*flam**2+8.499d8*flam**3-1.604d8*flam**4
       smq=0.666d0*((1.+2.045d0*flam)**-2.066d0)/(1.d0+rho/mue/temp)
       temp=(cv**2-ca**2)+n*(cvp**2-cap**2)
       temp=temp/((cv**2+ca**2)+n*(cvp**2+cap**2))
       bigq=0.5*((cv**2+ca**2)+n*(cvp**2+cap**2))*(1.d0-temp*smq)
       qphot=bigq*f*(rho/mue)*(flam**5)

       return
       end

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