La astronomía es sinónimo de imágenes asombrosas y siderales, literalmente fuera de este mundo. Composiciones a todo color que nos recuerdan que el universo está lleno de estructuras increibles, y conduce a la reflexión sobre nuestra realidad así como también del lugar que nos corresponde en el universo conocido. El correcto análisis de estas imágenes nos permite extraer información crítica desde la luz que emiten esos puntitos brillantes que habitan el cielo. A través de la física y sus artilugios matemáticos somos capaces de leer la luz entre lineas, para así comprender los secretos de quienes producen esta radiación: estrellas, galaxias, sus grupos...además de burbujas y marañas de gas y polvo, como las piedras que las orbitan.

La estética de la astronomía explica a sus multitudinarios y curiosos seguidores, despertando una de las dudas más comunes: ¿los colores de las imágenes son reales? ¿o son parte del malvado fotochop de la NASA con el fin de controlar nuestra vida?. Colorear imágenes es un proceso bien conocido, teniendo bases en el proceso biológico de cómo los humanos percibimos el color y cómo se entiende esta información de forma digital. Desde esta premisa, en las siguiente líneas haré el esfuerzo de unir todos estos cabos, cómo (des)componer una imagen y mostrales la forma de colorear una imagen astronómica.

El ojito humano y el color RGB.

El ojo humano es el resultado de millones de años de evolución en los cuales su sistema óptico se adaptó a la luz del Sol, nuestra bola de gas favorita. Este sistema óptico recibe y canaliza la luz hasta unas celulas especializadas fotosensibles, llamadas conos y bastones. Según wikipedia Los bastones se activan en condiciones de baja luminosidad. Los conos, por su lado, se ven estimulados en entornos con mayor iluminación y nuestros ojos están poblados por tres variedades diferentes: conos rojos, verdes y azules; sensibles a la luz visible con la cual se apellidan.

Esta sensibilidad particular a ciertos colores no es (taan) antojadiza, pasa que el sol emite muchos tipos de luz, pero la radiación que alcanza la tierra con mayor intensidad contiene estos colores, y a esto se le llama *espectro visible*, luz visible, o simplemente luz. Pero hay muchas otras variedades de luz las cuales son representadas en el rango electromagnético, que va desde las ondas de radio a los rayos gamma. Pasó que el sol emite mucha luz en el rango visible, el ojo humano se adaptó a esta radiación lo largo de su existencia y por eso vemos los colores del arcoiris. Esta es una idea muy poderosa, que nos lleva a pensar cómo se desarrollarían los detectores de radiacion (a.k.a. ojos) de seres que podrían habitan planetas con estrellas que emiten diferente luz, o también como la inexistencia de luz hace que no se desarrollen tales detectores. La siguente imágen muestra el rango electromagénico y su comparación de la luz que existe versus la luz que podemos observar, la luz visible.

Cuando la luz de una fuente luminosa incide sobre un objeto, absorbe ciertos colores y refleja otros. Es esa luz que hace rebotar un cierto objeto la que alcanza nuestros ojos y activa los diferentes conos en proporciones establecidas por los nervios ópticos, lo que nuestro cerebro interpreta como "un color". Simplificando todo, nuestro cerebro tiene 3 colores básicos y todas sus posibles proporciones para colorear nuestra realidad: rojo, verde, azul y sus conocidos mundialmente como RGB (del acrónimo inglés Red, Green, Blue). Estos colores clásicamente se conocen como colores primarios, ya que a partir de esta triada es posible obtener por completo el rango de colores del arcoiris.

Y ya que estamos aquí, es interesante señalar que de hecho vemos aún más colores que los que viven en la luz visible. El color Magenta es un claro ejemplo de esto, ya que no existe luz de este color con una longuitud de onda asociada en el mundo real. Esto significa que este color existe solo en nuestros cerebros, y es la interpretación de la superposición de dos diferentes longuitudes de onda de luz incidente. El siguiente video retrata de mejor manera lo que acabo de explicar:

Un mundo lleno de imágenes

Las imágenes son importantisimas en el día a día y ya son de uso cotidiano, al punto que ya no es necesario cuestionarse qué son en realidad y cómo pueden llegar a usarse para determinar misterios. Tan exitoso ha sido la integración de las imágenes a nuestro vertiginoso estilo de vida que muchas personas ignoran el hecho que son sólo un gran arreglo de números interpratados para que se parezcan lo más posible a lo que nosotros vemos: un instante congelado en el tiempo.

Las primeras imágenes se remontan al año 1839 gracias al daguerrotipo, el cual consistía en placas de cobre con películas de plata, que las convierte una superficie fotosensible, colectora de luz. Antes de esto había que ser buen dibujante para registrar los relieves y cráteres de la luna como hizo Galileo. Los métodos fotosensibles fueron perfecionandose a lo largo del tiempo, pasando por muchos formatos a lo largo de las décadas, como La placa fotográfica y cintas magnéticas. Hasta 1969, año en el que se desarrollaron los detectores electrónicos fotosensibles CCDs que cada uno de nosotros tiene en su teléfono, y que el 2009 (¡luego de 40 años!) les dio el novel de física a sus creadores Willard Boyle y George E. Smith.. Fue aquí donde se introdujo y popularizó el pixel como la undad básica de una imagen digital. Cada imagen, fotografía o contenido gráfico digital está hecho de pixeles, y de forma similar a los conos de nuestros ojos, cada uno de estos pixeles contiene 3 canales de color: rojo, verde y azul (o RGB). Combinando las intensidades de estos diferentes colores cada pixel puede motrar toda la gama de colores posibles, emulando así la lógica del ojo humano y creando la transición análogo/digital de la imagen.

