The Blackbody Radiation and the Cosmic Microwave Background (CMB) in Python: part II¶
Author: Eduardo Martín Calleja
In this post, continuation of this one, I will examine, using code written in Python, the graphical representation of the thermal radiation from a blackbody at a temperature of 2.725 K. Then I will include in this graph the results of the measurements of the spectrum of the cosmic microwave background (CMB) radiation by the FIRAS device of the Cosmic Background Explorer satellite COBE, to check the exquisite accuracy of the adjustment of the actual data to the spectrum of a black body at the temperature indicated.
One of the main difficulties when it comes to plot the radiation curves is, IMHO, the diversity of physical units in which these can be expressed. This is why in this post I will examine in some detail, and with the help of the Python package "quantities", the most common alternatives in terms of physical units for the presentation of the results.
In this post, continuation of this one, I will examine, using code written in Python, the graphical representation of the thermal radiation from a blackbody at a temperature of 2.725 K. Then I will include in this graph the results of the measurements of the spectrum of the cosmic microwave background (CMB) radiation by the FIRAS device of the Cosmic Background Explorer satellite COBE, to check the exquisite accuracy of the adjustment of the actual data to the spectrum of a black body at the temperature indicated.
One of the main difficulties when it comes to plot the radiation curves is, IMHO, the diversity of physical units in which these can be expressed. This is why in this post I will examine in some detail, and with the help of the Python package "quantities", the most common alternatives in terms of physical units for the presentation of the results.
Imports and references¶
In [1]:
%matplotlib inline
from __future__ import division
import quantities as pq
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# This IPython magic generates a table with version information
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information numpy, matplotlib, quantities, pandas
Out[1]:
In this entry we will use the same function defined in the previous post to calculate the spectral radiance according to the Plank's law.
In [2]:
def B(wl,T):
'''wl is an array of wavelengths with units of length
T is a temperature in Kelvin
the result is an array of s.r. with units W/(m**2 * nm * sr)
'''
I = 2 * pq.constants.h * (pq.c)**2 / wl**5 * \
1 / (np.exp((pq.constants.h*pq.c \
/ (wl*pq.constants.k*T)).simplified)-1)
return I.rescale(pq.watt/(pq.m**2 * pq.nm *pq.sr))
Blackbody radiation at T = 2.725K¶
Case 1. Units: \(\lambda\) in mm and B in \(W m^{-2} mm^{-1} sr^{-1}\)¶
As already said, the cosmic microwave background is a perfect example of radiation from a blackbody at a temperature of 2.725 K. Let's graph the spectral radiance curve that corresponds to this temperature. In the first place, and to get an idea of the range of magnitudes in which we move, we will compute the wavelength in which it reaches its maximum, for what we will use Wien's law:
In [3]:
TCMB = 2.725 *pq.K
lmax = pq.constants.b / TCMB
print lmax.simplified
I.e., the photons of the CMB have a wavelength of approximately one millimeter. With this information we can assign a range of values to the wavelengths for the calculation of the curve of the spectral radiance, and set the right units for B:
In [4]:
wl = np.arange(0.1,10,0.1) * pq.mm
I = B(wl,TCMB).rescale(pq.watt/(pq.m**2 * pq.mm * pq.sr))
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wl, I*1e7)
ax.set_title('Blackbody spectrum at T = 2.725 K \
\n (wavelength in mm)')
ax.title.set_fontsize(20)
ax.set_xlabel('Wavelength in mm')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Spectral Radiance ($10^{-7} W m^{-2} mm^{-1} sr^{-1}$)')
ax.yaxis.label.set_fontsize(15)
ax.grid()
Case 2. Units: \(\nu\) in GHz and B in \(W m^{-2} GHz^{-1} sr^{-1}\)¶
Quite often you will find graphics with B expressed as a function of the frequency \(\nu\) measured in \(Hz\) rather than depending on wavelength \(\lambda\) as we have been doing so far. In that case, Plank's law is written as:
\[B_\nu(T) = \frac{2 h \nu^3}{c^2} \frac{1}{e^{\frac{h \nu}{k_B T}}-1}\]
Where \(\nu\) es the radiation frequency, measured in \(Hz\)
\[B_\nu(T) = \frac{2 h \nu^3}{c^2} \frac{1}{e^{\frac{h \nu}{k_B T}}-1}\]
Where \(\nu\) es the radiation frequency, measured in \(Hz\)
Now we write the function to calculate B in accordance with the formula above:
In [5]:
def B_f(wf,T):
'''wf is an array of frequencies with units in GHz
T is a temperature in Kelvin.
It returns an array of spectral radiances
with units W/(m**2 x GHz x sr)
'''
wf = wf.rescale(pq.hertz)
I = 2 * pq.constants.h * wf**3 / pq.c**2 * 1 \
/ (np.exp((pq.constants.h*wf \
/ (pq.constants.k*T)).simplified)-1)
return I.rescale(pq.watt/(pq.m**2 * pq.gigahertz *pq.sr))
Let's test the function:
In [6]:
B_f(10**13*pq.hertz, 7000*pq.K)
Out[6]:
And we plot the graph:
In [7]:
wf = np.arange(0.1,1000,1)* pq.gigahertz
I = B_f(wf,TCMB)
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wf, I*10**9) # Scale the units
ax.set_title('Blackbody spectrum at T = 2.725K \
\n (Frequencies en GHz)')
ax.title.set_fontsize(20)
ax.set_xlabel(r'Frequence $ \nu \; (GHz)$')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Spectral Radiance $(10^{-9} W m^{-2} GHz^{-1} sr^{-1})$')
ax.yaxis.label.set_fontsize(15)
ax.grid()
Case 3. Units: \(\nu\) in \(cm^{-1}\) and B in \(W / (m^{2} \: cm^{-1} sr)\)¶
The fact is that the preferred unit of frequency when it comes to represent the CMB spectrum is the "wave number k", which gives the number of wavelengths per distance unit. It is normally expressed in "cycles / cm", or just \(cm ^ {-1} \) as the number of cycles is a magnitude without dimension. We will then see now how to make conversions from GHz to cycles/ cm
In [8]:
#Conversion of 160 GHz to cm**(-1)
frec = 160 * pq.gigahertz
wl = pq.c / frec
print 'wavelength in cm = ', wl.rescale(pq.cm)
k = 1/( wl.rescale(pq.cm))
print 'wavenumber in cycles per cm = ', k
In [9]:
# The conversion can be made easier
# just by dividing by the speed of light.
fec_in_wave_number = (frec/pq.c).rescale(1/pq.cm)
print fec_in_wave_number
As a curiosity, the frequency unit 1/cm is called, in the cgs system of units, a "kayser".
The formula of Plank's law changes when the frequency is entered in cycles per unit length (k). Its form in this case is the following:
\[B_k(T) = 2 h c^2 k^3 \frac{1}{e^{\frac{h c k}{k_B T}}-1}\]
Note: we should not confuse the wavenumber unit (cycles per length unit) with the angular wavenumber, which is usually also represented with the \(k\) symbol but whose value is 2 $ $ times the spacial wavenumber.
We will then define a function to calculate the spectral radiance according to the above formula.
Note: we should not confuse the wavenumber unit (cycles per length unit) with the angular wavenumber, which is usually also represented with the \(k\) symbol but whose value is 2 $ $ times the spacial wavenumber.
We will then define a function to calculate the spectral radiance according to the above formula.
In [10]:
def B_k(wk,T):
'''wk is an array of frequencies with units in cm^(-1)
T is a temperature in Kelvin.
It returns an array of spectral radiances
with units W/(m**2 x cm**(-1) x sr)
'''
wk = wk.rescale(1/pq.m)
I = 2 * pq.constants.h * pq.c**2 * wk**3 * 1 \
/ (np.exp((pq.constants.h*pq.c*wk \
/ (pq.constants.k*T)).simplified)-1)
return I.rescale(pq.watt/(pq.m**2 * pq.cm**(-1)*pq.sr))
A first test:
In [11]:
f = 5.35 * pq.cm**(-1)
B_k(f, TCMB)
Out[11]:
And now, let's plot the graph in the new units:
In [12]:
wk = np.arange(0.1,20,0.1) / pq.cm
I = B_k(wk,TCMB)
I = I*10**7 # We express the units in multiples of 10^(-7)
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wk, I)
ax.set_title('Blackbody spectrum at T = 2.725K \
\n (Frequencies in $cm^{-1}$)')
ax.title.set_fontsize(20)
ax.set_xlabel('Frequency (cycles / cm )')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Spectral Radiance $(10^{-7}W m^{-2} (cm^{-1})^{-1}sr^{-1})$')
ax.yaxis.label.set_fontsize(15)
ax.grid()
Case 4. Units: \(\nu\) in \(cm^{-1}\) and B in \(MJy / sr\)¶
There is yet another unit of measurement for frequencies, the Jansky, which is the one used with the data that we can download from the COBE, and is defined as:
\[ 1 Jy = 10^{-26} \frac{W}{m^2 Hz} \]
A multiple is the Megajansky = 10**6 janskys
Its use is particularly widespread in radio astronomy. To be able to work easily in Python with these units, we will define the Jansky unit through the package quantities check out this post:
\[ 1 Jy = 10^{-26} \frac{W}{m^2 Hz} \]
A multiple is the Megajansky = 10**6 janskys
Its use is particularly widespread in radio astronomy. To be able to work easily in Python with these units, we will define the Jansky unit through the package quantities check out this post:
In [13]:
Jy = pq.UnitQuantity('jansky', 1e-26*pq.watt/(pq.m**2*pq.hertz), symbol = 'Jy')
MJy = pq.UnitQuantity('megajansky', 10**6 * Jy, symbol = 'MJy')
An example follows of conversion of a spectral radiance to this units:
In [14]:
I = 3.34 * 10**(-18) * pq.W * pq.m**(-2) * pq.hertz**(-1) * pq.sr**(-1)
I.rescale(MJy)
Out[14]:
In the next plot we'll set the x-coordinate units to cycles/cm, and the y-coordinate to spectral radiance in \(MJy*sr^{-1}\)
In [15]:
wf = np.arange(0.1,1000,1)* pq.gigahertz
I = B_f(wf,TCMB)
I = I.rescale(MJy*pq.sr**(-1)) # conversion to MJy * sr**(-1)
wk = (wf/pq.c).rescale(1/pq.cm) # convert freq. to 1/cm
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wk, I)
ax.set_title('Blackbody spectrum at T = 2.725K \
\n spectral radiance in MJy / sr')
ax.title.set_fontsize(20)
ax.set_xlabel('Frequency (cycles/cm)')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Spectral radiance ($MJy \, sr^{-1}$)')
ax.yaxis.label.set_fontsize(15)
ax.grid()
CMB spectrum with data from COBE¶
Once explored the various types of units in which the graph of the spectrum of the cosmic microwave background can be represented, we will use the latter format, since it is precisely in these units that the COBE satellite data are provided.
It only remains for us to download the spectrum of the CMB data obtained with the COBE FIRAS instrument, which are available publicly in FIRAS CMB Monopole Spectrum. In order to read and manipulate the file with ease, we will use the Pandas Python package.
The first step will be downloading the file from the indicated address. It is the file "firas_monopole_spec_v1.txt". If you open it with a text editor, you will see that the columns are delimited by spaces, and that the first 18 rows are comments and must be skipped, although you need to read them since they describe the units used and the meaning of each column of data.
In this case I've downloaded the file in the directory "../datos". Now I will create a Pandas dataframe with all the 5 columns of data. As a separator between columns I use a regular expression which means "any number of spaces"
It only remains for us to download the spectrum of the CMB data obtained with the COBE FIRAS instrument, which are available publicly in FIRAS CMB Monopole Spectrum. In order to read and manipulate the file with ease, we will use the Pandas Python package.
The first step will be downloading the file from the indicated address. It is the file "firas_monopole_spec_v1.txt". If you open it with a text editor, you will see that the columns are delimited by spaces, and that the first 18 rows are comments and must be skipped, although you need to read them since they describe the units used and the meaning of each column of data.
In this case I've downloaded the file in the directory "../datos". Now I will create a Pandas dataframe with all the 5 columns of data. As a separator between columns I use a regular expression which means "any number of spaces"
In [16]:
df = pd.read_table('../datos/firas_monopole_spec_v1.txt',
skiprows=18, sep='\s+', header=None,
names =['freq', 'I', 'residual', 'uncert', 'poles'])
Let's show the ten first lines of the table:
In [17]:
df[0:10]
Out[17]:
Finally we will use the first two columns of the dataframe to superimpose these points of the meassured CMB spectrum to the blackbody curve at a temperature of 2.725 K. In the resulting graph, you will see that the measured CMB data fits perfectly to the theoretical spectrum, so we can conclude that the cosmic microwave background radiation follows a perfect blackbody emission pattern.
In [18]:
wf = np.arange(0.1,1000,1)* pq.gigahertz
I = B_f(wf,TCMB)
I = I.rescale(MJy*pq.sr**(-1)) # conversion to MJy * sr**(-1)
wk = (wf/pq.c).rescale(1/pq.cm) # convert freq. to 1/cm
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wk, I)
ax.set_title('Blackbody spectrum at T = 2.725K \
\n with data from the COBE satellite')
ax.title.set_fontsize(20)
ax.set_xlabel('Frequency (cycles/cm)')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Spectral Radiance ($MJy \, sr^{-1}$)')
ax.yaxis.label.set_fontsize(15)
ax.set_xlim(0,25)
ax.set_ylim(0,400)
ax.scatter(df['freq'], df['I'],c='red', s= 50)
ax.grid()
No hay comentarios:
Publicar un comentario