El diagrama de Hertzprung-Russell(HR) es de las herramientas más poderosas para estudiar la evolución de estrellas. Ya que permite clasificarlas en 5 grandes grupos según su ubicación en la gráfica (supergigantes rojas, gigantes azules y rojas, las de secuencia principal y enanans blancas). Para esto se registra datos de luminosidad, radio y temperatura y se ubican en un diagrama de temperatura vs luminosidad, y el tamaño de cada punto está en función de su radio. Dependiendo de la estructura interna de una estrella, su masa inicial y la forma en la que irradia energía, una estrella se desplaza dentro del diagrama HR, es por esto que entender la evolución en la posición de una estrella en este gráfico, permite determinar su estructura interna.
Primero importamos las librerías que vamos a utilizar para realizar el diagrama. Matplotlib y pyplot se van a utilizar para graficar, numpy para utilizar algunas de sus funciones que actúan sobre vectores y pandas para leer los archivos de datos.
import numpy as np
import matplotlib as plt
import pandas as pd
import matplotlib.pyplot
import matplotlib.animation
from IPython.display import HTML
Primero leemos los 4 archivos de datos los cuales proporcionan información sobre enanas blancas, estrellas de la secuencia principal, gigantes y supergigantes. Para esto utilizamos la herramientas de pandas pd.read_cvs, teniendo en cuenta que dos de los archivos están separados por espacios simples y los otros dos por comas. (Para decirle a pandas que interprete los espacios como separaciones utilizamos delimiter=' ').
enana = pd.read_csv("~/ejercicios-clase-03-datos/data/dwarfs.csv")
gigantes = pd.read_csv("~/ejercicios-clase-03-datos/data/giants.txt", delimiter=' ')
supergigantes = pd.read_csv("~/ejercicios-clase-03-datos/data/supergiants.txt", delimiter=' ')
ms = pd.read_csv("~/ejercicios-clase-03-datos/data/ms.csv")
Ahora utilizamos data.head() para entender los datos que hemos leído. Nos damos cuenta que tenemos 3 características de cada estrella. Luminosidad, temperatura y radio.
enana.head()
| lum | temp | radius | |
|---|---|---|---|
| 0 | 0.000109 | 5050.644696 | 7.096930 |
| 1 | 0.000128 | 5967.543450 | 4.583996 |
| 2 | 0.000230 | 6674.161524 | 4.151078 |
| 3 | 0.000269 | 7216.762974 | 3.491754 |
| 4 | 0.000472 | 7795.184395 | 3.472736 |
gigantes.head()
| lum | temp | radius | |
|---|---|---|---|
| 0 | 304.228573 | 3654.601099 | 145.483474 |
| 1 | 58.884366 | 3808.609875 | 66.642938 |
| 2 | 9.246982 | 3991.751692 | 27.603430 |
| 3 | 58.505945 | 4164.818180 | 50.832968 |
| 4 | 32.033176 | 4425.773883 | 33.290931 |
supergigantes.head()
| lum | temp | radius | |
|---|---|---|---|
| 0 | 359749.335156 | 3801.042587 | 278.055832 |
| 1 | 416869.383470 | 4398.962354 | 190.278395 |
| 2 | 1000000.000000 | 5465.163392 | 140.809113 |
| 3 | 920449.571753 | 7837.395137 | 46.187556 |
| 4 | 779830.110523 | 10200.701561 | 19.604244 |
ms.head()
| lum | temp | radius | |
|---|---|---|---|
| 0 | 0.000776 | 3577.003926 | 0.814703 |
| 1 | 0.002638 | 3691.168543 | 1.209778 |
| 2 | 0.006823 | 3793.506494 | 1.630027 |
| 3 | 0.019733 | 3862.471423 | 2.361574 |
| 4 | 0.040402 | 3963.530109 | 2.910924 |
Según el diagrama de referencia, vemos que se trata de un diagrama de dispersión donde el eje x corresponde a la temperatura y el eje y a la luminosidad. Entonces vamosa a graficar los 4 dataframes utilizando plt.pyplot.scatter, donde el eje x va a ser pd.data['temp'] y el eje y pd.data['lum'], para saber el orden de los ejes en la función pedimos información respecto a este comando, además nos proporciona ayuda respecto a las otras modificaciones que podemos hacer a la gráfica.
plt.pyplot.scatter?
plt.pyplot.scatter(enana['temp'],enana['lum'])
plt.pyplot.scatter(gigantes['temp'],gigantes['lum'])
plt.pyplot.scatter(supergigantes['temp'],supergigantes['lum'])
plt.pyplot.scatter(ms['temp'],ms['lum'])
<matplotlib.collections.PathCollection at 0x7f8db5311c50>
Claramente los ejes que estamos utilizando no permiten visualizar correctamente los puntos, es necesario que utilicemos la escala logarítmica en el eje y (plt.pyplot.yscale('log'))ya que la luminosidad de las estrellas supergigantes es mucho mayor que la de los otros grupos. Además es necesario invertir el eje x (plt.pyplot.gca().invert_xaxis())para que se asemeje más al diagrama HR.
plt.pyplot.scatter(enana['temp'],enana['lum'])
plt.pyplot.scatter(gigantes['temp'],gigantes['lum'])
plt.pyplot.scatter(supergigantes['temp'],supergigantes['lum'])
plt.pyplot.scatter(ms['temp'],ms['lum'])
plt.pyplot.yscale('log')
plt.pyplot.gca().invert_xaxis()
Aún quedan varios cambios por hacer, ahora vamos a darle tamaño a los puntos según el radio, esto se modifica en la función plt.pyplot.scatter(x,y,s=tamaño). En nuestro caso s=pd.data['radius'].
plt.pyplot.scatter(enana['temp'],enana['lum'],s=enana['radius'])
plt.pyplot.scatter(gigantes['temp'],gigantes['lum'],s=gigantes['radius'])
plt.pyplot.scatter(supergigantes['temp'],supergigantes['lum'],s=supergigantes['radius'])
plt.pyplot.scatter(ms['temp'],ms['lum'],s=ms['radius'])
plt.pyplot.yscale('log')
plt.pyplot.gca().invert_xaxis()
Ahora sigue modificar el color, vemos que todos los círculos tienen un borde gris delgado y que los marcadores de las enanas son blancos. Para esto cambiamos linewidths=0.5 y edgecolors='gray' para hacer el borde y c='white' para las enanas.
plt.pyplot.scatter(enana['temp'],enana['lum'],s=enana['radius'],c='white',linewidths=0.5,edgecolors='gray')
plt.pyplot.scatter(gigantes['temp'],gigantes['lum'],s=gigantes['radius'],linewidths=0.5,edgecolors='gray')
plt.pyplot.scatter(supergigantes['temp'],supergigantes['lum'],s=supergigantes['radius'],linewidths=0.5,edgecolors='gray')
plt.pyplot.scatter(ms['temp'],ms['lum'],s=ms['radius'],linewidths=0.5,edgecolors='gray')
plt.pyplot.yscale('log')
plt.pyplot.gca().invert_xaxis()
Para las otras estrellas, se necesita que los marcadores cambien su color según su posición en el eje x. Para esto es necesario modificar cmap en la función pyplot.scatter. En la ficha técnica de colormaps en matplotlib vemos todas las posibilidades combinaciones de colores para añadir a nuestro gráfico. Vemos entonces que plt.cm.RdYlBu es el que más se asemeja a lo que buscamos. Para aplicar el cmap correctamente es necesario normalizar la escala de colores para que 0 sea el valor más pequeño y 1 el más grande del eje x. Para eso aplicamos la función pyplot.normalize con el valor mínimo y máximo de temperatura.
normms=plt.pyplot.Normalize(ms['temp'].min(),ms['temp'].max())
plt.pyplot.scatter(enana['temp'],enana['lum'],s=enana['radius'],c='white',edgecolors='gray')
plt.pyplot.scatter(gigantes['temp'],gigantes['lum'],s=gigantes['radius'],c=gigantes['temp'],cmap=plt.cm.RdYlBu,norm=normms,linewidths=0.5,edgecolors='gray')
plt.pyplot.scatter(supergigantes['temp'],supergigantes['lum'],s=supergigantes['radius'],c=supergigantes['temp'],cmap=plt.cm.RdYlBu,norm=normms,linewidths=0.5,edgecolors='gray')
plt.pyplot.scatter(ms['temp'],ms['lum'],s=ms['radius'],c=ms['temp'],cmap=plt.cm.RdYlBu,norm=normms,linewidths=0.5,edgecolors='gray')
plt.pyplot.yscale('log')
plt.pyplot.gca().invert_xaxis()
Nos falta entonces, añadir título y etiquetar los ejes con las unidades correspondientes. Además hay que incluir texto dentro de la gráfica para identificar las zonas principales del diagrama HG. Y, finalmente aumentar el tamaño del diagrama y marcadores para que se distingan mejor.
# Normalizar el eje x para el cmap
normms=plt.pyplot.Normalize(ms['temp'].min(),ms['temp'].max())
# plot de enanas blancas
plt.pyplot.scatter(enana['temp'],enana['lum'],s=20*enana['radius'],c='white',edgecolors='gray')
# plot de estrellas gigantes
plt.pyplot.scatter(gigantes['temp'],gigantes['lum'],s=20*gigantes['radius'],c=gigantes['temp'],
cmap=plt.cm.RdYlBu,norm=normms,linewidths=0.5,edgecolors='gray')
# plot de estrellas supergigantes
plt.pyplot.scatter(supergigantes['temp'],supergigantes['lum'],s=20*supergigantes['radius'],c=supergigantes['temp'],
cmap=plt.cm.RdYlBu,norm=normms,linewidths=0.5,edgecolors='gray')
# plot de estrellas de la secuencia principal
plt.pyplot.scatter(ms['temp'],ms['lum'],s=20*ms['radius'],c=ms['temp'],
cmap=plt.cm.RdYlBu,norm=normms,linewidths=0.5,edgecolors='gray')
# Eje y en escala logarítimca
plt.pyplot.yscale('log')
# Invertir eje x
plt.pyplot.gca().invert_xaxis()
# Título del gráfico
plt.pyplot.title("Diagrama HR",fontsize=15)
# Etiquetar eje x
plt.pyplot.xlabel("Temperatura(K)",fontsize=15)
# Etiquetar eje y
plt.pyplot.ylabel("Luminosidad (L_sol)",fontsize=15)
# Añadir etiqueda 'Enanas Blancas'
plt.pyplot.text(8000,10**-3,'Enanas Blancas',fontsize=15)
# Añadir etiqueda 'Secuencia Principal'
plt.pyplot.text(8000,10**-0.8,'Secuencia Principal',fontsize=15)
# Añadir etiqueda 'Gigantes Rojas'
plt.pyplot.text(5500,10**3,'Gigantes Rojas',fontsize=15)
# Añadir etiqueda 'Enanas Azules'
plt.pyplot.text(11000,10**4.5,'Gigantes Azules',fontsize=15)
# Añadir etiqueda 'Supergigantes Rojas'
plt.pyplot.text(6400,10**4.7,'Supergigantes Rojas',fontsize=15)
# Aumentar el tamaño del gráfico y mejorar su calidad
plt.rcParams['figure.figsize'] = [20, 20]
plt.rcParams['figure.dpi'] = 100
Para animar la gráfica obtenida en el iteral anterior, es necesario definir 4 gráficas vacías que van a ir recibiendo los datos en x e y (ln1,ln2,ln3 y ln4) y, los van a ir ploteando en cada fotograma. Además se define una subgráfica que lleve las características de los ejes y etiquetas (ax). Finalmente se define una función de animación que reciba los 4 conjuntos de datos y los vaya almacenando en cada fotograma dentro de ln.
#Definimos las 4 gráficas vacías que van a recibir los datos.
fig, ax = plt.pyplot.subplots(figsize=(10,10))
ln1 = plt.pyplot.scatter([], [])
ln2 = plt.pyplot.scatter([], [])
ln3 = plt.pyplot.scatter([], [])
ln4 = plt.pyplot.scatter([], [])
# El número de fotogramas va a corresponder al número de datos del conjunto más largo en este caso (ms)
frame=len(ms['temp'])
#Ahora definimos una función que establezca las condiciones iniciales de ax, estas son simplemente una copia de los parámetros
# En las gráfica del iteral anterior
def init1():
# Establecer eje y como escala logarítmica
ax.set_yscale('log')
# Establecer los límites del eje x
ax.set_xlim(3000, 12000)
# Establecer los límites del eje y
ax.set_ylim(10**-5, 10**7)
# Invertir el eje x
ax.invert_xaxis()
# Añadir título y nombre a los ejes
ax.set_ylabel("Luminosidad (L_sol)",fontsize=15)
ax.set_xlabel("Temperatura(K)",fontsize=15)
ax.set_title("Diagrama HR",fontsize=15)
# Añadir etiqueda 'Enanas Blancas'
ax.text(8000,10**-3,'Enanas Blancas',fontsize=10)
# Añadir etiqueda 'Secuencia Principal'
ax.text(8000,10**-0.8,'Secuencia Principal',fontsize=10)
# Añadir etiqueda 'Gigantes Rojas'
ax.text(5500,10**3,'Gigantes Rojas',fontsize=10)
# Añadir etiqueda 'Enanas Azules'
ax.text(11000,10**4.5,'Gigantes Azules',fontsize=10)
# Añadir etiqueda 'Supergigantes Rojas'
ax.text(6400,10**4.7,'Supergigantes Rojas',fontsize=10)
return ln1,
def update1(frame):
data1 = np.hstack((ms['temp'][:frame,np.newaxis], ms['lum'][:frame, np.newaxis]))
ln1.set_offsets(data1)
data2 = np.hstack((enana['temp'][:frame,np.newaxis], enana['lum'][:frame, np.newaxis]))
ln2.set_offsets(data2)
data3 = np.hstack((gigantes['temp'][:frame,np.newaxis], gigantes['lum'][:frame, np.newaxis]))
ln3.set_offsets(data3)
data4 = np.hstack((supergigantes['temp'][:frame,np.newaxis], supergigantes['lum'][:frame, np.newaxis]))
ln4.set_offsets(data4)
return ln1, ln2, ln3, ln4
anim1 = matplotlib.animation.FuncAnimation(fig, update1, frames=frame,init_func=init1,blit=True)
anim1.save('HR.gif', writer='imagemagick')
MovieWriter imagemagick unavailable; using Pillow instead. /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:49: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version. Convert to a numpy array before indexing instead. /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:51: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version. Convert to a numpy array before indexing instead. /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:53: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version. Convert to a numpy array before indexing instead. /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:55: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version. Convert to a numpy array before indexing instead.