En términos computacionales, podremos entender una imagen como un cubo de datos de (Pixeles$_{alto}\times$ Pixeles$_{largo}\times 3$). Para entrar un poco más en detalle y hagamos una biopsia a una imagen, usaremos el lenguaje de programación Python en su versión 3, que ya estamos usando pero no te diste ni cuenta. Para que funcione como corresponde, es necesario importar las funcionalidades básicas y librerías con las cuales visualizaremos todo:

Con todo esto ya funcionando, leeremos una imagen de la estructura de una nube; imagen que fue capturada mientras el sol se estaba poniendo, dándole su color rojizo:

La información de las dimensiones arriba desplegada nos dice que la imagen tiene 2848 pixeles de alto, 4248 pixeles de largo y que está compuesta por tres arreglos de datos, que son los canales de color; cada uno de estos canales posee las mismas dimensiones antes estipuladas. Es así que la imagen desplegada, para cada valor Nubes[x,y] existen 3 valores de rojo, verde, y azul, los cuales en conjunto muestran el color. En el caso anterior supongo que el color rojo domina dadas las virtudes del arrebol de la tarde. Pero no nos quedemos en supuestos y veamos la contribución de cada canal en el color de nuestra nubecita roja. Seleccionaremos todas las filas de pixeles, y 900 columnas en la parte central de la imagen, así nuestra nueva imagen Nubecita[x,y] tendrá 2848x900 pixeles y sus 3 canales RGB. Seleccionaremos cada canal de color por separado para así visualizar la contribución de cada uno de estos colores ejerce en la imagen "Original":

En el esquema anterior, las imágenes de cada canal están "coloreadas falsamente", esto quiere decir que sus colores representan un código y una escala para que podamos distinguir sus contribuciones (en este caso sería en las partes...¿amarillas, verdes?, lo más brillante), o donde no la hay (las partes azules y oscuras). Normalmente estas imágenes se les representa en escala de grises, pero bueno, ese era el color por defecto, ná que hacerle ¯\_(ツ)_/¯.

Ahora, podemos ver cuánto aporta a la imagen cada canal RGB de forma individual. Es fácil darse cuenta que la contribución del canal rojo R está en todas partes, domina las estructuras y los valores de la imagen Nubecita. Por otro lado, el canal azul mayoritariamente no aporta, y solo está presente en la parte central, donde casi no hay nubes y el cielo está de fondo. Desde aquí debo destacar dos aspectos:

1. Desde la correcta suma de las imágenes R+G+B es posible obtener la imagen Nubecita de vuelta.
2. Las imágenes de canal son ejemplos perfectos de una imagen producida por un "filtro optico", que muestra la contribución de luz en un rango definido.
3. Hacer listas es mi pasión.

Una propiedad fundamental que se desprenden de los filtros es que ven cosas diferentes al observar la misma fuente de luz. Por ejemplo si quisieramos describir la forma de la nube imagen nubecita, sería mejor estudiar el canal rojo, ya que es más sensible a la luz que se emite desde las nubes. Por otro lado, podríamos usar el filtro azul para saber dónde no hay nube, por decir algo jajaja, el filtro malo.

Es por eso que una fuente luminosa se estudia en muchos filtos, ya que la cantidad y tipo de luz de cada filtro se pueden usar para resaltar distintos elementos de un mismo fenómeno.

¡Ya tenemos nuestras primeras herramientas! ¿Para qué podrían llegar a servir? Bueno, la manipulación de imágenes es una práctica recurrente en el día de hoy y muchas son sus aplicaciones, no solo para hacer memes. Lo importante aquí en mi humilde opinión, es darse cuenta que las imagenes sólo son números en un arreglo de datos dispuestos a ser ordenados como bloques, cual Legos. Entenderlos apropiadamente nos da muchas posibilidades, por ejemplo veamos una imagen de la luna que hice hace un tiempo cuando habia una conjunción con el planeta Júpiter, que se nota como un pequeño punto brillante al lado contrario de la luna y le haremos un corte (o un "Zoom digital" que le llaman a veces) al satélite natural:

O podríamos también crear nuevas imágenes a partir de otras. En gran parte de los casos esto será posible sólo si las imágenes a combinar poseen las mismas dimensiones. Si la condición anterios se cumple, podemos llegar y sumar las imágenes diréctamente. Por ejemplo, y sólo por capricho, agregaré otra imagen de la luna en fase menguante, además de rotar la imagen 1 en 180 grados, y sumaremos todo en una gran imagen. No se diga más.

Bueno, hay algunos detalles de los cuales mejorar en esta composición; como esos puntos blancos sobre la luna creciente que posiblemente están sobre-expuesta dada la combinación de las 3 imagenes. Pero en realidad da igual, programas de edicion de datos como Photoshop o Gimp se especializa en esto y sin tener que saber tanto sobre la composicion, estructura de imágenes, ni de filtros gaussianos, etc.

