Tomar imágenes siderales con un radiotelescopio requiere de una enorme precisión, pero factores como la rotación terrestre o las holguras de sus engranajes impiden el correcto apuntamiento del cielo, desplazando y distorsionando cualquier imagen del mismo. Poder calcular estos factores y su influencia permite que radiotelescopios DIY obtengan precisiones angulares en el orden de arcosegundos, más que suficientes para radioastrónomos aficionados. Este post será dedicado a explorar mi propuesta para un software de control que tenga en cuenta todas las peculiaridades de la observación por radiofrecuencia, como parte de un proyecto más grande: la construcción de un radiotelescopio casero.

Cómo tomar imágenes del cielo con una sola antena

Una imagen, normalmente, se produce cuando señales electromagnéticas de frecuencia entre \(430\) y \(790\) terahercios (luz) son recibidas por una matriz de antenas, donde cada una mide, en paralelo a las demás, la intensidad de la señal recibida. Como cada antena individual cubre una región de origen del espacio, las medidas de intensidad para un mismo instante en el tiempo pueden ser plasmadas en una nueva matriz de intensidad de señal por región del espacio, y es dicha matriz la que conforma lo que llamamos imagen. Así funcionan las cámaras digitales y el ojo humano, pero no los radiotelescopios.

El problema con los radiotelescopios es que no podemos permitirnos el lujo de tener miles o millones de pequeñas antenas actuando en paralelo, porque las frecuencias de radiación que miden son muy bajas (del orden de los megahercios, millones de veces inferior a la de la luz). Como el tamaño de la antena escala inversamente a la frecuencia de la señal objetivo, o directamente con la amplitud de esta, lo que requiere solo unos micrómetros para luz visible, puede necesitar metros para ondas de radio.

Una solución a este problema es usar interferómetros, una unión de radiotelescopios distribuidos por todo el planeta apuntando a un mismo lugar, como se hizo para tomar una imagen del agujero negro en el centro de nuestra galaxia, Sagitario A* (figura 1). Desgraciadamente, los medios para aplicar métodos de interferometría escapan al radioastrónomo amateur, incapaz de montar una red de potentes radiotelescopios que cubra un hemisferio completo.

Pero existe otra forma de tomar imágenes del cielo con una sola antena: moverla en zig-zag cubriendo toda el área de observación, y realizar las mediciones de frecuencia e intensidad durante ese movimiento. Este no es el método ideal, porque la imagen no representará un solo instante fijo en el tiempo (cada píxel tendrá un momento temporal distinto) y la observación tomará tanto tiempo como el que tarde el radiotelescopio en cubrir el área con su movimiento. Aun así, los objetos siderales cambian tan poco en espacios cortos de tiempo, que podemos ignorar dicho cambio, siendo solo un problema el tiempo prolongado de observación.

Sin embargo, tomar imágenes con este método no es trivial. Podríamos pensar que simplemente habría que cubrir el área de observación con N puntos (digamos una cuadrícula de 16x16, por ejemplo), apuntar a cada uno de esos puntos mientras se mide la intensidad de la señal, y devolver los datos. Lo cierto es que durante la observación, los puntos que antes tenían una posición \(\mathbf p_i\) se habrán desplazado ligeramente por la rotación de la Tierra, lo que significa que, si no tenemos en cuenta este efecto, distorsionaremos la imagen. Cuantos más puntos observemos, más tiempo tardemos en apuntarlos, o más tiempo pasemos en cada uno, más se desplazarán las posiciones de los siguientes.

Para poder tomar las imágenes, propongo un sistema similar al de una impresora 3D: en vez de apuntar al primer punto (\(\mathbf p_1\)) y después ir haciendo los cálculos de desplazamiento de los siguientes puntos on-the-fly, podemos precalcular dónde acabará cada uno y reducir así la carga lógica del telescopio, dejando tres fases definidas:

1. Preparación: el usuario establece el área a escanear del cielo, su PC determina el camino a seguir y los datos son enviados al telescopio.
2. Observación: el telescopio observa el cielo de forma autónoma, siguiendo los pasos que se le han indicado previamente, devolviendo los datos al PC.
3. Procesamiento: los datos en el PC son procesados y convertidos en una imagen. El telescopio no hace nada.

Este sistema tiene dos ventajas: primero, permite tener un radiotelescopio tonto, sin mucha potencia de cómputo, ya que solo sigue órdenes predefinidas por el PC. Segundo, los cálculos no son hechos durante la observación, lo que interferiría con la misma, porque el tiempo de cálculo es tamibién tiempo durante el cual el cielo sigue rotando y distorsionando los datos, creando una situación de huevo-gallina. La única desventaja de este método es, como veremos, su complejidad matemática.

