subroutine photo(flam,rho,mue,cv,cvp,ca,cap,n,qphot) implicit double precision (a-h,o-z) dimension a(0:2),b(3) double precision mue xi=(((rho/mue)/1e9)**(1./3.))/flam a(0)=4.886d10 a(1)=7.580d10 a(2)=6.023d10 b(1)=6.290d-3 b(2)=7.483d-3 b(3)=3.061d-4 c=1.5654d0 call qcalc(a,b,c,xi,flam,f) qphot=(rho/mue)*(flam**5.)*f temp=1.875d8*flam+1.654d8*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) bigq=bigq*0.893e9*flam**8*((1+143.8d0*flam**3.555d0)**0.3516d0) qphot=bigq*xi**3*dexp((0.556d0*xi**4.48d0)/(150.d0+xi**3.48d0))*f return end ************************************************************************