Pero yo les prometí imágenes astronómicas, además en astronomía/astrofísica, a casi a nadie le importa la luna. Es más, nos cuidamos de su prescensia y apuntamos los telescopios cuando no está para que su terrible luz no opaque los escasos fotones que provienen de nuestras fuentes astronómicas favoritas. Esto es todo lo que tengo que decir con respecto a la luna.

El color de las imágenes astronómicas

Por si no lo sabían, la cantidad de información astronómica disponible en este momento es brutal, onda en el orden de Petabytes de información (1.000.000.000.000.000 bytes) se producen a diario, y en el futuro cercano esta tasa irá al alza, o bueno esa era la proyección antes de la pandemia y la desolación mundial del 2020. Mucha de estas imágenes hecha por telescopios en Tierra y el espacio, como la información que ha sido obtenido a lo largo de su análisis ya es de libre acceso, y algunos sets de datos están especialmente orientados a la divulgación científica, como las históricas imágenes de la Nebulosa del Aguila.

Hechas con la cámara WFPC2 montada en el telescopio espacial Hubble, las imágenes alojadas en la página de la NASA muestran los llamados "Pilares de la creación", las cuales componen parte de la nebulosa del Aguila, y son columnas de polvo desde donde la gravedad hace lo suyo y ensambla nuevas estrellas. Ya el hecho de que tenga nombre propio habla de la popularidad e impacto que tuvo esta imágen alrededor del mundo. Formalmente a este tipo de estructuras se les conoce como "pilares moleculares fríos", como también "*trompas de elefante*".

Estas imágenes están hechas usando filtros de banda estrecha (narrow filters), los cuales se encargan de medir la luz emitida por un elemento en particular, el cual emite radiación en una longuitud de onda específica del rango electromagnético. Las imágenes disponibles están en formato FITS, que es el formato por defecto en la astronomía actualmente.

Entonces, las imágenes disponibles para componer los pilares de la creación son:

1. Filtro estrecho de Oxígeno doblemente ionizado (OIII), medido a 502 nanómetros.
2. Filtro estrecho de Hidrógeno Alfa (H$_\alpha$), la cual posee $\lambda=656$ nanómetros, por supuesto.
3. Y por ultimo, filtro estrecho de Azufre ionizado (SII), medido a 673 nanómetros

Los elementos anteriores se encuentran presentes en el medio estelar, y son emitidos por la composición e interacción de nebulosas como la que nos interesa. Entonces, descargaremos las imágenes, las leeremos de la misma forma como lo hicimos hace unos momentos más arriba, y luego usaremos la funcion sqrt disponible en astrobetter para ayudar en la visualización de las imágenes:

Es posible ver las diferencias, sin dudas (¿verdad?). Por un lado la imagen de AzufreII remarca mucho más los bordes de las nubes, dado que este elemento se encuentra cuando gases a grandes velocidades colisionan. Por otro lado el OxigenoIII y H$_\alpha$ son elementos que se emiten y mapean la materia que está en nuestra querida nebulosa. El OxigenoIII también se ve difuso y marca mucho la absorción de dos de los tres pilares.

Ahora ¿Podemos crear una imágen con estas tres imágenes? Por supuesto, pero será una IMAGEN DE COLOR FALSO. Lo que quiere decir que imágenes realizadas con filtros que no tienen que ver con un color en sí, serán asignados al canal de color para que nuestros ojos puedan entender sus caracteristicas y complejas composiciones fuera de nuestro rango de visión.

Con esto en mente, usaremos python para crear un cubo de datos, y asignaremos cada imagen de filtro estrecho a un canal de color RGB: rojo al mapeo del AzufreII, verde para la emisión de Hidrógeno alfa y el canal azul será para la luz del Oxígeno doblemente ionizado. Con esto, podemos obtener la siguiente imagen.

Desde el punto de vista astronómico, la Fotometría es la medición, calculo y cuantificación de la luz emitida por diferentes objetos astronómicos. Etimilogicamente, proviene de la voz griega “φωτος” (photos,luz) y del sufijo “μετρια” (metria, medida).

De esta forma, la fotometría es la medición del flujo integrados $F$ (cuentas/tiempo x area) de una fuente en el cielo, y se define matemáticamente de la siguiente forma:

Donde $F_{\lambda}$ es el flujo específico proveniente de la fuente, y $S_{\lambda}$ la función de transmisión, que toma en cuenta la respuesta del detector, la condiciones atmosféricas y el filtro que se ha usado. Lógicamente, y por cómo está planteada la medición, mientras mayor sea el tiempo $t$, mayor flujo integrado se obtiene.

Magnitudes y flujos

La magnitud astronómica es, por definición, una cantidad relativa. Esto quiere decir que su valor es definido con respecto a una fuente específica (como por ejemplo, el sistema fotométrico basado en la estrellas Vega). Dos variaciones populares de la magnitud astronomómica son las siguientes:

• Magnitud Bolométrica: $$ \displaystyle m_1 - m_2 = 2.5\log_{10}\left(\frac{F_2}{F_1}\right) $$

• Magnitud específica: $$ m_{\lambda,1} - m_{\lambda,2} = 2.5\log_{10}\left(\frac{F_{\lambda,2}}{F_{\lambda,1}}\right) $$

Entonce, necesitamos los Flujos de, al menos, dos estrellas relativamente cercana. Hay tres métodos principales para medir dicho flujo $F$:

