Felipe Reyes Osorio (@reyesf)

Ejercicio 1

Para comenzar, importamos numpy, principalmente para poder usar np.loadtxt para leer los datos. También importamos matplotlib.pyplot como herramienta principal para graficar.

import numpy as np

from matplotlib.pyplot import *

Leemos todos los datos, saltando la primera línea, pues contiene los títulos, y especificando los separadores apropiados.

También definimos la función tRange para encontrar $T_{max}$ y $T_{min}$ de todo el set de datos; esto con el objetivo de calibrar la escala de colores apropiadamente. El mecanismo de la función es sencillo: toma una lista de listas que contiene las temperaturas. Un loop se encarga de hallar los máximos y mínimos de cada elemento, y al final retorna los máximos y mínimos globales.

giants = np.loadtxt("./data/giants.txt",skiprows=1)
supGiants = np.loadtxt("./data/supergiants.txt",skiprows=1)
dwarfs = np.loadtxt("./data/dwarfs.csv",skiprows=1,delimiter=",")
ms = np.loadtxt("./data/ms.csv",skiprows=1,delimiter=",")

def tRange(listOfLists):
    ans=[]
    for i in range(len(listOfLists)):
        ans.append([min(listOfLists[i]), max(listOfLists[i])])
    ans = np.array(ans)
    return min(ans[:,0]),max(ans[:,1])

bigListT=[ms[:,1],giants[:,1],supGiants[:,1],dwarfs[:,1]]
temps = tRange(bigListT)

Ahora definimos la figura y los ejes en los que se va a graficar. Cuadramos el tamallo de letra, los nombres de los ejes, la escala logarítmica para el eje y, y los rangos de los ejes. Estos últimos fueron elegidos de manera manual, experimentando con diferentes valores hasta hallar unos útiles y estéticos. Una de las características principales del diagrama de Hertzprung-Russel es que la temperatura en el eje horizontal disminuye hacia la derecha; incluimos el comando invert_xaxis() por este motivo.

fig1, ax1= subplots(figsize=(9,6))
ax1.tick_params(axis='both', which='major', labelsize=14)

vmin, vmax = tRange(bigListT)

ylabel("Luminosity [L$_\odot$]", fontsize=20)
ax1.set_yscale('log')
ax1.set_ylim(0.3e-4, 5e7)

xlabel("$T$ [K]", fontsize=20)
ax1.invert_xaxis()
ax1.set_xlim(vmax+500, vmin-900);

Ahora definimos la paleta de colores que usaremos (en este caso va del rojo al amarillo al azul), y graficamos los datos de cada tipo de estrella. Se usa scatter en lugar de plot para poder cambiar las características de los puntos. Estas son:

• c: color de los puntos. Se pasa la columna correspondiente a la temperatura. Como las estrellas más frías son más rojas, la escala de color es apropiada.
• vmin/vmax: calibra los colores. De esta manera cada tipo de estrella tiene la misma escala de color, en vez de cada una ser independiente. Si no especificamos esto, la gigante más caliente sería roja, aunque globalmente sigue siendo azul.
• s: tamaño de los puntos. Se pasa la $1.7r^{\frac{2}{3}}$. La potencia 2/3 hace que los tamaños no sean tan disparejos.
• cmap: escala de colores.
• lw: ancho de línea de los puntos. Hace que los círculos tengan borde negro.
• ec: color de los bordes.

cm = "RdYlBu"

ax1.scatter(ms[:,1],                   # Esta línea grafica las estrellas
            ms[:,0],                    # de la secuencia principal. Se 
            c=ms[:,1],                # muestran de esta manera para 
            vmin=vmin,                # mejorar la legibilidad. Las otras
            vmax=vmax,                    # estrellas sí se hacen en una sola
            s=1.7*ms[:,2]**1.5,        # línea.
            cmap=cm, 
            lw=0.8, 
            ec="k")

ax1.scatter(giants[:,1], giants[:,0], s=1.7*giants[:,2]**1.5, c=giants[:,1], vmin=vmin, vmax=vmax, cmap=cm, lw=2, ec="k")
ax1.scatter(supGiants[:,1], supGiants[:,0], s=1.7*supGiants[:,2]**1.5, c=supGiants[:,1], vmin=vmin, vmax=vmax, cmap=cm, lw=2, ec="k")
ax1.scatter(dwarfs[:,1], dwarfs[:,0], s=1.7*dwarfs[:,2]**1.5, c=dwarfs[:,1], vmin=vmin, vmax=vmax, cmap=cm, lw=0.8, ec="k");

Como las enanas resultan ser muy pequeñas en la gráfica, definimos nuevos ejes axins que servirán para hacer un zoom en esta zona. Las coordenadas son coordenadas relativas que le dicen a python donde poner los ejes dentro de la gráfica.

axins = ax1.inset_axes([0.03, 0.03, 0.38, 0.4])

Ahora, en los ejes axins graficamos los datos de las enanas de nuevo. Para que no se vea tan llena la imagen eliminamos los números de los ticks. La manera en que se hizo es un poco extraña pero en este caso fue necesario porque la escala logarítmica en un rango tan pequeño no permitió otro método.

axins.scatter(dwarfs[:,1], dwarfs[:,0], s=10*dwarfs[:,2]**1.5, c=-dwarfs[:,1], vmin=-vmax, vmax=0, cmap=cm, lw=1, ec="k")
axins.set_yscale('log')
axins.set_xticklabels("")
axins.yaxis.set_major_formatter(NullFormatter())
axins.yaxis.set_minor_formatter(NullFormatter())

