El objetivo es calcular las coordenadas locales (ángulo horizontal y vertical desde el observador) de los objetos celestes a partir de las coordenadas ecuatoriales que nos ofrece Stellarium (si no lo hiciste ya, échale un vistazo a este otro post).
Coordenadas horizontales y ecuatoriales (wikipedia).
Para los cálculos nos basaremos en el método matricial desarrollado por el ingeniero japonés Toshimi Taki en Noviembre de 2002, publicado en su página web personal. En contraposición a los usuales métodos de cálculo de coordenadas celestes basados en trigonometría esférica, el método matricial ofrece una vía de cálculo alternativa, más adecuada para su implementación en microcontroladores.
En concreto, nos basaremos en las ecuaciones propuestas por Toshimi Taki para implementar a partir de ellas un sistema de cálculo de coordenadas en C++, para la placa Arduino (el código esta probado en un ATmega328). El código fuente resultante será empaquetado en una librería estándar con el objetivo de proporcionar una implementación reusable de este método para Arduino.
El método se basa en el cálculo de una matriz de transformación, para el cambio de sistema de referencia de vectores en coordenadas cartesianas, obtenidas a partir de su representación equivalente en coordenadas polares.
En el gráfico de la derecha, el punto O representa el origen de coordenadas. El vector OR muestra un vector unitario con dirección hacia el objeto celeste. La posición del objeto es expresada en coordenadas polares (ξ, ζ), donde:
- ξ representa el ángulo horizontal medido en el sentido anti-horario desde el eje X, en el plano XY (en radianes).
- ζ ángulo que forma el vector OR medido hacia arriba desde el plano XY (en radianes).
El vector OR también se puede expresar en coordenadas cartesianas (L, M, N), que serán las usadas para realizar las transformaciones.
Por un lado tendremos un vector que representará las coordenadas del objeto tomando el sistema ecuatorial como referencia, la matriz de transformación resultante permitirá obtener el vector equivalente para el sistema local de coordenadas, considerando además el factor del tiempo en las ecuaciones, y permitiendo así compensar la influencia del giro de la tierra en el segundo sistema respecto del primero. En el siguiente esquema se muestra cómo obtener el vector en coordenadas cartesianas (L, M, N) de cada sistema de referencia:
Para el cálculo de la matriz de transformación, necesitamos conocer las coordenadas en ambos sistemas de referencia de al menos tres objetos. Si bien, el tercer objeto se puede aproximar mediante el producto vectorial de los dos primeros. En el siguiente gráfico se muestran las ecuaciones implicadas en el proceso. En mayúsculas se representan los vectores de coordenadas cartesianas en el sistema de referencia ecuatorial (L, M, N), y en minúsculas los vectores en el sistema de referencia local (l, m, n). En amarillo, abajo a la izquierda, se muestra la ecuación del producto vectorial para obtener la tercera coordenada local (en minúsculas); por supuesto, la misma ecuación para las coordenadas ecuatoriales (en mayúsculas..) sería análoga.
A partir de las ecuaciones de transformación, agrupando y despejando podemos obtener la expresión de la matriz. Para hacer la transformación del vector en el sentido contrario (de coordenadas horizontales a ecuatoriales), necesitamos obtener también la matriz inversa de T.
En el recuadro de la derecha, las operaciones para obtener de nuevo las coordenadas horizontales o locales (ac = acimut, alt= altitud) y ecuatoriales (ar = ascensión recta, dec = declinación) a partir de los vectores cartesianos.
Aunque no es absolutamente necesario, para aumentar la precisión en los cálculos es recomendable elegir los objetos de referencia con una separación cercana a los 90º. Entre 60º y 120º los resultados son bastante aceptables.
La librería desarrollada para Arduino implementa las operaciones anteriores, proveyendo de algunos métodos públicos para el establecimiento de los objetos de referencia y la obtención de las coordenadas equivalentes en ambos sistemas.
Los métodos públicos reciben y devuelven coordenadas ecuatoriales u horizontales, según el caso, expresadas en radianes. El cálculo de los vectores y las operaciones con matrices se hacen de forma transparente. El siguiente fragmento de código muestra el uso básico de la librería, de forma conceptual:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
#include "CoordsLib.h" CoordsLib Coords = CoordsLib(); //Definición de la librería como variable global. //... void loop(){ //... Coords.setTime(t0); //Establecemos el tiempo inicial //Definir los dos objetos de referencia, a partir de sus coordenadas locales y ecuatoriales: Coords.setRef_1(ar_1, dec_1, t1, ac_1, alt_1); Coords.setRef_2(ar_2, dec_2, t2, ac_2, alt_2); //Obtener las coordenadas locales (ac, alt, t) a partir de las ecuatoriales y el tiempo (ar, dec): Coords.getHCoords(ar, dec, t, &ac_res, &alt_res); //Obtener las coordenadas ecuatoriales (ar, dec) a partir de las locales y el tiempo (ac, alt, t): Coords.getECoords(ac, alt, t, &ar_res, &dec_res); |
Para obtener los valores en radianes de los parámetros a partir de las coordenadas o el tiempo, podemos usar las funciones definidas en el archivo coords.py incluido en este post sobre las comunicaciones con Stellarium. La función hourStr_2_rad nos servirá para la ascensión recta (ar) y el tiempo (hora en que se toma cada medida), ambos en formato “HhMmSs“, mientras que degStr_2_rad nos servirá para la declinación (dec), el acimut (ac) y la altura (alt), todos en formato “DºM’S”“.
Por último, el código de la librería, comentado y documentado al estilo Doxygen:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 |
#ifndef CoordsLib_h #define CoordsLib_h #include <math.h> /** * \brief Biblioteca de cálculos relacionados con el cambio entre sistemas de coordenadas. * * Dispone de los métodos necesarios para el establecimiento de los tres objetos de referencia necesarios, * generación de la matriz de paso entre los sistemas de coordenadas ecuatorial y horizontal, y transformación * de vectores entre ambos sistemas. */ class CoordsLib{ private: /** * Constante de multiplicación que relaciona el tiempo solar con el sidéreo. */ float _k; /** * Timestamp con el momento inicial de las medidas. */ float _t0; /** * Indicadores de definición de los objetos de referencia. */ bool _isSetR1, _isSetR2, _isSetR3; /** * Matrices auxiliares. */ float _lmn1[3], _LMN1[3], _lmn2[3], _LMN2[3], _lmn3[3], _LMN3[3]; /** * Matriz de paso del sistema ecuatorial al horizontal. */ float _T[3][3]; /** * Matriz de paso inversa, para la transformación del sistema horizontal al ecuatorial. */ float _iT[3][3]; /** * Si se han establecido los objetos de referencia, establece la matriz de paso a partir de éstos. */ void _setT(); /** * Calcula un vector en notacion polar, a partir de las coordenadas ecuatoriales y el tiempo. * * \param ar Ascensión Recta. * \param dec Declinación. * \param t Timestamp con el momento de la medición. * \param *EVC Parámetro por referencia: Contendrá el vector de tres dimensiones en notación polar. */ void _setEVC(float ar, float dec, float t, float* EVC); /** * Calcula un vector en notacion polar, a partir de las coordenadas ecuatoriales y el tiempo. * * \param ac Acimut. * \param alt Altura. * \param t Timestamp con el momento de la medición. * \param *HVC Parámetro por referencia: Contendrá el vector de tres dimensiones en notación polar. */ void _setHVC(float ac, float alt, float* HVC); /** * Método para calcular la inversa de una matriz de 3x3. * * \param m[3][3] Matriz de entrada. * \param res[3][3] Parámetro de salida: Contendrá la inversa de m al finalizar. */ void _inv(float m[3][3], float res[3][3]); /** * Método para calcular el producto de dos matrices de 3x3. * * \param m1[3][3] Matriz de entrada. * \param m2[3][3] Matriz de entrada. * \param res[3][3] Parámetro de salida: Contendrá el resultado de la multiplicación m1xm2 al finalizar. */ void _m_prod(float m1[3][3], float m2[3][3], float res[3][3]); public: /** * Constructor de la clase. */ CoordsLib(); /** * Establece la hora de comienzo en la obtención de los objetos de referencia. * * Este parámetro es usado para considerar la influencia del transcurso del tiempo en el * sistema de coordenadas horizontal. * * \param t0 Timestamp Unix que representa el momento de inicio de las mediciones. */ void setTime(float t0); /** * Establece el primer objeto de referencia a partir de las coordenadas de un objeto * en ambos sistemas. * * \param ar Ascensión recta (sistema ecuatorial). * \param dec Declinación (sistema ecuatorial). * \param t Timestamp Unix del momento en que se toma la medida (sistema horizontal). * \param ac Acimut (coordenadas horizontales). * \param alt Altura (coordenadas horizontales). */ void setRef_1(float ar, float dec, float t, float ac, float alt); /** * Establece el segundo objeto de referencia a partir de las coordenadas de un objeto * en ambos sistemas. * * \param ar Ascensión recta (sistema ecuatorial). * \param dec Declinación (sistema ecuatorial). * \param t Timestamp Unix del momento en que se toma la medida (sistema horizontal). * \param ac Acimut (coordenadas horizontales). * \param alt Altura (coordenadas horizontales). */ void setRef_2(float ar, float dec, float t, float ac, float alt); /** * Establece el tercer objeto de referencia a partir de las coordenadas de un objeto * en ambos sistemas. * * \param ar Ascensión recta (sistema ecuatorial). * \param dec Declinación (sistema ecuatorial). * \param t Timestamp Unix del momento en que se toma la medida (sistema horizontal). * \param ac Acimut (coordenadas horizontales). * \param alt Altura (coordenadas horizontales). */ void setRef_3(float ar, float dec, float t, float ac, float alt); /** * Establece el tercer objeto de referencia de forma automática. * * A partir de los dos primeros objetos de referencia, calcula el producto vectorial de éstos * en ambos sistemas coordenadas, obteniendo así el tercer vector en cada uno de ellos. * Para que este cálculo sea efectivo, los dos primeros objetos deben distar entre sí aproximadamente * 90º. */ void autoRef_3(); /** * Calcula las coordenadas horizontales a partir de las ecuatoriales y el tiempo. * * \param ar Ascensión recta (sistema ecuatorial). * \param dec Declinación (sistema ecuatorial). * \param t Timestamp Unix del momento en que se toma la medida (sistema horizontal). * \param *ac Parámetro por referencia: Contendrá el Acimut al terminar (coordenadas horizontales). * \param *alt Parámetro por referencia: Contendrá la Altura al terminar (coordenadas horizontales). */ void getHCoords(float ar, float dec, float t, float *ac, float *alt); /** * Calcula las coordenadas ecuatoriales a partir de las horizontales y el tiempo. * * \param ac Acimut (coordenadas horizontales). * \param alt Altura (coordenadas horizontales). * \param t Timestamp Unix del momento en que se toma la medida (sistema horizontal). * \param *ar Parámetro por referencia: Contendrá la Ascensión recta al terminar (sistema ecuatorial). * \param *dec Parámetro por referencia: Contendrá la Declinación al terminar (sistema ecuatorial). */ void getECoords(float ac, float alt, float t, float *ar, float *dec); /** * Indica si se han establecido los tres puntos de referencia. * * \return Booleano. True significa que los tres puntos se han establecido. */ bool isConfigured(); }; #endif |
Para las funciones trigonométricas, incluimos el módulo math.h de la librería estándar AVR Libc (#include <math.h>). Sin embargo aún no existen utilidades para el trabajo con matrices en Arduino, por lo que se han tenido que implementar las operaciones necesarias (multiplicación de matrices, cálculo de la inversa, etc):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 |
#include "WProgram.h" #include "CoordsLib.h" CoordsLib::CoordsLib(){ _t0 = 0; _k = 1.002737908; // Constante.. Relación entre el tiempo solar (M) y el sidéreo (S): (S = M * 1.002737908) _isSetR1 = false; _isSetR2 = false; _isSetR3 = false; } /* * Calcula la inversa de una matriz m[3x3], y la devuelve en el * segundo parámetro. */ void CoordsLib::_inv(float m[3][3], float res[3][3]){ float idet; //La inversa del determinante idet = 1/( (m[0][0]*m[1][1]*m[2][2]) + (m[0][1]*m[1][2]*m[2][0]) + (m[0][2]*m[1][0]*m[2][1]) - (m[0][2]*m[1][1]*m[2][0]) - (m[0][1]*m[1][0]*m[2][2]) - (m[0][0]*m[1][2]*m[2][1]) ); res[0][0] = ((m[1][1]*m[2][2]) - (m[2][1]*m[1][2]))*idet; res[0][1] = ((m[2][1]*m[0][2]) - (m[0][1]*m[2][2]))*idet; res[0][2] = ((m[0][1]*m[1][2]) - (m[1][1]*m[0][2]))*idet; res[1][0] = ((m[1][2]*m[2][0]) - (m[2][2]*m[1][0]))*idet; res[1][1] = ((m[2][2]*m[0][0]) - (m[0][2]*m[2][0]))*idet; res[1][2] = ((m[0][2]*m[1][0]) - (m[1][2]*m[0][0]))*idet; res[2][0] = ((m[1][0]*m[2][1]) - (m[2][0]*m[1][1]))*idet; res[2][1] = ((m[2][0]*m[0][1]) - (m[0][0]*m[2][1]))*idet; res[2][2] = ((m[0][0]*m[1][1]) - (m[1][0]*m[0][1]))*idet; } /* * Multiplica dos matrices m1[3x3] y m2[3x3], y devuelve el * resultado en el tercer parámetro. */ void CoordsLib::_m_prod(float m1[3][3], float m2[3][3], float res[3][3]){ for(int i=0; i<3; i++) for(int j=0; j<3; j++){ res[i][j] = 0.0; for(int k=0; k<3; k++) //Producto de fila por columna res[i][j] += m1[i][k] * m2[k][j]; } } /* * Establece el "vector cosines" (EVC) a partir de coordenadas ecuatoriales (ar, dec, t). */ void CoordsLib::_setEVC(float ar, float dec, float t, float* EVC){ EVC[0] = cos(dec)*cos(ar - _k*(t-_t0)); EVC[1] = cos(dec)*sin(ar - _k*(t-_t0)); EVC[2] = sin(dec); } /* * Establece el "vector cosine" (HVC) a partir de coordenadas horizontales (ac, alt). */ void CoordsLib::_setHVC(float ac, float alt, float* HVC){ HVC[0] = cos(alt)*cos(ac); HVC[1] = cos(alt)*sin(ac); HVC[2] = sin(alt); } /* * Establece el tiempo de observación inicial. */ void CoordsLib::setTime(float t0){ _t0 = t0; } /* * Si ya están definidos los tres objetos, calcula T e iT. */ void CoordsLib::setRef_1(float ar, float dec, float t, float ac, float alt){ _setEVC(ar, dec, t, _LMN1); _setHVC(ac, alt, _lmn1); _isSetR1 = true; _isSetR3 = false; if(_isSetR1 && _isSetR2 && _isSetR3) _setT(); } /* * Establece el segundo de los objetos de referencia necesarios para calcular T e iT. * Si ya están definidos tres objetos, calcula T e iT. */ void CoordsLib::setRef_2(float ar, float dec, float t, float ac, float alt){ _setEVC(ar, dec, t, _LMN2); _setHVC(ac, alt, _lmn2); _isSetR2 = true; _isSetR3 = false; if(_isSetR1 && _isSetR2 && _isSetR3) _setT(); } /* * Establece el tercero de los objetos de referencia necesarios para calcular T e iT. * Si ya están definidos tres objetos, calcula T e iT. */ void CoordsLib::setRef_3(float ar, float dec, float t, float ac, float alt){ _setEVC(ar, dec, t, _LMN3); _setHVC(ac, alt, _lmn3); _isSetR3 = true; if(_isSetR1 && _isSetR2 && _isSetR3) _setT(); } /* * Calcula el tercero de los objetos de referencia a partir del producto vectorial de los dos primeros. * Si ya están definidos los tres objetos, calcula T e iT. * Si no están definidos los dos primeros, no hace nada. */ void CoordsLib::autoRef_3(){ float sqrt1, sqrt2; if(_isSetR1 && _isSetR2){ sqrt1 = (1/( sqrt( pow(( (_lmn1[1]*_lmn2[2]) - (_lmn1[2]*_lmn2[1])),2) + pow(( (_lmn1[2]*_lmn2[0]) - (_lmn1[0]*_lmn2[2])),2) + pow(( (_lmn1[0]*_lmn2[1]) - (_lmn1[1]*_lmn2[0])),2)) )); _lmn3[0] = sqrt1 * ( (_lmn1[1]*_lmn2[2]) - (_lmn1[2]*_lmn2[1]) ); _lmn3[1] = sqrt1 * ( (_lmn1[2]*_lmn2[0]) - (_lmn1[0]*_lmn2[2]) ); _lmn3[2] = sqrt1 * ( (_lmn1[0]*_lmn2[1]) - (_lmn1[1]*_lmn2[0]) ); sqrt2 = (1/( sqrt( pow(( (_LMN1[1]*_LMN2[2]) - (_LMN1[2]*_LMN2[1])),2) + pow(( (_LMN1[2]*_LMN2[0]) - (_LMN1[0]*_LMN2[2])),2) + pow(( (_LMN1[0]*_LMN2[1]) - (_LMN1[1]*_LMN2[0])),2)) )); _LMN3[0] = sqrt2 * ( (_LMN1[1]*_LMN2[2]) - (_LMN1[2]*_LMN2[1]) ); _LMN3[1] = sqrt2 * ( (_LMN1[2]*_LMN2[0]) - (_LMN1[0]*_LMN2[2]) ); _LMN3[2] = sqrt2 * ( (_LMN1[0]*_LMN2[1]) - (_LMN1[1]*_LMN2[0]) ); _isSetR3 = true; if(_isSetR1 && _isSetR2 && _isSetR3) _setT(); } } /* * Cálculo de la matriz de paso y su inversa (T e iT, respectivamente). */ void CoordsLib::_setT(){ float subT1[3][3], subT2[3][3], aux[3][3]; subT1[0][0] = _lmn1[0]; subT1[0][1] = _lmn2[0]; subT1[0][2] = _lmn3[0]; subT1[1][0] = _lmn1[1]; subT1[1][1] = _lmn2[1]; subT1[1][2] = _lmn3[1]; subT1[2][0] = _lmn1[2]; subT1[2][1] = _lmn2[2]; subT1[2][2] = _lmn3[2]; subT2[0][0] = _LMN1[0]; subT2[0][1] = _LMN2[0]; subT2[0][2] = _LMN3[0]; subT2[1][0] = _LMN1[1]; subT2[1][1] = _LMN2[1]; subT2[1][2] = _LMN3[1]; subT2[2][0] = _LMN1[2]; subT2[2][1] = _LMN2[2]; subT2[2][2] = _LMN3[2]; _inv(subT2, aux); _m_prod(subT1, aux, _T); _inv(_T, _iT); } /* * Multiplica la matriz T por el "Vector cosines" de las Coordenadas Ecuatoriales (EVC), y a partir * del "Vector cosines" de las Coordenadas Horizontales (HVC), calcula ac y alt (rac, ralt). * En caso de que no esté definido el tercer objeto de referencia, se calcula a partir del producto * vectorial entre los dos primeros antes de realizar la conversión. */ void CoordsLib::getHCoords(float ar, float dec, float t, float *ac, float *alt){ float HVC[3]; float EVC[3]; _setEVC(ar, dec, t, EVC); if(!_isSetR3) autoRef_3(); for(int i=0; i<3; i++) HVC[i] = 0.0; for(int i=0; i<3; i++) for(int j=0; j<3; j++) HVC[i] += _T[i][j] * EVC[j]; (*ac) = atan2(HVC[1], HVC[0]); (*alt) = asin(HVC[2]); } /* * Multiplica la matriz iT por el "Vector cosines" de las Coordenadas Horizontales (HVC), y a partir * del "Vector cosines" de las Coordenadas Ecuatoriales (EVC), calcula ar y dec (rar, rdec). * En caso de que no esté definido el tercer objeto de referencia, se calcula a partir del producto * vectorial entre los dos primeros antes de realizar la conversión. */ void CoordsLib::getECoords(float ac, float alt, float t, float *ar, float *dec){ float HVC[3]; float EVC[3]; _setHVC(ac, alt, HVC); if(!_isSetR3) autoRef_3(); for(int i=0; i<3; i++) EVC[i] = 0.0; for(int i=0; i<3; i++) for(int j=0; j<3; j++) EVC[i] += _iT[i][j] * HVC[j]; (*ar) = atan2(EVC[1], EVC[0]) + (_k*(t-_t0)); (*dec) = asin(EVC[2]); } /** * Indica si se han establecido los tres puntos de referencia. */ bool CoordsLib::isConfigured(){ return (_isSetR1 && _isSetR2 && _isSetR3); } |
En otros posts comentaré cómo conectar con Arduino mediante el puerto serie para pasarle los parámetros directamente desde el PC. Para las coordenadas ecuatoriales tenemos las comunicaciones con Stellarium, sin embargo para las coordenadas locales tendríamos que disponer de un mecanismo físico que podamos apuntar a una dirección concreta. El proceso de construcción de este mecanismo también lo pienso publicar más adelante, pero mientras tanto podemos usar las coordenadas acimutales (Az/Alt) mostradas en Stellarium -aunque las tendremos que copiar a mano- para calibrar y hacer las pruebas de conversión entre ambos sistemas.
Y bueno, esto es todo por el momento amigos 🙂
Salud! y saludos!