1. Fotometría de apertura.
2. Fotometría diferencial.
3. Fotometría de función de dispersión de punto, o point-spread-function (PSF).

Cada una de estas técnicas posee diferentes filosofías, ventajas y desventas, pero todos tienen en común algunos parámetros. Partamos con la primera:

Fotometría de apertura

Esta técnica se basa en considerar todas las cuentas (ADU) dentro de un radio $R$ definido, el cual es llamado "apertura", teniendo como regla que el radio de apertura debe ser, a lo menos, menor que el seeing $R<seeing$ (visión o visibilidad astronómica).

Hay que recordar que una imagen (astronómica, fotográfica, etc) es un arreglo de datos de $(n \times n)$ píxeles, entonces los valores $I_{ij}$ describiran las coordenadas físicas de la imagen. Naturalmente, el flujo de luz proveniente de una fuente estelar, por un tiempo de antemano conocido, puede ser estimado sumando cada pixel $I_{ij}$ dentro de la circunferencia caracterizada por el radio R. En palabras más precisas, dentro de la apertura.

Luego, debemos sustraer la contribución del cielo a la medición usando alguna estimación directa del valor de "cielo". Podemos pensar en esta contribución como el ruido de fondo de la imagen, y será cuantificado como el valor promedio de una imagen en lugares donde no existen objetos astronómicos. Para estimar el cielo, se define una región llamada "Annulus", que es una estructura creada con dos circulos que posee forma de donut. La idea es medir la contribución local (el área circundante) de la imagen para obtener el flujo usando algún estadístico, digamos, el promedio.

La Figura 1 representa de mejor manera este concepto. Por ejemplo, medimos todos los píxeles $I_{ij}$ dentro del radio interior (Inner Annulus marcado en la imagen) y el exterior (Outer Annulus), obtenemos un valor promedio y será considerado el cielo.

Bajo esta interpretación, el número total de cuentas $I$ dentro del radio de apertura $R$ será la suma de cada uno de los píxeles $I_{ij}$ de la imagen dentro del radio R considerado, menos la contribución del cielo $I_{cielo}$ en cada píxel:

donde:

1. $I$: Cuentas totales de la fuente dentro del radio de apertura.
2. $I_{ij}$: Cuentas específicas en cada pixel dentro del radio de apertura.
3. $n_{pixeles}$: Número de píxeles dentro de la apertura.
4. $I_{cielo}$: Valor de la contribución estimada del cielo usando el annulus. Usualmente es llamado 'Cielo'.

Generalmente, el valor del cielo se obtiene desde la distribución del número de cuentas dentro de la imagen.

Finalmente, la magnitud instrumental viene dado por:

El termino constante depende de las variables asociadas a las condiciones de observación (por ejemplo la extincion atmosferica, la masa de aire, el seeing), y básicamente es lo que se modifica para poder calibrar la fotometría en cierca escala de interés.

La gran desventaja de este método es la alta densidad estelar. No es posible usar fotometría de apertura cuando la fuente bajo estudio no está aislada, que es el caso de grupos de estrellas como los Cúmulos globulares, o en el centro de nuestra galaxia. De alguna forma, el veloz avance de la resolución espacial de las nuevas generaciones de telescopios ha demostrado que la fotometría de apertura es la primera aproximación a la medición del flujo estelar (o gáctico).

Fotometría Diferencial

La fotometría diferencial consiste en usar estrellas de comparación dentro de la misma imagen para estudiar una estrella en particular, y así ignorar hasta cierto punto las variables atmosféricas y el tiempo de integración. Este método está enfocado en cuantificar la variabilidad estelar de la fuente de interés, y para ser implementada es necesario elegir unas 4 o 10 estrellas de comparación C, medir su magnitud instrumental $m_C$ a través de varias observaciones (o épocas) y calcular el diferencial de magnitud $\Delta m_C$ de cada estrela de comparación, con el fin de cuantificar el cambio de brillo de la estrella de interés.

Este método ofrece muy buenos resultados cuando la cadencia de observación es alta, siendo uno de los favoritos del AAVSO, Kepler, entre otros, pero desconozco si es un método que se siga usando de forma masiva, pero es una interesante forma de medir la luz de las estrellas.

Fotometría PSF

Este método se basa en la idea que todas las fuentes en una imagen tienen exactamente la misma forma. En otras palabras, la distribución espacial de la luz capturada sobre una imagen sigue alguna regla general, que puede ser usada para estimar el flujo F.

PSF es un acrónimo en ingles Point spread functionLa aproximación más usual es asumir que las fuentes poseen un perfil Gaussiano con la misma desviación estándar. Con esta aproximación, el "tamaño" o apertura para una fuente estará caracterizada por el seeing atmosférico. Dependiendo del tamaño de la imagen en la cual trabajemos, podemos considerar que tenemos, en promedio, las mismas perturbaciones atmosféricas. Bajo estas consideraciones, todas las estrellas deberían ser descritas por la misma función puntual de dispersión (point spread function, PSF).

La forma radial de una Gaussiana bidimensional se describe como:

donde:

