Deriving Luminosity from an SED

This notebook demonstrates how to use the dust_emissivity module to compute the luminosity of a dusty source.

In [1]: 
from __future__ import print_function
from dust_emissivity.blackbody import modified_blackbody,integrate_sed,blackbody
from astropy import units as u
from astropy import constants
import pylab as pl
In [2]: 
frequencies = np.logspace(2, 4, 1000)*u.GHz
fmin,fmax = frequencies.min(), frequencies.max()
print("Range: ",fmax.to(u.um, u.spectral()), " to ", fmin.to(u.um, u.spectral()))

Range:  29.9792458 um  to  2997.92458 um
In [3]: 
temperature = 25*u.K
In [4]: 
# modified blackbody returns normal bb units
modified_blackbody(fmin, temperature, beta=0).to(u.erg/u.s/u.cm**2/u.Hz)
Out[4]:
$1.30397\times 10^{-16} \; \mathrm{\frac{erg}{Hz\,s\,cm^{2}}}$
In [5]: 
pl.loglog(frequencies, modified_blackbody(frequencies, temperature, beta=0, column=1e22*u.cm**-2), color='b')
pl.loglog(frequencies, modified_blackbody(frequencies, temperature, beta=1, column=1e22*u.cm**-2), color='r')
pl.loglog(frequencies, modified_blackbody(frequencies, temperature, beta=2, column=1e22*u.cm**-2), color='g')
pl.loglog(frequencies, modified_blackbody(frequencies, temperature, beta=0, column=1e21*u.cm**-2), color='b', linestyle='--')
pl.loglog(frequencies, modified_blackbody(frequencies, temperature, beta=1, column=1e21*u.cm**-2), color='r', linestyle='--')
pl.loglog(frequencies, modified_blackbody(frequencies, temperature, beta=2, column=1e21*u.cm**-2), color='g', linestyle='--')
Out[5]:
[<matplotlib.lines.Line2D at 0x10893f750>]
In [6]: 
integrated_sed = integrate_sed(fmin, fmax, function=modified_blackbody, temperature=25*u.K, beta=1.75)
print(integrated_sed)

0.168799551811 erg / (cm2 s)

To compute the luminosity, we need to measure the emitting area. For resolved, dusty objects, this is simply the angular area on the sky times the distance. For example, a 25"x30" arcsecond Gaussian blob at a distance of 3.4 kpc:

In [7]: 
# Area of a gaussian is 2*pi*a*b
angular_area = 2*np.pi*(25*u.arcsec*30*u.arcsec)
distance = 3.4*u.kpc
# Convert the angles to radians and treat them as unitless
physical_area = (angular_area*distance**2).to(u.cm**2, u.dimensionless_angles())
#physical_area = (angular_area_in_arcseconds / arcsec_per_radian**2) * distance_in_cm**2
print("Source area: ",physical_area," and effective physical radius: ", (physical_area**0.5).to(u.pc))

Source area:  1.21912935864e+37 cm2  and effective physical radius:  1.13155156532 pc
In [8]: 
luminosity = integrated_sed * physical_area
print(luminosity.to(u.L_sun))

535.071475137 solLum

Now let's go back the other way and determine how bright the source would be. We are aiming to determine its brightness in MJy/sr.

In [9]: 
flux_density = modified_blackbody((500*u.um).to(u.GHz,u.spectral()), temperature=25*u.K, beta=1.75, column=1e22*u.cm**-2)
print(flux_density,flux_density.to(u.Jy))

3.71537524828e-15 erg / (cm2 Hz s) 371537524.828 Jy
In [10]: 
# Total energy emitted at the source location in a spectral bin per second
e_per_bin = flux_density*physical_area
print(e_per_bin)

4.52952304354e+22 erg / (Hz s)
In [11]: 
# Dilute that by the distance**2
received_e_per_bin = (e_per_bin / distance**2).to(u.erg/u.s/u.Hz/u.cm**2)
print(received_e_per_bin)

4.11522465561e-22 erg / (cm2 Hz s)
In [12]: 
# Then divide out the amount of sky that energy occupies
surface_brightness = received_e_per_bin/angular_area
print(surface_brightness)

8.73277794499e-26 erg / (arcsec2 cm2 Hz s)
In [13]: 
# Convert that to the desired unit
print(surface_brightness.to(u.MJy/u.sr))

371.537524828 MJy / sr
In [14]: 
# Put that all together:
surface_brightness = (flux_density * physical_area / distance**2 / angular_area).to(u.MJy/u.sr)
surface_brightness = ((flux_density / distance**2 * physical_area) / angular_area).to(u.MJy/u.sr)
print(surface_brightness)

371.537524828 MJy / sr
In [15]: 
(flux_density/u.sr).to(u.MJy/u.sr)
Out[15]:
$371.538 \; \mathrm{\frac{MJy}{sr}}$

Vega

Tests with Vega to see that we're sane

In [16]: 
b_vega = blackbody((0.5*u.um).to(u.Hz,u.spectral()), 10000*u.K)
b_vega.to(u.Jy)
Out[16]:
$1.89515\times 10^{+19} \; \mathrm{Jy}$
In [17]: 
r_vega = 2.818*constants.R_sun
d_vega = (25.05*u.lyr).to(u.pc)
angular_area_vega = (r_vega/d_vega)**2 * np.pi * u.sr
In [18]: 
(b_vega * (r_vega**2 * np.pi) / d_vega**2).to(u.Jy)
Out[18]:
$4072.07 \; \mathrm{Jy}$
In [19]: 
(b_vega * angular_area_vega).to(u.Jy, u.dimensionless_angles())
Out[19]:
$4072.07 \; \mathrm{Jy}$
In [20]: 
((b_vega * angular_area_vega).to(u.Jy, u.dimensionless_angles()) / angular_area_vega).to(u.MJy/u.sr)
Out[20]:
$1.89515\times 10^{+13} \; \mathrm{\frac{MJy}{sr}}$

An example specific to The Snake

In [21]: 
T = 24*u.K
angular_area = (8*u.arcsec)**2
distance = 3.6*u.kpc
physical_area = (angular_area.to(u.rad**2)*distance**2).to(u.cm**2, u.dimensionless_angles())
beta = 1.75
In [22]: 
nu = (np.array([70,160,250,350,500])*u.um).to(u.Hz,u.spectral())
column = 1e23*u.cm**-2
In [23]: 
flux_density = modified_blackbody(nu, T, column=column, beta=beta)
In [24]: 
surface_brightness = ((flux_density / distance**2 * physical_area) / angular_area).to(u.MJy/u.sr)
surface_brightness2 = ((flux_density / (u.sr))).to(u.MJy/u.sr)
surface_brightness,surface_brightness2
Out[24]:
(<Quantity [ 12069.25897604, 39742.21507926, 20739.83151157,
              9405.08410352,  3427.14928482] MJy / sr>,
 <Quantity [ 12069.25897604, 39742.21507926, 20739.83151157,
              9405.08410352,  3427.14928482] MJy / sr>)
In [25]: 
pl.loglog(nu, surface_brightness)
pl.xlabel("Frequency")
pl.ylabel("Brightness (MJy/sr)")
Out[25]:
<matplotlib.text.Text at 0x1089e73d0>