# Python imports, bringing in all the useful tools we need
# some are abbreviated for convenience
# some we import specific modules that we know we want to use
import astropy
import numpy as np
from scipy import integrate
from astropy import units as u, constants as c
from astropy.modeling.blackbody import blackbody_lambda, blackbody_nu, BlackBody1D
from matplotlib import pyplot as plt
import pickle as pkl
import pandas
# this is a 'magic' function specific to jupyter notebooks, making plots appear inline
%matplotlib inline
# run with ctrl-c right away to make sure these all import without errors
We'll use astropy's units and constants tools along with their blackbody tools to find the luminosity of a blackbody using the Planck function.
The Planck function gives us the intensity at a given wavelength (or frequency). The total intensity is the intensity over all frequencies.
First, let's explore what we have to work with.
u.
c.
blackbody_nu?
Signature: blackbody_nu(inx, temperature) Docstring: Calculate blackbody flux per steradian, :math:`B{\nu}(T)`.
.. note::
Use `numpy.errstate` to suppress Numpy warnings, if desired.
.. warning::
Output values might contain ``nan`` and ``inf``.
Parameters¶
in_x : number, array-like, or ~astropy.units.Quantity
Frequency, wavelength, or wave number.
If not a Quantity, it is assumed to be in Hz.
temperature : number, array-like, or ~astropy.units.Quantity
Blackbody temperature.
If not a Quantity, it is assumed to be in Kelvin.
Returns¶
flux : ~astropy.units.Quantity
Blackbody monochromatic flux in
:math:erg \; cm^{-2} s^{-1} Hz^{-1} sr^{-1}.
Astropy can be used to make a Quantity object, which has both a value and a unit. To do this, you "multiply" the value by the unit object, accessed through the imported units module, like this:
example_frequency = 50.7 * u.GHz
example_frequency.value
example_frequency.unit
This also enables conversions to other units like this:
example_frequency.to(u.MHz)
More complex conversions can be done using equivalencies, like the spectral equivalency for converting between wavelength, frequency, and wavenumber
# This does not work, it cannot naively convert GHz to nanometer
example_frequency.to(u.nm)
example_frequency.to(u.nm, equivalencies=u.spectral())
example_frequency.to(u.k, equivalencies=u.spectral())
(30. * u.k).to(u.GHz,equivalencies=u.spectral())
Let's use the blackbody_nu function to generate the frequency spectrum of a 3 Kelvin blackbody
nu = np.linspace(30.,300.,1000) * u.GHz
bb3K_nu = blackbody_nu(in_x = nu, temperature = 3. * u.Kelvin )
Inspect the first element to see what it is, and its units
bb3K_nu[0]
plt.plot(nu,bb3K_nu)
We can use $LaTeX$ markup to add nice-looking axis labels.
We enclose $LaTeX$ markup text in dollar signs, within a string r'\$ ... \$'.
The r before the open-quote denotes that the string is "raw," and backslashes are treated literally.
plt.plot(nu, bb3K_nu)
plt.xlabel(r'$\nu$, [{0:latex_inline}]'.format(nu.unit))
plt.ylabel(r'$I_{\nu}$, ' + '[{0:latex_inline}]'.format(bb3K_nu.unit))
plt.title('Planck function in frequency')
We can do a simple integration using numpy's trapz function to get the integrated flux over the frequency bins that we specified
np.trapz(x=nu, y=bb3K_nu).to('erg s-1 cm-2 sr-1')
blackbody_lambda could be used to do the same thing in wavelength units
Blackbody1D is a built in model for blackbodies (slightly redundant with what's shown above)
bb = BlackBody1D(temperature = 3. * u.K)
bb.lambda_max
bb_spectrum = bb(nu)
plt.plot(nu, bb_spectrum)
plt.xlabel(r'$\nu$, [{0:latex_inline}]'.format(nu.unit))
plt.ylabel(r'$I_{\nu}$, ' + '[{0:latex_inline}]'.format(bb_spectrum.unit))
Let's compare this with the COBE FIRAS data
# The name of the file we want to open.
# We have to be careful about protected variable names so don't call your file "file", be more original
# This is a relative path to the file, starting from wherever Jupyter notebook was launched.
# If you're unsure about the relative path, you can always try tab-complete,
# or use the full absolute path
fname = 'DataSets/FIRAS/firas_monopole_spec_v1.txt'
# there are many different ways to load the data contained in a textfile into python
# since this file only contains numbers, it is well suited to numpy's built in loadtxt function.
# I need to open the text file and look at it to figure out how many lines are headers and what separates columns.
data = np.loadtxt(fname, skiprows=18)
print(data)
# take all rows of the 0th column of the array 'data', and give it units of 'k' (inverse cm)
frequency = data[:,0] * u.k
# print it out to examine and make sure its the right part of the dataset
print(frequency)
spectrum = data[:,1] * (u.MJy / u.steradian)
uncertainty = data[:,3] * (u.kJy / u.steradian)
uncertainty_MJy = uncertainty.to(u.MJy/u.steradian)
plt.plot(frequency, spectrum)
plt.errorbar(frequency.value, spectrum.value, yerr = 1000*uncertainty_MJy.value, label = r'$Uncertainty \times 1000$, [{0:latex_inline}]'.format(uncertainty_MJy.unit))
plt.xlabel(r'$Wavenumber$, [{0:latex_inline}]'.format(frequency.unit))
plt.ylabel(r'$Intensity$, [{0:latex_inline}]'.format(spectrum.unit))
plt.legend(loc='best')
plt.title('COBE/FIRAS CMB Monopole Spectrum')
plt.savefig('COBE_FIRAS_monopole.png')
To compare this with our model, we'll need to do some unit conversion gymnastics
freq_ghz = frequency.to(u.GHz, equivalencies = u.spectral())
spectrum_mjy = bb3K_nu.to(u.MJy / u.steradian, equivalencies = u.spectral_density(100. * u.GHz))
plt.plot(freq_ghz, spectrum, label = 'COBE/FIRAS Monopole')
plt.plot(nu, spectrum_mjy, label = "3 K Blackbody Model (Planck's Law)")
plt.legend(loc='best')
plt.xlabel('Frequency (GHz)')
plt.ylabel('Intensity, MJy/steradian')
plt.title('Comparison of COBE Data to Blackbody Model')
plt.savefig('FIRAS_vs_model.png')
bb2p7K = blackbody_nu(in_x = nu, temperature = 2.7 * u.Kelvin ).to(u.MJy / u.steradian, equivalencies = u.spectral_density (100. * u.GHz))
plt.plot(nu, bb2p7K)
plt.plot(freq_ghz, spectrum)
