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 ************************************************************************