dic 11 2007

Espectropolarimetría con Python

Tag: Herramientas y Programaciónfranwerst @ 7:21

Aparece hoy en arXiv.org un documento que se publicará en el libro de proceedings del congreso "Solar Polarization Workshop #5" (editores S. Berdyugina, K.N. Nagendra, R. Ramelli). Se trata de un artículo de F. Paletou, R. Rezaei y L. Léger en el que se muestra un ejemplo de reducción y análisis de datos espectropolarimétricos en el campo de la física solar. Los datos que sirven para el ejemplo fueron tomados con el telescopio franco-italiano THéMIS, en el Observatorio del Teide, Tenerife. La novedad es que se utiliza para ello software completamente basado en Python y algunas de sus librerías científicas, como PyFITS, NumPy/Scipy y matplotlib.

Reproduzco aquí el resúmen del documento:

 

Most of the solar physicists use very expensive software for data reduction and visualization. We present hereafter a reliable freeware solution based on the Python language. This is made possible by the association of the latter with a small set of additional libraries developed in the scientific community. It provides then a very powerful and economical alternative to other interactive data languages. Although it can also be used for any kind of post-processing of data, we demonstrate the capabities of such a set of freeware tools using THeMIS observations of the second solar spectrum.

Algunos enlaces relacionados:


dic 06 2007

Estrellas que fuman

Tag: Astronomía y Astrofísicajapp @ 12:44

Existen muchos tipos de estrellas variables, según el origen de variación de brillo y de la manera en que cambia y generalmente están relacionados con cosas que pasan en el interior de la propia estrella, como ocurre con las variables cefeídas o por motivos externos, como las variables tipo Algol, que son binarias eclipsantes. Las R Corona Borealis (R CrB), son un tipo de variable muy peculiar, estrellas gigantes unas 50 veces mayores que el Sol, que varían de forma irregular y repentina, disminuyendo varias magnitudes en pocas semanas y recuperándose en unos meses o años.

Siempre se pensó que esa caída de brillo se debe a polvo liberado por la propia estrella, que oculta su luz en la línea de visión hacia la Tierra y que luego se dispersa poco a poco volviendo a su brillo natural. En realidad, no se tenía ninguna prueba directa de esta explicación hasta hace poco, cuando los astrónomos Patrick de Laverny y Djamel Mékarni utilizaron una cámara infrarroja con óptica adaptativa (NaCo, NAOS-CONICA) del telescopio VLT de la ESO, para observar la estrella RY Sagittarii. Haciendo interferometría con dos de los cuatro telescopios VLT (de 8.2m cada uno), exploraron una zona a 110 unidades astronómicas de la estrella donde observaron el envoltorio de polvo que la rodea; además, fueron capaces de ver las variaciones de las nubes polvo en pocos meses.

RY Sgt - NaCo/VLT

Este descubrimiento confirma la teoría de porqué las estrellas R CrB varían, aunque aún quedan muchos de detalles sobre ellas que resolver, como dónde se forman las nubes de polvo.


dic 03 2007

Medir sobre figuras con python

Tag: Herramientas y Programaciónahc @ 11:42

Actualización: Hay una nueva versión de Srip disponible en la página web del proyecto.

Cuando en una publicación -del tipo que sea- queremos mostrar un conjunto de datos numéricos que guardan cierta relación entre sí, la forma habitual de hacerlo es mediante una tabla o un gráfico. A menudo, el gráfico es preferible porque la distribución de los valores salta a la vista, y permite una comprensión cualitativa de la información representada mucho más rápida y eficiente que la lectura de valores en una tabla. Por otro lado, cualquier análisis cuantitativo requerirá disponer de los valores numéricos tabulados.

 

Representar los datos gráficamente, si se dispone de una tabla, es una tarea trivial y rápida si se utiliza un ordenador. El proceso inverso -el de recuperar los valores numéricos representados en una figura- suele ser tedioso e impreciso, porque en general el ordenador no puede hacerlo y es necesario que el usuario lea los valores manualmente, uno por uno, y a menudo con la dificultad añadida de escalas logarítmicas en los ejes o pocas marcas de referencia.
captura de sripPor desgracia, es muy habitual que necesitemos utilizar información numérica que sólo está disponible en forma de gráficos. En astrofísica, casos comunes son las curvas de transmisión de filtros, de eficiencia de detectores, espectros, cartas de localización de astros, diagramas color-magnitud, etc etc.

 

Para facilitar la tarea del investigador que necesita medir valores sobre un gráfico, he escrito un programa en python que, aunque aún está en fase de desarrollo, ya ha facilitado la vida a más de uno.

 

El mecanismo de funcionamiento es muy simple: se calibra la imagen que contiene el gráfico para establecer una correspondencia entre las coordenadas en píxeles de la imagen, y las unidades físicas representadas en los ejes del gráfico. A continuación, el usuario simplemente marca los puntos de la imagen que quiere medir y obtiene sus valores en las unidades físicas.

 

La calibración consiste en marcar una serie de puntos (bastan dos para la coordenada X y otros dos para la Y) para los que conocemos su valor en unidades físicas (por ejemplo, las marcas en los ejes del gráfico), e indicarle al programa cuál es este valor.

 

Si los píxeles p0 = (px0,py0) y p1 = (px1,py1) tienen valores físicos v0 = (vx0,vy0) y v1 = (vx1,vy1), entonces para cualquier otro pixel p la correspondencia será:

 

v = v0 + (p - p0)(v1-v0)/(p1-p0)

si la escala de los ejes es lineal, y:

 

v = v0 exp[(p - p0) ln(v1/v0)/(p1-p0)]