Breve nota sobre hardware y notación

Sobre la notación que voy a usar en este post quiero decir que en negrita marcaré los vectores, todos siempre de dos componentes (coordenadas esféricas). Cuando escriba \(\mathbf v_i\), me refiero al \(i\)-ésimo vector \(v\), mientras que \((\mathbf v_i)_q\) será el componente \(q\) del \(i\)-ésimo vector \(v\).

Sobre el hardware del radiotelescopio, se asume que tiene dos ejes de libertad, uno por coordenada, controlados por motores paso a paso. No se hace ninguna otra suposición.

Ahora sí, pasemos con los algoritmos.

Cálculo del camino de observación

Una vez conocemos el área de la esfera sideral que queremos observar, podemos dividirla en una cuadrícula de puntos \(\mathbf p_i = \begin{pmatrix}\theta_i, \varphi_i\end{pmatrix}^T\) y predecir sus posiciones durante la observación:

Donde \(\mathbf p'_i\) será la posición final del punto \(\mathbf p_i\) cuando vaya a ser apuntado, y \(t_j\) será el tiempo que pasó desde la observación de \(\mathbf p'_{i-1}\) hasta la de \(\mathbf p'_i\), durante el cual la esfera sideral se habrá desplazado a una velocidad constante \(\boldsymbol \omega^s\) (la velocidad de rotación de la Tierra). \(t_j\) es la suma de varios componentes:

• El tiempo que tardará en moverse de un punto al otro, dado como el máximo entre la distancia en cada eje entre la velocidad de giro del telescopio en ese eje más la holgura angular, si la hubiese, para ese eje (\((\mathbf h_i)_q\)).
• El tiempo que el telescopio pasa parado observando cada punto, \(t_o\).
• El tiempo que tarda en realizar los cálculos necesarios para la observación, \(t^c\).
• El tiempo de espera hasta comenzar la observación desde que se le indican las instrucciones, \(t^e\) (sólo para el primer punto).

El resultado es que \(t_j\) será definido de la siguiente manera:

De esta fórmula cabe destacar que \(t^e\) y \(t^o\) son valores introducidos por el usuario, mientras que todos lo demás han de ser obtenidos empíricamente.

A continuación se adjunta un simulador (simulador 1) que permite visualizar la transformación de los puntos \(\mathbf p_i\) en sus contrapartes \(\mathbf p'_i\) usando las fórmulas de esta sección:

Obteniendo los valores empíricos

Veamos ahora, uno a uno, cómo obtener los valores necesarios para resolver \(t_j\).

Holgura de los engranajes \(\mathbf h_j\)

La holgura de los engranajes sucede porque, cuando hay un cambio de sentido en la rotación del motor, sus dientes dejan de hacer contacto por un instante a causa de imprecisiones milimétricas de fabricación. Su valor es dado por el vector \(\mathbf H\) que, en cada una de sus dos componentes, almacena cuantos grados deja de girar el radiotelescopio (en ese eje) por causa de la holgura. De esta forma sabemos que, si cambiamos la dirección de uno de los ejes, debemos sumar \(\mathbf H\) a las coordenadas que queremos apuntar para deshacernos de esta imprecisión (como se refleja en la ecuación \(1\)).

Para saber si ha sucedido un cambio de sentido, introducimos dos nuevos vectores, \(\mathbf c^\rightarrow_j\) y \(\mathbf c^\leftarrow_j\), que comprueban, para cada eje, si ha cambiado. Después, unimos todo en el vector \(\mathbf h_j\), que contiene el valor correspondiente de \(\mathbf H\) o \(0\), dependiendo de si un cambio de sentido ha sucedido en el eje, permitiendo lógica condicional:

Lo único que resta es obtener el valor de \(\mathbf H\), para lo que primero definiremos el algoritmo de home, que coloca ambos ejes del telescopio en una posición sin holguras y a un paso de tocar los finales de carrera iniciales (como en una impresora 3D):

1. Nos movemos en el eje \(q\) en la dirección del final de carrera deseado.
2. Cuando hagamos contacto, nos desplazamos en sentido contrario paso a paso hasta dejar de tenerlo.
3. Repetimos con el otro eje.

En esta posición no hay holguras porque, se mueva adonde se mueva el telescopio, solo podrá ser dentro del rango delimitado por el final de carrera, en la misma dirección en la que lo dejamos.

Ahora que conocemos el algoritmo de home, pasemos con el necesario para calcular \(\mathbf H\):

1. Completamos el algoritmo de home.
2. Nos movemos hacia el final de carrera más próximo.
3. \(a^q_i\) será la cantidad de pasos que ha tardado en activarse el final de carrera. Si no hubiese holguras, sería \(1\).
4. Nos movemos en sentido contrario hasta desactivar el final de carrera, almacenamos el número de pasos en \(b^q_i\). Si no hubiese holguras, sería \(1\) también.
5. Repetimos \(n - 1\) veces los pasos 2, 3 y 4.
6. Repetimos todo el proceso con el otro eje.

La holgura angular en cada eje vendrá dada por:

Donde \(\mathbf R\) es la resolución angular (pasos por ángulo) del radiotelescopio, a la que le dedicaremos la siguinte sección.

Resolución angular \(\mathbf R\) y velocidad angular \(\boldsymbol \omega^t\)

Para obtener \(\mathbf R\) y \(\boldsymbol \omega^t\) hay dos prerrequisitos. El primero es conocer la amplitud angular de cada eje del telescopio \(\mathbf A\) (cuántos grados puede cubrir en ese eje). Esta viene dada por el diseño del propio radiotelescopio, así que nos es conocida. El segundo es calcular la cantidad de pasos que debe dar el radiotelescopio para cubrir esa amplitud, \(\mathbf S(\mathbf A)\). Obtener este último valor es tan fácil como ir a la posición de home y contar los pasos que tardamos en llegar desde ahí al final de carrera opuesto. También conviene restar a este valor un pequeño factor de corrección \(c \in [0, 1)\), porque la activación del final de carrera puede haber sucedido durante el último paso, no necesariamente al terminar.

Con \(\mathbf A\) y \(\mathbf S(\mathbf A)\) podemos calcular ya \(\boldsymbol \omega^t\) yendo a la posición de home y haciendo un barrido por eje de \(\left( \mathbf S(\mathbf A) \right)_q\) pasos, midiendo el tiempo que tardamos en hacerlo, \(\left( \mathbf T(\mathbf S(\mathbf A)) \right)_q\), y aplicando la siguiente fórmula:

Sobre \(\mathbf R\), su cálculo es directo mediante la siguiente fórmula (pasos entre grados):

Tiempo de cómputo \(t^c\)

El tiempo de cómputo es el valor empírico más fácil de obtener, por ser obtenido solo mediante software. Para ello, medimos primero el tiempo que tardamos en medir el tiempo (sí, esta frase tiene sentido), y después hacemos una observación falsa (solo los cálculos, si mover el radiotelescopio) para medir el tiempo que tardamos en hacerlo. El resultado es la resta del primer valor al segundo (tiempo bruto que se tarda en realizar la medición, sin contar el tiempo que hemos perdido por hacer la medición):

using namespace std::chrono;

microseconds::rep get_measurement_time() {
  auto b = steady_clock::now();
  auto e = steady_clock::now();
  return duration_cast<microseconds>(e - b).count();
}

microseconds::rep get_compute_time() {
  auto b = steady_clock::now();
  // Realizamos el cómputo de p'_i + t^e*w^s
  // Leemos datos de la RTL_SDR
  // Guardamos los datos
  auto e = steady_clock::now();

  auto compute_time = duration_cast<microseconds>(e - b).count()
  return compute_time - measurement_time;
}

Conclusiones

Los métodos aquí descritos permiten la toma de imágenes de objetos celestes de forma extremadamente precisa, por lo menos en la teoría, y son casi independientes del diseño del radiotelescopio sobre el que se apliquen. Aun así, quedan muchas cosas por hacer, entre ellas diseñar y construir el hardware; finiquitar el software (del que se ha dado una maqueta aquí); y hacer un algoritmo que convierta la nube de puntos en una imagen RGB. Son todas tareas bastante complejas, que ameritarán sus propios artículos en un futuro, a medida que las vaya abordando.

Sobre la implementación real de estos algoritmos, se podría pensar que realizar todos los pasos para obtener los valores empíricos (ir a home, calcular amplitud, volver a home calcular holguras, etc.) tomaría demasiado tiempo, pero muchos de ellos pueden ser realizados simultáneamente (cálculo de holguras al la vez que cálculo de resolución) y en paralelo (ambos ejes a la vez). Dependiendo de la velocidad de giro del radiotelescopio, el proceso total optimizado podría durar menos de un minuto.

Y no sé qué más escribir sobre esto. El telescopio está diseñado y escribiré pronto sobre ello, así que voy a darle una pensada a la conversión de imágenes.