1. $I(r)$ Es la intensidad para algún pixel de coordenadas $(x,y)$.
2. $r=\sqrt{x^2+y^2}$ es la componente radial de una Gaussiana 2-D.
3. $I_{0}$ Es el peak de intensidad del perfil.
4. $\sigma$ Es la desviación estándar del perfil (el ancho del perfil Gaussiano).

La simpleza y utilidad de la fotometría PSF recae en que todas las estrellas pueden ser modeladas usando, en este caso, un perfil Gaussiano con la misma desviación estándar (seeing) $\sigma$, y cada una con diferentes intensidades $I_0$, lo que reduce el problema a elegir un set de parámetros para cada imagen a analizar.

A continuación, se presentan 3 estrellas de diferentes magnitudes, pero que presentan la misma forma:

¿Cuál es la diferencia? La forma de la fuente se mantiene, pero el nivel de ruido del fondo crece. En general, si el peak de intensidad $I_0$ está por sobre el límite del ruido, podrá ser detectada y analizada.

Para entender mejor la idea, simularemos una estrella usando una distribución gaussiana y estudiaremos sus parámetros con más detalles. Le añadiremos ruido para hacer todo un poco más realista.

#Agregamos las librerias fundamentales para trabajar.
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import os
import sys

#Creamos el dominio:
Rango_imagen=100
numero_filas    = 100
numero_columnas = 100

x   = np.linspace(-Rango_imagen, Rango_imagen, numero_columnas)
y   = np.linspace(-Rango_imagen, Rango_imagen, numero_filas)
X,Y = np.meshgrid(x, y)

#calculo de la intensidad
r0 = 0                                    #Desplazamiento
sigma = 20.                               #Desviación estándar.
I0= 1000.                                 #Peack de intensidad.
R = np.sqrt((X-r0)**2 + (Y-r0)**2)        #Componente radial.

#Calculamos el Perfil.
Perfil = I0*np.exp(-0.5*(R/sigma)**2)     

#Simulamos ruido que emule la contribución de fondo:
Valor_piso = 100.
desviacion_estandar_ruido = 150
Dimension_imagen = Perfil.shape
background = np.random.normal(Valor_piso, desviacion_estandar_ruido, Dimension_imagen)

#componemos todo:
Perfil_ruido = Perfil + background

#Hacemos el plot:
plt.figure(figsize=(12,6))

plt.subplot(121)
plt.title('Perfil Gaussiano')
plt.imshow(Perfil, cmap=cm.jet)
plt.colorbar(orientation='vertical', fraction=0.046) 
plt.xlabel('pixel X',fontsize=15)
plt.ylabel('pixel Y',fontsize=15)

plt.subplot(122)
plt.title('Perfil Gaussiano con ruido')
plt.imshow(Perfil_ruido, cmap=cm.jet)
cbar = plt.colorbar(orientation='vertical', fraction=0.046) 
cbar.ax.get_yaxis().labelpad = 20
cbar.ax.set_ylabel('Intensidad $I_0$', rotation=270, fontsize=15)
plt.xlabel('pixel X',fontsize=15)
plt.show()

#Creamos dos perfiles con la misma desviación estándar y diferente intensidad:
centro_imagen = len(X[:,0])/2

I1 = I0
I2 = 2700.                          #Peack de intensidad.

#Creamos un segundo perfil para comparar:
Perfil_2 = I2*np.exp(-0.5*(R/sigma)**2)     

#Simulamos (otro) ruido que emule la contribución de fondo:
background = np.random.normal(Valor_piso, desviacion_estandar_ruido, Dimension_imagen)

#componemos el segundo perfil:
Perfil_ruido_2 = Perfil_2 + background

plt.figure(figsize=(10,7))
plt.plot(X[centro_imagen],Perfil_ruido[centro_imagen],'ro',ms=7,label='$I_1 = %1.0f$' %I1)
plt.plot(X[centro_imagen],Perfil_ruido_2[centro_imagen],'kH',ms=6,label='$I_2 = %1.0f$' %I2)
plt.xlabel('pixel X',fontsize=15)
plt.ylabel(r'Intensidad',fontsize=15)
plt.plot([min(X[centro_imagen]),max(X[centro_imagen])], Valor_piso*np.ones(2),'g-',lw=4,label="Piso constante")
plt.legend(fontsize=15)
plt.show()

#Importamos todas las librerías necesarias para nuestro propósito:
from photutils.detection import IRAFStarFinder
from photutils.psf import IntegratedGaussianPRF, DAOGroup
from photutils.background import MMMBackground, MADStdBackgroundRMS
from astropy.modeling.fitting import LevMarLSQFitter
from astropy.stats import gaussian_sigma_to_fwhm

from astropy.table import Table
from photutils.datasets import (make_random_gaussians_table,
                                make_noise_image,
                                make_gaussian_sources_image)

#Ajustamos algunos valores para empezar
sigma_psf = 3.0
tshape = (32, 32)   #Dimensiones de la imagen

sources = Table()   #crea una tabla para mapear los parametros de 3 fuentes
sources['flux'] = [700, 1900, 1000]
sources['x_mean'] = [14, 16, 16]
sources['y_mean'] = [15, 15, 17]
sources['x_stddev'] = sigma_psf*np.ones(3)
sources['y_stddev'] = sources['x_stddev']
sources['theta'] = [0, 0, 0]
sources['id'] = [1, 2, 4]

