Autor: Eduardo Martín Calleja
El objeto de este post es construir un diagrama de Hertsprung-Russel y posicionar en el mismo diferentes clases de estrellas
Importaciones y referencias¶
%matplotlib inline
from __future__ import division
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, pandas
Una interesante actividad educativa acerca de los diagramas de Hertsprung-Russel es The Hertzsprung-Russell Diagram en el SLOAN DIGITAL SKY SURVEY III (SDSS III).
Algo de teoría¶
Los hechos relevantes para la elaboración de un diagrama de Hertsprung-Russel (diagrama H-R en breve) son los siguientes:
Por un lado existe una relación entre la magnitud absoluta de una estrella y su luminosidad, la cual ha sido descrita en este post: El color de los objetos celestes: parte I
También existe una relación entre la clase espectral de la estrella (O B A F G K M), su temperatura superficial y su índice de color (por ej. B-V). Esta relación entre índice de color y temperatura se puede consultar también en el mismo post anterior. En cuanto a la relación entre clases espectrales y temperaturas superficiales, puede consultarse en mi post El color de los objetos celestes: parte II
El objeto de los diagramas H-R es mostrar la relación existente entre ambos grupos de magnitudes. El diagrama H-R "clásico" muestra la clase espectral de las estrellas (o su correspondiente índice de color B-V) en el eje horizontal y las magnitudes absolutas en el eje vertical.
La posición de una estrella en el diagrama H-R nos dice algo acerca de sus carácterísticas. Así por ejemplo si una estrella presenta una clase espectral que indica una temperatura superficial "fria", por ejemplo clases K o M, pero su magnitud absoluta / luminosidad es elevada, esto apunta a que el astro posee una gran superficie, con lo que la elevada superficie radiativa "compensa" la baja intensidad de su radiación debido a su baja temperatura. Es decir, nos hallaríamos ante una estrella "gigante"o "supergigante". Como los índices de color correspondientes a temperaturas frias tienden al rojo, hablaríamos en este caso por ejemplo de gigantes rojas.
El catálogo HIPPARCOS¶
La mayor dificultad para construir un diagrama H-R es determinar la luminosidad o magnitud absoluta de las estrellas, ya que ello supone conocer sus distancias. En este sentido el catálogo de objetos celestes del SDSS que he venido utilizando hasta ahora en artículos anteriores no es el más adecuado ya que su principal propósito es la fotometría de objetos extragalácticos, disponiendo de poca información de estrellas de nuestra galaxia, y desde luego careciendo de información precisa de la distancia que nos separa de ellas.
Por esta razón, el catálogo de estrellas en el que vamos a basarnos es HIPPARCOS, ya que el objetivo priomordial de esta misión fue precisamente realizar mediciones precisas de cerca de 120000 estrellas. Además, la información en dicho catálogo incluye la clase espectral de la gran mayoría de estos objetos, por lo que es la fuente más adecuada para realizar un diagrama H-R.
Una descripción de los campos en el catálogo puede descargarse en formato pdf en la página siguiente: Hipparcos Overview. En particular, serán relevantes para este artículo los campos siguientes:
- H1 identificador en el catálogo (HIP)
- H5 magnitud en el filtro V del sistema Johnson (Vmag)
- H11 paralaje trigonométrico en milisegundos de arco (Plx)
- H37 índice de color (B-V)
- H76 Spectral type (SpType)
Las magnitudes absolutas se calcularán a partir de los valores de V. Para ello se puede partir de la relación entre la magnitud absoluta \(M\), la magnitud aparente \(m\) y la paralaje \(p\) en arcosegundos (ver por ejemplo el artículo Absolute magnitude):
\[ M = m + 5 (1 + \log_{10} p )\]
Cuando la paralaje se expresa en milisegundos de arco (Plx), y tomando como magnitud aparente el valor \(V\) del catálogo, obtenemos:
\[ M = Vmag + 5 (1 + \log_{10} \frac{Plx}{1000}) \]
Que, simplificando queda:
\[ M = Vmag + 5 * log_{10} \frac{Plx}{100} \]
Para descargar el catálogo se debe acceder a la página VizieR Search Page del catálogo I/239 del Centro de Datos Astronómicos de Estrasburgo (CDS). Seleccionar solo los campos indicados anteriormente: HIP, Vmag, Plx, B-V, SpType. En el panel izquierdo, en "Preferences" seleccionar "unlimited" y "; -Separated Values". Pulsar el botón "Submit" y descargar el fichero, que yo he llamado "I_239_selection.tsv" y solo ocupa 3.8 MB. Abriendolo con un editor de textos se comprueba que se trata de un fichero de texto con unas líneas de cabecera explicativas de su contenido, seguidas de 118322 entradas con los cinco campos seleccionados separados por el carácter ;. Al final del fichero hay una línea en blanco que debemos omitir al leer el fichero.
A partir de este momento supondremos pues descargado el fichero "I_239_selection.tsv" en el directorio "../Hipparcos" A continuación se va a leer el fichero creando un dataframe de Pandas:
filename = '../Hipparcos/I_239_selection.tsv'
df = pd.read_table(filename, skiprows=44, sep=';', header=None, index_col=0,
names = ['HIP', 'Vmag', 'Plx', 'B-V', 'SpType'],
skipfooter=1, engine='python')
df.head()
df.tail()
Preparación del dataframe con los datos¶
Vamos ahora a trabajar sobre este dataframe. En primer lugar, se analizará si hay entradas en el catálogo sin datos válidos (por ejemplo, con valores NaN)
df.describe()
El valor de "count" indica que de las 118322 entradas disponibles en el catálogo, faltan valores en cada una de las tres columnas, que corresponderán a valores no existentes (NaN). Vamos a identificar y eliminar esas entradas.
Previamente vamos a cambiar cualquier campo con espacios a valores NaN, ya que estos valores no los considera Pandas como NaN, pero no obstante deben ser eliminados ya que con ellos no se podrá realizar cálculos (el método a continuación está tomado de la siguiente respuesta en Stackoverflow
df_clean = df.applymap(lambda x: np.nan if isinstance(x, basestring) and x.isspace() else x)
df_clean.describe()
Se puede comprobar que los "count" son menores que los anteriores, lo que significa que habia valores, sobre todo en la columna Plx, con espacios, los cuales no fueron considerados como NaN. Ahora están marcados todos como NaN y finalmente, se eliminarán todas las filas que contienen algún valor NaN
df_clean= df_clean.dropna()
df_clean.describe()
df_clean.shape
Se comprueba que el contador de valores no nulos ahora coincide con el número de líneas del dataframe. A partir de ahora se trabajará con el dataframe df_clean.
A continuación se va a calcular, mediante la fórmula anteriormente indicada, la magnitud absoluta en el filtro V, a la que llamaremos 'M_V'. Se añadirá una columna al dataframe, pero previamente se van a cambiar las columnas del dataframe con las que se va a operar, convirtiendolas de strings a valores numéricos.
df_clean['Vmag'] = df_clean['Vmag'].astype(np.float)
df_clean['Plx'] = df_clean['Plx'].astype(np.float)
df_clean['B-V'] = df_clean['B-V'].astype(np.float)
df_clean['M_V'] = df_clean['Vmag'] + 5 * np.log10(df_clean['Plx']/100.)
df_clean.head()
A continuación se van a estudiar los tipos espectrales presentes en el dataframe. Realmente la única información que se necesita de momento de la columna SpType son los dos primeros caracteres. Además solo interesan aquellos casos en los que el primer carácter es alfabético y el segundo numérico, como G8 o F5. Una vez eliminadas las anomalías, se creará una nueva columna con solo los dos primeros caracteres:
# Se eliminan las filas que no cumplen la condición alfa + num
f = lambda s: (len(s) >= 2) and (s[0].isalpha()) and (s[1].isdigit())
i = df_clean['SpType'].apply(f)
df_clean = df_clean[i]
# Se crea una columna nueva con los dos primeros caracteres de 'SpType'
f = lambda s: s[0:2]
df_clean['SpType2'] = df_clean['SpType'].apply(f)
df_clean.shape
df_clean.head()
Veamos cuales son las clases espectrales presentes
f = lambda s: s[0]
clases = df_clean['SpType'].map(f)
clases.value_counts()
Se van a eliminar las líneas con las clases especiales C, N, R y S para quedarnos solo con las de la secuencia OBAFGKM
f = lambda s: s[0] in 'OBAFGKM'
df_clean = df_clean[df_clean['SpType'].map(f)]
Comprobación:
f = lambda s: s[0]
clases = df_clean['SpType'].map(f)
clases.value_counts()
El objetivo que se persigue es realizar un gráfico donde en el eje horizontal figuren las clases espectrales en el sistema de Morgan–Keenan (MKK): OBAFGKM, seguidas de un dígito 0-9. Pero además se pretende que las clases espectrales aparezcan ordenadas precisamente en este orden (B5 antes de A0...). Por ello se van a sustituir las letras por dígitos en secuencia creciente en la columna SpType2 con el fin de poder ordenar por esa columna.
orden = {'O':'0', 'B':'1', 'A':'2', 'F':'3', 'G':'4', 'K':'5', 'M':'6'}
f = lambda s: orden[s[0]]+s[1]
df_clean['SpType2'] = df_clean['SpType2'].apply(f)
df_clean.head()
Representación gráfica¶
fig, ax = plt.subplots(figsize=(8,10))
ax.set_xlim(0, 70)
ax.set_ylim(15, -10)
ax.grid()
ax.set_title(u'Diagrama H-R \n (Catálogo Hipparcos)')
ax.title.set_fontsize(20)
ax.set_xlabel('Clase espectral')
ax.xaxis.label.set_fontsize(20)
ax.set_ylabel('Magnitud absoluta')
ax.yaxis.label.set_fontsize(20)
ax.scatter(df_clean['SpType2'].astype(np.int), df_clean['M_V'],
s=50, edgecolors='none', alpha=0.015, c='k')
ax.set_xticks(range(5,75,10))
ax.set_xticklabels(['O', 'B', 'A', 'F', 'G', 'K', 'M'])
ax.tick_params(axis='both', labelsize=14)
El diagrama anterior no es muy satisfactorio, en el sentido de que, al corresponder el eje horizontal a un conjunto discreto de valores (las diferentes clases espectrales), se forman bandas verticales por acumulación de observaciones en esos valores. Sería preferible utilizar una escala continua de valores decimales. Por este motivo se va a repetir el diagrama utilizando en esta ocasión el índice de color B-V en lugar de las clases espectrales. Además, ello permitirá comprobar que ambas representaciones presentan la misma forma confirmando la relación entre clases espectrales e índices de color.
fig, ax = plt.subplots(figsize=(8,10))
ax.set_xlim(-0.5, 2.5)
ax.set_ylim(15, -10)
ax.grid()
ax.set_title(u'Diagrama H-R \n (Catálogo Hipparcos)')
ax.title.set_fontsize(20)
ax.set_xlabel(u'Índice de color B-V')
ax.xaxis.label.set_fontsize(20)
ax.set_ylabel('Magnitud absoluta')
ax.yaxis.label.set_fontsize(20)
ax.scatter(df_clean['B-V'], df_clean['M_V'],
# s=50, edgecolors='none', alpha=0.015, c='k')
s=1, edgecolors='none', c='k')
ax.tick_params(axis='both', labelsize=14)
Clases de luminosidad¶
Además del tipo espectral OBAFGKM, muchas estrellas en el catálogo Hipparcos contienen, en la codificación del tipo de espectro (SpType), la clase de luminosidad en números romanos: I a VII, que nos indican el tipo de estrella de que se trata. Puede leerse una explicación de este sistema de clasificación en el artículo de wikipedia: Stellar classification.
A continuación se va a analizar la presencia de estas clases de luminosidad en el catálogo:
f = lambda s: 'VII' in s
b = df_clean['SpType'].map(f)
print "Clase VII: enanas blancas, hay %d estrellas" %sum(b)
f = lambda s: ('VI' in s) and ('VII' not in s)
b = df_clean['SpType'].map(f)
print "Clase VI: subenanas, hay %d estrellas" %sum(b)
f = lambda s: ('V' in s) and ('VI' not in s) and ('IV' not in s)
b = df_clean['SpType'].map(f)
print "Clase V: secuencia principal, hay %d estrellas" %sum(b)
f = lambda s: 'IV' in s
b = df_clean['SpType'].map(f)
print "Clase IV: subgigantes, hay %d estrellas" %sum(b)
f = lambda s: 'III' in s
b = df_clean['SpType'].map(f)
print "Clase III: gigantes, hay %d estrellas" %sum(b)
f = lambda s: ('II' in s) and ('III' not in s) and ('VII' not in s)
b = df_clean['SpType'].map(f)
print "Clase II: gigantes brillantes, hay %d estrellas" %sum(b)
f = lambda s: ('I' in s) and ('II' not in s) and ('V' not in s)
b = df_clean['SpType'].map(f)
print "Clase I: supergigantes, hay %d estrellas" %sum(b)
Como se puede ver a continuación, aproximadamente la mitad de las estrellas no tienen codificada una clase de luminosidad en el campo SpType:
f = lambda s: ('I' not in s) and ('V' not in s)
b = df_clean['SpType'].map(f)
print sum(b)
No obstante, existe un importante número de entradas en el catálogo que sí contienen información de la clase estelar, por lo que a continuación se va a mostrar la posición que ocupa cada una de estas clases de luminosidad en el diagrama H-R
def plot_lum_class(b,c, label):
''' b: Series booleana para hacer la selección
c: Color
label: etiqueta para la leyenda
'''
x = df_clean['B-V'][b]
y = df_clean['M_V'][b]
ax.scatter(x, y, c = c, s=6, edgecolors='none', label = label)
fig = plt.figure(figsize=(8,10))
ax = fig.add_subplot(111, axisbg='0.8')
ax.set_xlim(-0.5, 2.5)
ax.set_ylim(15, -15)
ax.grid()
ax.set_title(u'Diagrama H-R \n (Catálogo Hipparcos)')
ax.title.set_fontsize(20)
ax.set_xlabel(u'Índice de color B-V')
ax.xaxis.label.set_fontsize(20)
ax.set_ylabel('Magnitud absoluta')
ax.yaxis.label.set_fontsize(20)
f = lambda s: 'VII' in s
b = df_clean['SpType'].map(f)
plot_lum_class(b,'white', 'VII: enanas blancas')
f = lambda s: ('VI' in s) and ('VII' not in s)
b = df_clean['SpType'].map(f)
plot_lum_class(b,'blue', 'VI: subenanas')
f = lambda s: ('V' in s) and ('VI' not in s) and ('IV' not in s)
b = df_clean['SpType'].map(f)
plot_lum_class(b,'black', 'V: secuencia principal')
f = lambda s: 'IV' in s
b = df_clean['SpType'].map(f)
plot_lum_class(b,'grey', 'IV: subgigantes')
f = lambda s: 'III' in s
b = df_clean['SpType'].map(f)
plot_lum_class(b,'green', 'III: gigantes')
f = lambda s: ('II' in s) and ('III' not in s) and ('VII' not in s)
b = df_clean['SpType'].map(f)
plot_lum_class(b,'orange', 'II: gigantes brillantes')
f = lambda s: ('I' in s) and ('II' not in s) and ('V' not in s)
b = df_clean['SpType'].map(f)
plot_lum_class(b,'yellow', 'I: supergigantes')
ax.tick_params(axis='both', labelsize=14)
legend = ax.legend(scatterpoints=1,markerscale = 6, shadow=True)
frame = legend.get_frame()
frame.set_facecolor('0.90')
Aquí concluye este artículo. No dudeis en enviarme vuestros comentarios. ¡Hasta la próxima!
¡Lo incomprendo todo!
ResponderEliminar¡huff!
ResponderEliminarExcelente trabajo! Muchísimas gracias por compartirlo!
ResponderEliminarTengo interés en hacer algo parecido con 2MASS, pero no me deja bajar el catálogo completo como vos lo hiciste con Hipparcos a través de VizieR. Tenés alguna idea al respecto?
Saludos,
Germán.