def Prad(T,mu):
return (old_div(radiation_constant,<NUM_LIT>)) * T**<NUM_LIT:4><EOL>
P = a/3 * T**4

Parameters
----------
a : float
    Radiation constant [erg cm^-3 K^-4].
T : float
    Temperature [K].