#Crea la imagen, que será la suma de 3 fuentes, 
#y agregando ruido basado en una distribución de Poisson
# y ruido normal:
image = (make_gaussian_sources_image(tshape, sources) +
         make_noise_image(tshape, type='poisson', mean=8.,
                          random_state=1) +
         make_noise_image(tshape, type='gaussian', mean=0.,
                          stddev=5., random_state=1))

#ploteamos la imagen resultante:
plt.figure(figsize=(5,6))
plt.imshow(image, cmap='viridis', aspect=1, interpolation='nearest') 
plt.title('Imagen simulada') 
plt.colorbar(orientation='vertical', fraction=0.046, pad=0.04) 
plt.show()

#importamos más funciones especificas:
from photutils.detection import IRAFStarFinder
from photutils.psf import IntegratedGaussianPRF, DAOGroup
from photutils.background import MMMBackground, MADStdBackgroundRMS
from astropy.modeling.fitting import LevMarLSQFitter
from astropy.stats import gaussian_sigma_to_fwhm
from photutils.psf import IterativelySubtractedPSFPhotometry

#Medimos el cielo usando el estadistico MAD:
bkgrms = MADStdBackgroundRMS()
#usamos la imagen creada como dato de entrada:
std = bkgrms(image)

#Método para encontrar iterativamente fuentes en la imagen.
iraffind  = IRAFStarFinder(threshold=3.5*std, \
            fwhm=sigma_psf*gaussian_sigma_to_fwhm, \
            minsep_fwhm=0.01, roundhi=5.0, roundlo=-5.0, \
            sharplo=0.0, sharphi=2.0)

daogroup = DAOGroup(2.0*sigma_psf*gaussian_sigma_to_fwhm)

#estimamos la contribucion del cielo (background):
mmm_bkg = MMMBackground()

#eligen el modelo PSF:
psf_model = IntegratedGaussianPRF(sigma=sigma_psf)

#Creamos el modelo a aplicar en nuestra imagen:
photometry = IterativelySubtractedPSFPhotometry(finder=iraffind,
                                                group_maker=daogroup,
                                                bkg_estimator=mmm_bkg,
                                                psf_model=psf_model,
                                                niters=1, fitshape=(11,11))
#Realizamos la fotometría:
resultados = photometry(image=image)

#revisamos los parámetros principales de la fotometría, el flujo:
list(result_tab)

[<Row index=0 masked=True>
       flux_0            x_fit               x_0               y_fit               y_0              flux_fit        id  group_id      flux_unc           x_0_unc             y_0_unc       iter_detected
      float64           float64            float64            float64            float64            float64       int64  int64        float64            float64             float64           int32    
 ----------------- ------------------ ------------------ ------------------ ------------------ ------------------ ----- -------- ----------------- ------------------- ------------------- -------------
 2320.577813609559 15.489564216228956 15.531749294180742 15.775939905519783 15.318938753350247 2451.8336398655583     1        1 37.86912299152522 0.06859803592598557 0.06901161861692053             1]

#obtenemos la imagen residual de la ejecución anterior.
residual_image = photometry.get_residual_image()

plt.figure(figsize=(9,6))

plt.subplot(121)
plt.imshow(image, cmap='viridis', aspect=1, interpolation='nearest')
plt.title('Imagen simulada')

plt.subplot(122)
plt.imshow(residual_image, cmap='viridis', aspect=1,
           interpolation='nearest', origin='lower')
plt.title('Imagen residual')
plt.show()

En el 2017, un objeto extraño y sin precedentes fue ubicado por el telescopio Pan-STARRS (localizado en Hawaii) entre las miles de millones de fuentes estelares que habitan en nuestra galaxia. En el momento del hallazgo, esta cuerpo fue catalogado premilimarmente como un cometa, dado que cumplía alguno de los parámetros esperados para este tipo de objetos celestes; pero luego de la revisión rutinaria se empezaron a develar las exóticas propiedades de este cuerpo: Su brillo cambiaba a lo largo del tiempo, una velocidad abrumadora y una trayectoria completamente inusual, tal que no parecía orbitar nuestro sol en absoluto. Por toda la expectación que esto generó, el objeto en cuestión tenía más un perfil de ser un asteroide, y fue bautizado "Oumuamua", que significa "Explorador" en Hawaiano.

Un objeto fuera del sistema solar.

La velocidad y trayectoria con la que Oumuamua atravesaba el sistema solar (velocidad estimada de 315.431 $\frac{Km}{hora}$ en su punto más cercano al sol [1]) sugería que no poseía una ligadura gravitacional a nuestro sistema solar, por lo que debía provenir desde otro sistema estelar ajeno al nuestro, algo sin precedentes hasta la fecha. Además, al terminar la visita de nuestro vecindario cósmico, Oumuamua seguiría su trayectoria a lo largo de el espacio interestelar, tomando un nuevo rumbo gracias a modificación de la trayectoria causada por la cercanía a nuestra estrella huésped. El siguiente video ilustra la trayectoria calculada del misterioso objeto, y se compara con el movimiento observado. También es posible ver la dirección desde dónde provino y hacia dónde fue.

Pero eso no lo es todo, había otro dato observacional del cual podía extraerse información valiosa: El asteroide poseía una fluctuación periódica de su brillo superficial, es decir, que su luz cambiaba siguiendo un ciclo aparentemente definido, algo que es propio de las estrellas y cuerpos suficientemente energéticos para producir luz propia. Como un objeto de esta clase no posee las condiciones mínimas para generar reacciones nucleares en su interior, la variabilidad de su brillo debe producirse por otro tipo de características, como la forma del mismo.

