       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

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