Para terminar de graficar, especificamos la zona en la cual estamos haciendo zoom, y con el comando indicate_inset_zoom se agregan las líneas que lo indican. La línea comentada sirve para decargar la imagen.

axins.set_xlim(4800, 8550)
axins.set_ylim(0.00009, 0.0007)

ax1.indicate_inset_zoom(axins, lw=2);

Finalmente, ponemos labels usando el comando text. El parámetro transform es para ubicar el texto con coordenadas relativas a la gráfica, no con las coordenadas marcadas por los ejes. Las posiciones fueron escogidas manualmente.

ax1.text(0.55, 0.75, "Super-giants", transform=ax1.transAxes, fontsize=16)
ax1.text(0.85, 0.38, "Giants", transform=ax1.transAxes, fontsize=16)
ax1.text(0.35, 0.58, "Main Sequence", transform=ax1.transAxes, fontsize=16)
axins.text(0.1, 0.6, "Dwarfs", transform=axins.transAxes, fontsize=16)

fig1

fig1.savefig("fullDiagram.jpg")    #En caso de querer guardar

Ejercicio 2

Para hacer la animación primero organizamos la figura, reusando casi todo el código del Ejercicio 1. Primero volvemos a crear la figura para tener un lienzo vacío.

fig1, ax1= subplots(figsize=(9,6))
ax1.tick_params(axis='both', which='major', labelsize=14)

vmin, vmax = tRange(bigListT)

cm = "RdYlBu"

ylabel("Luminosity [L$_\odot$]", fontsize=20)
ax1.set_yscale('log')
ax1.set_ylim(0.3e-4, 5e7)

xlabel("$T$ [K]", fontsize=20)
ax1.invert_xaxis()
ax1.set_xlim(vmax+500, vmin-900);

La primera adición es esta celda. Importamos de la librería celluloid la clase Camera.

from celluloid import Camera

cam = Camera(fig1)

Ahora corremos cuatro bucles, uno para cada tipo de estrella. La idea es que graficamos una porción de los datos (usando np.delete) y tomamos una foto virtual, usando la función Camera.snap(). Luego, incluimos una nueva estrella en los datos antes de tomar la siguiente foto. En las estrellas de la secuencia principal tomamos un snap cada 3 estrellas para que la animación fuera más rápida.

dMs = len(ms[:,0])
for j in range(dMs):
    msTrash = np.delete(ms, slice(j+1, dMs-1), 0)
    ax1.scatter(msTrash[:,1], msTrash[:,0], c=msTrash[:,1], vmin=vmin, vmax=vmax, s=1.7*msTrash[:,2]**1.5, cmap=cm, lw=0.8, ec="k")
    if j%3==0: cam.snap()
        
dGs = len(giants[:,0])
for j in range(dGs):
    gsTrash = np.delete(giants, slice(j+1, dMs-1), 0)
    ax1.scatter(ms[:,1], ms[:,0], c=ms[:,1], vmin=vmin, vmax=vmax, s=1.7*ms[:,2]**1.5, cmap=cm, lw=0.8, ec="k")
    ax1.scatter(gsTrash[:,1], gsTrash[:,0], c=gsTrash[:,1], vmin=vmin, vmax=vmax, s=1.7*gsTrash[:,2]**1.5, cmap=cm, lw=2, ec="k")
    cam.snap()
    
dSg = len(supGiants[:,0])
for j in range(dSg):
    sgTrash = np.delete(supGiants, slice(j+1, dMs-1), 0)
    ax1.scatter(ms[:,1], ms[:,0], c=ms[:,1], vmin=vmin, vmax=vmax, s=1.7*ms[:,2]**1.5, cmap=cm, lw=0.8, ec="k")
    ax1.scatter(giants[:,1], giants[:,0], c=giants[:,1], vmin=vmin, vmax=vmax, s=1.7*giants[:,2]**1.5, cmap=cm, lw=2, ec="k")
    ax1.scatter(sgTrash[:,1], sgTrash[:,0], c=sgTrash[:,1], vmin=vmin, vmax=vmax, s=1.7*sgTrash[:,2]**1.5, cmap=cm, lw=2, ec="k")
    cam.snap()
    
dDw = len(dwarfs[:,0])
for j in range(dDw):
    dwTrash = np.delete(dwarfs, slice(j+1, dMs-1), 0)
    ax1.scatter(ms[:,1], ms[:,0], c=ms[:,1], vmin=vmin, vmax=vmax, s=1.7*ms[:,2]**1.5, cmap=cm, lw=0.8, ec="k")
    ax1.scatter(giants[:,1], giants[:,0], c=giants[:,1], vmin=vmin, vmax=vmax, s=1.7*giants[:,2]**1.5, cmap=cm, lw=2, ec="k")
    ax1.scatter(supGiants[:,1], supGiants[:,0], c=supGiants[:,1], vmin=vmin, vmax=vmax, s=1.7*supGiants[:,2]**1.5, cmap=cm, lw=2, ec="k")
    ax1.scatter(dwTrash[:,1], dwTrash[:,0], c=dwTrash[:,1], vmin=vmin, vmax=vmax, s=1.7*dwTrash[:,2]**1.5, cmap=cm, lw=0.8, ec="k")
    cam.snap()

Finalmente animamos todas las fotos que tomamos, y guardamos la animación a la carpeta local.

animation = cam.animate()
animation.save('./animation.gif', fps=10)

MovieWriter ffmpeg unavailable; using Pillow instead.

La animación se muestra en el notebook, siguiendo el método propuesto por @carrilloj en el mattermost. También se encuentra en el archivo animation.gif