Un asteroide con forma...¿de cigarro?

Para internarse en el problema de su brillo, observaciones de diferentes telescopios fueron combinadas para dar claridad en los detalles. Varios equipos de astrónomos usaron estos datos y encontraron periodo de fluctuación de brillo de este "explorador" era de algunas horas (entre 6.9 a 8.3 horas)[2], lo que podía explicarse si es que el asteroide estuviera girando sobre si mismo en una especie de tambaleo. Una de las consecuencias de esa afirmación, es que Oumuamua poseía una forma preferentemente cilíndrica y alargada, con algunos cientos de metros de largo. Una forma como la señalada es enormemente inusual para un asteroide. El hecho de sugerir esta forma, y dados todos los otros detalles antes señalados, encendió rápidamente las alarmas y la imaginación, dado que anteriormente equipos de trabajo que han desarrollado investigación en transporte espacial de largas distancias habían sugerido que la forma de aguja o "cigarro" era la más adecuada para una nave interestelar, dado que la geometría minimizaría la fricción y daños provocados por el polvo, gas u otras estructuras interestelares. Esta geometría para una nave es tan intuitiva que se ve reflejada en las naves espaciales de algunas sagas de ciencia ficción. Muchas comparaciones fueron hechas, y acá hay un ejemplo de esto. Con todo esto en mente, la primera impresión que uno podría tener de Oumuamua se refleja en la Figura 1.

Aún así, debemos considerar que la posibilidad de que un objeto tenga una forma inusual, como por ejemplo la luna "Pan" de Saturno, que posee la peculiar forma de un 'ravioli' (Figura 2) no es desestimable. Varios astrónomos han señalado diferentes escenarios donde se podría obtener estas particulares formas, como la colisión de dos o más asteroides. Aún así, la idea de que pudiera ser un artefacto alienigena es atractiva y, por ende, difícil de quitar del imaginario colectivo.

Sin embargo, nada más en en sistema solar se ve de esta forma alargada cilíndrica, y mientras pasaban los días, más astrónomos reportaban diferentes periodos de cambio de luminosidad, lo que podía sugerir que si estaba efectivamente rotando sobre un eje, quizá había más de un eje de rotación[3], lo que producía un periodo más caótico. El hecho de asumir más de un eje de rotación nos conduce a volver a pensar cómo sería la forma de Oumuamua, y se piensa que la siguiente representación artística sea más acertada:

Buen viaje, Oumuamua.

En septiembre del 2017, el asteroide pasó por su punto más cercano al Sol, pasando por dentro de la órbita de Mercurio; y emprendió el viaje de retirada hacia otros horizontes de la galaxia. En enero de este año (2019) el asteroide ya sobrepasó la órbita de Saturno, y le tomará otros varios miles de años para escapar completamente del sistema solar. En su corta estadía en la parte interna de nuestro pequeño sistema planetario, Oumuamua ha dejado muchas más dudas que certidumbre, iluminando nuevamente los ojos de aquellos que buscan evidencia de vida extraterrestre. Aunque gran parte de su fenomenología ha conseguido ser explicada a través de fenómenos físicos que han sido observados en otros cuerpos, Oumuamua es un objeto asombroso y único en su especie, por lo que quedamos con la expectativa de que otro objeto similar pueda ser visualizado en una época temprana para su completo estudio, entendimiento y caracterización.

Referencias:

• Breakthrough listen to oberve interestellar object ‘OUMUAMUA: https://breakthroughinitiatives.org/news/14 [1]
• What and Whence 1I/‘Oumuamua: A Contact Binary from the Debris of a Young Planetary System? https://arxiv.org/pdf/1712.06721.pdf [2]
• The tumbling rotational state of 1I/‘Oumuamua https://arxiv.org/pdf/1711.11530.pdf [3]

Hace bastante tiempo quería explicar el por qué creo que programar es importante como una manera de estructurar e implementar el pensamiento lógico. Comprender sus aplicaciones nos habilita a entender diferentes aspectos de nuestra realidad. Con sus herramientas uno podría, por sólo dar ejemplos, construir y reproducir experimentos científicos, estudiar bases de datos para llegar a conclusiones y exhibirlas de una manera comprensible.

Dominar en algún nivel sus herramienta es fundamental en la actualidad, no sólo para las personas relacionadas a la ciencias e ingeniería sino también para aquellos que se interesan en la búsqueda, visualización e interpretación de diversas bases de datos, como periodistas y/o psicólogos.

Creo que al adquirir el punto de vista de la implementación y flujo de cálculos matemáticos, tomé perspectiva del potencial de programar un sistema. Además si nos ponemos románticos, podríamos decir que programar es como crear un pequeño mundo, el cual sigue las reglas asignadas por el escritor sin excepción alguna.

Gran parte de estas técnicas de análisis se desarrollaron originalmente usando lógica y cálculos matemáticos complejos, los cuales fueron potenciados con la aparición de computadores capaces de optimizar operaciones, automatizar procesos y realizar cálculos sostenidos a través del tiempo. Hoy en día, el acceso a este tipo de análisis es bastante transversal en comparación al siglo pasado, básicamente es necesario una computadora personal o smartphone, acceso a internet, los datos a analizar y simpatía por el idioma inglés.