si es logarítmica.
captura de srip
Con la imagen calibrada, las mediciones se realizan simplemente haciendo click sobre la imagen, y las coordenadas físicas de los puntos marcados van apareciendo en una tabla. La incertidumbre en la medida es del equivalente a un píxel en las unidades físicas -si se tiene buen pulso al marcar- y puede reducirse haciendo zoom sobre la imagen.

 

El código, disponible aquí, está escrito en python y usa la librería Numeric para los cálculos y pyGTK para la interfaz gráfica. Requiere ImageMagik (convert) para rotar/escalar la imagen original, y puede utilizarse con prácticamente cualquier formato de imagen.


dic 02 2007

Ajustes 1-D con Python + Scipy

Tag: Herramientas y Programaciónfranwerst @ 12:30

En ocasiones es necesario ajustar una serie de puntos a una función dada, por ejemplo cuando queremos determinar parámetros físicos de una determinada línea espectral o cuando queremos modelar la transmisión de un determinado instrumento mediante un polinomio, etc...

Los ajustes unidimensionales en Python se pueden realizar de modo sencillo utilizando el módulo optimize del paquete matemático scipy. La función leastsq implementa una modificación del algoritmo Levenberg-Marquardt para minimizar la suma de los cuadrados de M ecuaciones con N incógnitas, a partir de una estimación inicial de la solución. Para obtener más información sobre el funcionamiento de leastsq con los parámetros de entrada o salida que permite se puede simplemente recurrir a la ayuda de Python:

>>> from scipy import optimize
>>> help(optimize.leastsq)

Aquí vamos a mostrar un ejemplo de su uso aplicándolo a un problema concreto como es el ajuste de una línea espectral a dos funciones: una gaussiana y a una lorentziana. El código correspondiente a este ejemplo se puede descargar desde aquí. Primero definimos las dos funciones de Python, una por cada perfil:


def gauss(params,x):
     """ Funcion gaussiana"""
     A=params[0]
     x0=params[1]
     sigma=params[2]
     g=A*exp(-0.5*power((x-x0)/sigma,2))
     return g

def lorentz(params,x):
     """ Funcion lorentziana"""
     A=params[0]
     x0=params[1]
     gamma=params[2]
     lorentz=1.0/(power(x-x0,2)+power(gamma,2))
     lorentz_norm=lorentz/max(lorentz)
     return A*lorentz_norm

Ambas fuciones dependen de tres parámetros que se introducen mediante la lista "params", y que representan la altura máxima de la curva (A), la abscisa correspondiente al máximo (x0) y una medida de la anchura (sigma en la gaussiana o gamma en la lorentziana). El vector x representa las abscisas del dominio de la función.

También necesitaremos otras dos funciones que dado el vector "y" con los puntos experimentales y un conjunto de parámetros "params" (que representan a una única curva) nos devuelva las diferencias entre los puntos y la función, es decir, calcula los residuos:

def res_gauss(params,y,x):
     """Residuos de la gaussiana"""
     return y-gauss(params,x)

def res_lorentz(params,y,x):
     """Residuos de la lorentziana"""
     return y-lorentz(params,x)

Ahora veamos cómo leastsq utiliza estas funciones para calcular el mejor ajuste. Pero primero vamos a crear una línea sintética que tenga una forma gaussiana con un cierto ruido. Definimos un vector de abscisas con 200 puntos entre x=-20 y x=20:

x=(arange(200)-100)/5.0

Creamos también un vector de ruido mediante el uso de la función random(). La intensidad del ruido será al máximo de 0.3, y para que el vector se mantenga en torno a cero le restamos 0.15 a cada punto.

random_vector=array([0.3*random()-0.15 for i in x])

Creamos la línea sintética problema con los parámetros A=-1, x0=2 y sigma=3 y le añadimos el ruido que hemos creado anteriomente.

linea=gauss([-1.0,2.0,3.0],x)+random_vector

Comencemos con los ajustes. Necesitamos unos parámetros iniciales con los que comenzar a iterar en la búsqueda de la solución óptima. Probemos con todos los parámetros iguales a 1. El algoritmo puede no converger adecuadamente si escogemos unos parámetros iniciales poco afortunados.

guess=[1.0,1.0,1.0]

Por fin realizamos los ajustes con leastsq:

ajuste_g=optimize.leastsq(res_gauss,guess,args=(linea,x))
ajuste_l=optimize.leastsq(res_lorentz,guess,args=(linea,x))

Finalmente mostramos los coeficientes resultantes de los ajustes:

print ajuste_g[0]
print ajuste_l[0]

Vemos que leastsq requiere como parámetros de entrada la función de python que calcula los residuos (esto es lo que vamos a minimizar), un conjunto inicial de parametros y después como argumentos las abscisas y las ordenadas correspondientes a la línea a ajustar. Como salida leastsq devuelve una lista cuya primera componente son los coeficientes de la solución. En nuestro ejemplo los coeficientes obtenidos son:

[-0.9927015   1.99923233 -3.01312331]   ->  para la gaussiana
[-1.08437082  2.02586623 -2.76829756]   ->  para la lorentziana

En la gráfica se pueden ver los resultados de los ajustes. En negro vemos la línea sintética que hemos creado, en rojo la gaussiana ajustada y en azul la lorentziana. Como es lógico el ajuste gaussiano es mejor ya que la línea es en realidad una gaussiana. Para verlo podemos, por ejemplo, echar un vistazo a la suma de los residuos en ambos casos:

print sum(res_gauss(ajuste_g[0],linea,x))
print sum(res_lorentz(ajuste_l[0],linea,x))
0.189616782797  -> para la gaussiana
5.68377589754   -> para la lorentziana

Algunos enlaces relacionados: