Autor: Eduardo Martín Calleja
En este post, continuación del anterior examinaremos con ayuda de código escrito en Python la representación gráfica del espectro de radiación térmica de un cuerpo negro a la temperatura de 2.725K. A continuación incluiremos en este gráfico los resultados de las mediciones del espectro de la radiación del fondo cósmico de microondas (CMB) por el dispositivo FIRAS del satélite Cosmic Background Explorer COBE, para comprobar la exquisita precisión del ajuste de los datos reales al espectro de un cuerpo negro a la temperatura indicada.
En particular trataré de aclarar al máximo mediante ejemplos un aspecto que personalmente me ha resultado algo complejo. Se trata de los diferentes sistemas de unidades en los que se puede presentar el gráfico del espectro de la radiación térmica del CMB
Importaciones y referencias¶
%matplotlib inline
from __future__ import division
import quantities as pq
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# Generar un cuadro con versiones de las librerías utilizadas en este notebook
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information numpy, matplotlib, quantities, pandas
En esta entrada vamos a utilizar la función que definimos en el anterior post para calcular la radiancia espectral según la fórmula de Plank
def B(wl,T):
'''wl es un array de longitudes de onda con unidades de longitud
T es una temperatura expresada en Kelvin
el resultado es un array de r.e. con unidades 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))
Radiación del cuerpo negro para T = 2.725K¶
Caso 1. Unidades: \(\lambda\) en mm y B en \(W m^{-2} mm^{-1} sr^{-1}\)¶
Ya se ha avanzado que el fondo cósmico de micrrondas constituye un ejemplo perfecto de radiación de cuerpo negro correspondiente a una temperatura de 2.725K. Vamos a representar gráficamente la curva de la radiancia espectral correspondiente a esta temperatura para hacernos cargo de las magnitudes. En primer lugar, y para hacernos idea del rango de magnitudes en que nos movemos, calculemos la longitud de onda en la que alcanza su máximo, para lo que utilizaremos la ley de Wien:
TCMB = 2.725 *pq.K
lmax = pq.constants.b / TCMB
print lmax.simplified
Es decir, los fotones del CMB poseeran una longitud de onda de aproximadamente un milímetro. Con esta información podemos asignar un rango de valores a las longitudes de onda para el cálculo de la curva de la radiancia espectral, y ajustar las unidades de B:
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('Espectro del cuerpo negro para T = 2.725 K \
\n (Longitudes de onda en mm)')
ax.title.set_fontsize(20)
ax.set_xlabel('Longitud de onda en mm')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Radiancia espectral ($10^{-7} W m^{-2} mm^{-1} sr^{-1}$)')
ax.yaxis.label.set_fontsize(15)
ax.grid()
Caso 2. Unidades: \(\nu\) en GHz y B en \(W m^{-2} GHz^{-1} sr^{-1}\)¶
A menudo nos encontraremos con gráficos en que B viene expresado en función de la frecuencia \(\nu\) medida en \(Hz\) en lugar de obtenerla en función de la longitud de onda \(\lambda\) como hemos estado haciendo hasta ahora. La fórmula de Plank se escribe en este caso:
\[B_\nu(T) = \frac{2 h \nu^3}{c^2} \frac{1}{e^{\frac{h \nu}{k_B T}}-1}\]
Donde \(\nu\) es la frecuencia de la radiación, medida en \(Hz\)
Escribamos la función para calcular B de acuerdo con la fórmula anterior:
def B_f(wf,T):
'''wf es un array de frecuencias con unidades en GHz
T es una temperatura expresada en Kelvin
el resultado es un array de radiancias espectrales
con unidades 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))
Hagamos una prueba de la función:
B_f(10**13*pq.hertz, 7000*pq.K)
Y construyamos la gráfica:
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) # Escalamos las unidades
ax.set_title('Espectro del cuerpo negro para T = 2.725K \
\n (Frecuencias en GHz)')
ax.title.set_fontsize(20)
ax.set_xlabel('Frecuencia ($GHz)$')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Radiancia espectral $(10^{-9} W m^{-2} GHz^{-1} sr^{-1})$')
ax.yaxis.label.set_fontsize(15)
ax.grid()
Caso 3. Unidades: \(\nu\) en \(cm^{-1}\) y B en \(W / (m^{2} cm^{-1} sr)\)¶
El hecho es que la unidad preferida de frecuencias cuando se trata de representar los valores de las frecuencias del CMB es el "número de onda k", que expresa el número de longitudes de onda en una unidad de distancia. Se expresa normalmente en unidades de "ciclos / cm", o bien \(cm^{-1}\) al ser el número de ciclos una magnitud sin dimensión. Veremos a continuación como se pueden hacer las conversiones de GHz a ciclos / cm
# Conversión de 160 GHz en cm**(-1)
frec = 160 * pq.gigahertz
wl = pq.c / frec
print 'longitud de onda en cm = ', wl.rescale(pq.cm)
k = 1/( wl.rescale(pq.cm))
print u'número de onda en ciclos por cm = ', k
# La conversión se puede hacer de forma más abreviada
# simplemente dividiendo por la velocidad de la luz:
fec_en_wave_number = (frec/pq.c).rescale(1/pq.cm)
print fec_en_wave_number
Como curiosidad, la unidad de frecuencia 1/cm es denominada, en el sistema cgs de unidades, un "kayser".
La ley de Plank se escribe de forma diferente de la anterior cuando la frecuencia se introduce en ciclos por unidad de longitud (\(k\)). Su forma en este caso es la siguiente:
\[B_k(T) = 2 h c^2 k^3 \frac{1}{e^{\frac{h c k}{k_B T}}-1}\]
Aviso: no se debe confundir la unidad número de onda (ciclos por unidad de longitud) con la unidad de número de onda angular, que se suele representar también con el símbolo \(k\) pero cuyo valor es \(2 \pi\) veces la unidad anterior.
Definiremos a continuación una función para calcular la radiancia espectral según la fórmula anterior.
def B_k(wk,T):
'''wk es un array de frecuencias con unidades en cm^-1
T es una temperatura expresada en Kelvin
el resultado es un array de radiancias espectrales
con unidades 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))
Una primera comprobación:
f = 5.35 * pq.cm**(-1)
B_k(f, TCMB)
Y a continuación obtendremos el correspondiente gráfico en las nuevas unidades
wk = np.arange(0.1,20,0.1) / pq.cm
I = B_k(wk,TCMB)
I = I*10**7 # Expresamos las unidades en múltiplos de 10^(-7)
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wk, I)
ax.set_title('Espectro del cuerpo negro para T = 2.725K \
\n (Frecuencias en $cm^{-1}$)')
ax.title.set_fontsize(20)
ax.set_xlabel('Frecuencia (ciclos / cm )')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Radiancia espectral $(10^{-7}W m^{-2} (cm^{-1})^{-1}sr^{-1})$')
ax.yaxis.label.set_fontsize(15)
ax.grid()
Caso 4. Unidades: \(\nu\) en \(cm^{-1}\) y B en \(MJy / sr\)¶
Pues bien, Existe aún otra unidad de medida de frecuencia, el Jansky, que es en la que están codificados los datos que podemos descargar del COBE, y que está definido como: \[ 1 Jy = 10^{-26} \frac{W}{m^2 Hz} \]
Un múltiplo suyo será el Megajansky = 10**6 janskys
Su uso está especialmente extendido en radioastronomía. Para poder trabajar con facilidad en Python con estas unidades definiremos dichas unidades mediante el paquete quantities
Jy = pq.UnitQuantity('jansky', 1e-26*pq.watt/(pq.m**2*pq.hertz), symbol = 'Jy')
MJy = pq.UnitQuantity('megajansky', 10**6 * Jy, symbol = 'MJy')
Veamos un ejemplo de conversión de una radianza espectral a estas unidades:
I = 3.34 * 10**(-18) * pq.W * pq.m**(-2) * pq.hertz**(-1) * pq.sr**(-1)
I.rescale(MJy)
Vamos a formar ahora el gráfico poniendo en abcisas frecuencias en ciclos/cm, y en ordenadas el flujo radiante en \(MJy*sr^{-1}\)
wf = np.arange(0.1,1000,1)* pq.gigahertz
I = B_f(wf,TCMB)
I = I.rescale(MJy*pq.sr**(-1)) # convertimos a MJy * sr**(-1)
wk = (wf/pq.c).rescale(1/pq.cm) # convertimos frecuencias a 1/cm
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wk, I)
ax.set_title('Espectro del cuerpo negro para T = 2.725K \
\n radiancia espectral en MJy / sr')
ax.title.set_fontsize(20)
ax.set_xlabel('Frecuencia (ciclos/cm)')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Radiancia espectral ($MJy \, sr^{-1}$)')
ax.yaxis.label.set_fontsize(15)
ax.grid()
Espectro del CMB según datos del COBE¶
Una vez explorados los distintos tipos de unidades en que puede presentarse el gráfico del espectro de radiación del fondo cósmico de microondas, utilizaremos el último formato, ya que es precisamente en estas unidades en que se pueden obtener los datos del satélite COBE.
Solo nos resta descargar los datos del espectro del CMB obtenidos con el instrumento FIRAS del COBE, los cuales están disponibles públicamente en FIRAS CMB Monopole Spectrum. Para leer el fichero y poderlo manipular con facilidad, vamos a utilizar la librería Pandas de Python (a la cual dedicaremos un post en este blog más adelante)
El primer paso va a ser descargar el fichero que se encuentra en la dirección indicada. Se trata del fichero "firas_monopole_spec_v1.txt" en el que, abriendolo con un editor de texto, se observa que las columnas están delimitadas por espacios, y en el cual las 18 primeras filas son comentarios y deben saltarse, si bien debemos leerlas ya que en ellas vienen descritas las unidades utilizadas y el significado de cada columna de datos.
Una vez descargado el fichero en algún subdirectorio de nuestro disco, creamos ahora un dataframe de Pandas con las 5 columnas de datos. El separador utilizado es una expresión regular que significa utilizar como separador entre columnas cualquier número de espacios
#En este caso he descargado el fichero en el directorio ../datos
df = pd.read_table('../datos/firas_monopole_spec_v1.txt',
skiprows=18, sep='\s+', header=None,
names =['freq', 'I', 'residual', 'uncert', 'poles'])
Mostraremos unas cuantas líneas de la tabla obtenida:
df[0:10]
Finalmente vamos a utilizar las dos primeras columnas del dataframe para superponer estos puntos del espectro del CMB medidos por la misión COBE a la anterior curva del espectro de un cuerpo negro a la temperatura de 2.725K. En el gráfico resultante se podrá apreciar que los datos del espectro del CMB se ajustan perfectamente al espectro teorico, pudiendose concluir que la radiación del fondo de microondas cósmico sigue un patron perfecto de cuerpo negro.
wf = np.arange(0.1,1000,1)* pq.gigahertz
I = B_f(wf,TCMB)
I = I.rescale(MJy*pq.sr**(-1)) # convertimos a MJy * sr**(-1)
wk = (wf/pq.c).rescale(1/pq.cm) # convertimos frecuencias a 1/cm
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(wk, I)
ax.set_title(u'Espectro del cuerpo negro para T = 2.725K \
\n con datos del espectro del CMB del satélite COBE')
ax.title.set_fontsize(20)
ax.set_xlabel('Frecuencia (ciclos/cm)')
ax.xaxis.label.set_fontsize(15)
ax.set_ylabel('Radiancia espectral ($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()