Hay toda clase de tutoriales en internet con los cuales uno puede aprender de manera autodidacta sobre variadas temáticas, lenguajes de programación y futuras aplicaciones. Yo aquí recolectaré y daré un espacio a lo interesante que todo esto podría llegar a ser.

Python como primer lenguaje¶

Como sistema de comunicación, la lógica computacional puede ser descrita y comunicada a través de diferentes idiomas, los cuales son denominados lenguajes de programación. El criterio para elegir uno de ellos pueden ser variados, como la velocidad de ejecución u entorno de desarrollo.

Python es un lenguaje apropiado para aprender de cero a implementar código de manera simple, ya que posee una sintaxis simple y directa lo que conlleva una fácil interpretación para el usuario. Varías estadísticas revelan que Python ha crecido exponencialmente en términos de usuarios en los últimos años. Esto alienta a la creación de librerías especializadas y de carácter libre.

Daré algunos ejemplos de lo que se puede hacer:

El atractor de Lorenz:¶

El sistema de Lorenz es un set de ecuaciones diferenciales acopladas las cuales tiene soluciones caóticas dados ciertos parámetros $\sigma$, $\rho$ y $\beta$ con las condiciones iniciales $x(0)$, $y(0)$ y $z(0)$ evaluadas en la posición inicial de nuestra partícula P.

$$\displaystyle \frac{\rm dx}{\rm dt} = \sigma(y - x)$$

$$\displaystyle \frac{\rm dy}{\rm dt} = x(\rho - z) - y$$

$$\displaystyle \frac{\rm dz}{\rm dt} = xy - \beta z$$

El truco acá será discretizar los diferenciales del set de ecuaciones, asumiendo que el infinitesimal temporal ${\rm d}t$ puede ser aproximada usando una diferencia finita $\Delta t$ suficientemente pequeña. Lo mismo será asumido para las variables espaciales, así por ejemplo para la variable i-ésima $x_{i}$, tenemos que $\Delta x_{i+1} = x_{i+1} - x_{i}$.

Dicho esto, nuestras ecuaciones toman la siguiente forma:

$${\rm dx} = \sigma(y - x){\rm dt}\ \ \ \ \ \ \ \Rightarrow\ \ \ \ \ \ \ x_{i+1}=x_{i}+\sigma (y_{i}-x_{i}) \Delta t$$

$${\rm dy} = (\rho x - z x - y){\rm dt}\ \ \ \ \ \ \ \ \ \Rightarrow\ \ \ \ \ \ \ \ y_{i+1}=y_{i}+(\rho x_{i} - z_{i} x_{i} - y_{i}) \Delta t$$

$${\rm dz} = (xy - \beta z){\rm dt}\ \ \ \ \ \ \ \ \Rightarrow\ \ \ \ \ \ \ \ z_{i+1}=z_{i}+(x_{i}y_{i}-\beta z_{i})\Delta t$$

Así, pondremos una partícula de prueba en este mundo gobernado por estas tres leyes y evaluaremos su comportamiento usando 50 mil pequeños pasos $\Delta t$ para ver cómo evoluciona a través de 50 segundos.

No se diga más!

Entonces, la partícula roja se ha desplazado acorde a las leyes de Lorenz, y ha trazado el comportamiento característico de este set de ecuaciones diferenciales, y por lo cual es frecuentemente llamado el Atractor de Lorenz. El comportamiento individual de cada variable a través del tiempo puede ser consultado en los gráficos en la parte superior.

Un mapa sismológico de Sudamérica¶

Si vives en Chile, probablemente te habrás dado cuenta que es un país muy sísmico. Tan solo en ésta década hemos tenidos varios terremotos con una magnitud sobre 7 en escala de magnitud local. Cada vez que hay un temblor, los memes vuelan e imploran que se presente Marcelo Lagos a dar el reporte, y hasta un trago se llama Terremoto. Entonces: ¿será verdad que en otras partes del mundo, esto de los sismos no es normal?

Esta pregunta podría ser resuelta si es que visualizamos cada evento registrado desde alguna fecha arbitraria y ver cómo estos se localizan a través de Sudamérica. Para saber esto, hemos consultado la base de datos de la pagina de monitoreo sismológico de la Universidad de Chile, entre los años 2003 y 2015, obteniendo una lista con más de 58 mil sismos distribuidos a lo largo del continente.

Estos datos han sido recolectados en este post (usando otro lenguaje de programación :p) y crearemos un mapa con la distribución de sismos en Sudamérica entre el 2003 y 2015.

Así nomás.

Con todo esto, podemos ver que 8 de cada 10 sismos, prefieren Chile :p

Esto es sólo una muestra de todo lo que es posible hacer cuando tienes datos disponibles, herramientas para su análisis, y tiempo para aprender.

Referencias:¶

1) Sobre el atractor de Lorenz: https://scipython.com/blog/the-lorenz-attractor/

2) Datos sismológicos: http://benjad.github.io/2015/08/21/base-de-datos-sismos-chile/

3) Basemaps y los plots de mapas: https://jakevdp.github.io/PythonDataScienceHandbook/04.13-geographic-data-with-basemap.html