Octave y el PRN

Y ahora los detalles escabrosos. En Octave se puede generar la secuencia PRN del código C/A de los satélites GPS con la siguiente función

function [xca,g1,g2] = xcacode(s1, s2, lon)
 xca=zeros(1,lon);
 g1=zeros(1,lon);
 g2=zeros(1,lon);
 r1=ones(1,10);
 r2=ones(1,10);
 for ind = [1:lon]
  g2i=bitxor(r2(s1),r2(s2));
  g1(ind)=r1(10);
  g2(ind)=r2(10);
  xca(ind)=bitxor(g1(ind),g2i);
  new1=bitxor(r1(3),r1(10));
  new2=bitxor(r2(2),r2(3));
  new2=bitxor(new2,r2(6));
  new2=bitxor(new2,r2(8));
  new2=bitxor(new2,r2(9));
  new2=bitxor(new2,r2(10));
  r1=[new1 r1(1:9)];
  r2=[new2 r2(1:9)];
 endfor
endfunction

La secuencia de interés se obtiene en el primer elemento de retorno, xca.
r1 y r2 son los dos registros de desplazamiento de bits, con diez bits cada uno. Los parámetros de entrada s1 y s2 determinan la secuencia concreta que se devuelve, indicando lon la longitud. Para los satélites 05, 19, y 25 sus PRN se obtienen ( Interface Specification IS-GPS-200, pg. 7 y 8) así

xca05=xcacode(1,9,1023);
xca19=xcacode(3,6,1023);
xca25=xcacode(5,7,1023);
...
xca19_ = de01a1m1(xca19);
...

Para pasar de la versión "lógica" de 0 y 1 a la versión de señal -1 y +1 se usa esta función bien corta (gracias a la sugerencia de Misra&Enge):

function v1m1 = de01a1m1(v01)
 v1m1 = (v01==0)*(1) + (v01==1)*(-1);
endfunction

Las siguientes funciones desplazan (shift) circularmente (circ) a la izquierda o a la derecha una secuencia ini un número de pasos n

function  res = shiftcircizq(ini,n)
 siz = max(size(ini));
 if (abs(n)>= siz )
  res=ini;
 elseif (n>=0)
  res = ini([n+1:siz 1:n]);
 else
  n=-n;
  res = ini([siz-n+1:siz 1:siz-n]);
 endif
endfunction

function  res = shiftcircder(ini,n)
 res = shiftcircizq(ini,-n);
endfunction

Por ejemplo

> shiftcircizq([0 0 1 0 1 1 1 0],2)
ans =
1   0   1   1   1   0   0   0
> shiftcircder([0 0 1 0 1 1 1 0],2)
ans =
1   0   0   0   1   0   1   1

La función de correlación circular de dos vectores x e y se calcula con la siguiente función

function res = correlacir(x,y)
 siz=max(size(x));
 res=zeros(1,siz);
 for i = [1:siz]
  res(i) = sum(x.*shiftcircizq(y,i-1));
 endfor
endfunction

que corresponde a esta definición matemática

La primera línea es la definición, y la siguiente describe propiedades de simetría al intercambiar el orden de las funciones a correlacionar. Según esta definición para la correlación con n=10 de x con y, multiplico término a término los elementos de x con una copia de y desplazada a la izquierda 10 posiciones. x(1) se multiplica con y(1+10), etc.

No es por tanto lo mismo la correlación de x con y que la de y con x, para un mismo valor de n. Pero sí es lo mismo la correlación de x con y para n=10, que la de y con x para n=-10.

Hay que andarse con cuidado a la hora de poner en el orden adecuado los vectores en la llamada a la función de correlación. En la anterior entrada del blog, la gráfica de la correlación con el pico en el paso 350 se hacía mostrando

cc19compruido=correlacir(xca19_, xcomp+ruido);

xca19_representa al PRN19 generado en el receptor(con valores +1 y -1), xcomp+ruido es la suma de las señales de los tres satélites y el ruido aleatorio, la señal total recibida. A la vista del orden, lo que se hace es coger el vector PRN19, luego el xcomp+ruido al que se desplaza a la izquierda entre 0 y 1022 pasos, obteniendo la función de autocorrelación para n=0 hasta n=1022. Pero resulta que eso es lo mismo que dejar xcomp+ruido quieto, y desplazar a la derecha xca19_.

Sin embargo si se calcula la correlación con correlacir(xcomp+ruido, xca19_), el pico del máximo no sale en el paso 350, sino en el 673 (es decir, 1023-350). ESto significa lo mismo pues por la periodicidad en 1023 de la función de correlación el máximo en 673 también está en -350, y desplazar a la izquierda -350 pasos el segundo parámetro, xca19_ , es justamente desplazarlo +350 a la derecha. Es lo mismo, pero si se interpreta bien.

En cuanto al ruido, se "genera" usando una función de octave ruido=4*randn(1,1023);, y hay que recordar que cada vez que se invoca la función se obtienen diferentes valores.

Para hacer las gráficas de un vector v se puede usar en octave el plot(v) sin más, y se pueden usar símbolos discretos para cada valor, como se hizo en la figura de autocorrelación de las PRN 05, 19 y 25. Si en vez de símbolos se usa una línea tendremos que esta consiste en segmentos rectilineos que unen los diferentes valores. Eso no sirve directamente para sacar un vector de valores +1 o -1 como pulsos rectangulares. Hay un tipo de gráfica en octave, el bar(...) que tampoco es lo que necesitaba, de modo que lo mejor es construir a partir de un vector de valores, otro con los puntos adecuados que unidos den los pulsos rectangulares deseados, una función como

#{
Una secuencia de valores en un vector v se quiere representar en forma de pulsos rectangulares,
que se inicien en x0 y de anchura xdelta. Por cada valor de v se emplean dos puntos desde el inicio al final
del intervalo correspondiente, a la altura del valor v. El primer punto se asume que es (x0,0)
#}
function [x,y] = graficabin(v, x0, xdelta)
 n=max(size(v));
 lon=2*n+1;
 x=zeros(1,lon);
 y=zeros(1,lon);
 x(1)=x0;
 y(1)=0;
 for i = [2:lon]
  y(i)=v(floor(i/2));
  x(i)=x0+xdelta*floor((i-1)/2);
 endfor 
endfunction

function plotbin(v,x0, xdelta)
 [x,y]=graficabin(v, x0, xdelta);
 margenx=(max(x)-min(x))/100;
 margeny=(max(y)-min(y))/50;
 plot(x,y);set (gca (), "xlim", [min(x)-margenx, max(x)+margenx]);set (gca (), "ylim", [min(y)-margeny, max(y)+margeny]);
endfunction

La segunda función es útil como sustituto rápido del plot.

Octave sigue siendo una herramienta inmejorable en relación calidad/precio para cacharrear con números, y hacer deberes. Aunque a veces hasta que consigues hacer algo pasas un buen (o mal) rato, como la simple tarea de comprobar los 10 primeros dígitos de las secuencias PRN del código C/A en octal, o los 16 primeros en hexadecimal. Pero al final es algo tan sencillo como
convertir de número a texto, quitar los espacios, pasar de binario a decimal y de decimal a octal o hexadecimal

> xca19(1:16)
ans =
1 1 1 0 0 1 1 0 1 1 0 1 0 1 1 0
> dec2base(bin2dec(strrep(num2str(xca19(1:10))," ","")),8)
ans = 1633
> dec2hex(bin2dec(strrep(num2str(xca19(1:16))," ","")))
ans = E6D6


La magia del PRN

La magia del PRN consigue que un cacharro "gepeese" sea capaz de escuchar a los satélites que allá arriba le mandan sus señales.

Un satélite GPS gira en una órbita casi circular a unos 20.180 km de altura sobre la superficie terrestre. En enviar su señal no consume demasiado, unos 27 vatios, como una bombillita de poca luz. Emplea una buena antena, que enfocada hacia la Tierra distribuye algo más de energía hacía las zonas más alejadas, respecto de las que le quedan justo debajo. Esto compensa en parte el hecho de que el debajo "sólo" está a 20.180km, y la parte rasante a unos 25.700km. Con todo, dada la lejanía del satélite respecto de un receptor típico, a este le llegan aproximadamente 10^-16 vatios. La señal es tan rematadamente débil que está inmersa en medio del ruido radioeléctrico existente debido a todo tipo de fuentes humanas presentes en la propia superficie.

Es como estar al fondo de una larga sala repleta de gente y esperar que nuestras palabras sean entendidas por alguién en el otro extremo, rodeado de personas en animada conversación.

Obviando muchos detalles, el caso es que podemos considerar que además de otra información, un satélite GPS repite sin cesar, cada milisegundo, una secuencia binaria (de ceros y unos) de 1023 dígitos. Cada satélite emplea una secuencia única, diferente a la empleada por el resto. El orden de los ceros y unos en la secuencia parece a primera vista aleatorio, sin una pauta clara, por lo que se les conoce como secuencias PRN.

PRN son las siglas inglesas de Pseudo Random Noise, ruido pseudo aleatorio, aunque puede usarse Number en vez Noise, entendiendo que no es pseudo aleatorio el número sino la secuencia. La señal formada por la repetición de esas secuencias de 1023 dígitos es lo que se llama el código C/A, pensado para uso civil, y encargado de ayudar a que el receptor pueda Adquirir (la A de C/A) la señal del satélite, además de hacer con él una estimación burda (la C en C/A es de Coarse) de la pseudodistancia. Hay otros códigos con secuencias PRN mucho más largas, como el P (de Preciso, en contraste a burdo) pensado para uso militar, o su versión encriptada Y (encrYpted).

El "Pseudo" indica que realmente la secuencia de ceros y unos no es al azar sino que se determina exactamente mediante ciertas reglas. Sin embargo tiene unas propiedades muy similares a las secuencias genuinamente aleatorias, y eso es lo que importa para su aplicación. Claro que el código C/A nos puede dar el pego de aleatorio como mucho 1 milisegundo, pues luego la secuencia se repite machaconamente.

¿Cuales son las reglas de creación de cada secuencia única? Se pueden explicar en términos de registros de desplazamiento de bits, y de operaciones lógicas. El caso es que esas reglas que usa el satélite, son conocidas y empleadas por los receptores. Y por cualquiera interesado si se lee por ejemplo la Interface Specification IS-GPS-200.

Yo he preferido seguir un excelente libro de texto sobre el GPS, el "Global Positioning System. Signals, Measurements, and Performance", segunda edición, de Pratap Misra y Per Enge, y hacer los "deberes" del final del capítulo 2, "GPS in 2005: An overview". Allí se pide calcular las secuencias PRN de 3 satélites distintos, el 05, el 19 y el 25, según las reglas explicadas. Las reglas son comunes a todos los satélites, cada satélite se distingue por usar unos parámetros con valores específicos. Una etiqueta como por ejemplo PRN05 designa tanto una secuencia PRN concreta, como el único satélite GPS que la emite.

Usando de nuevo octave, las gráficas de las secuencias de esos 3 satélites tienen esta pinta:



La secuencia de 1023 ceros y unos lógicos se muestra como una señal digital de pulsos, usando un valor de la señal de +1 para representar los ceros, y una señal con valor -1 para representar los unos. Cada dígito se convierte así en lo que se llama un "chip", un pulso rectangular de amplitud positiva o negativa. Hay 1023 chips en cada secuencia repetida. Visto como en la imagen anterior parece que se trata de un código de barras, pero acercandonos para ver sólo los 200 primeros chips de cada satélite, se aprecia mejor el sube y baja de la señal, sobre todo si se abre la imagen en una ventana para ella sola.



En estas figuras en el eje horizontal se representa el tiempo de emisión del chip. Si el primer chip se emite en t=0, el segundo se emite en t=1, ..., y el último chip de la secuencia, el 1023, se emite en t=1022, todo expresado en unidades de chip (de duración de chip), donde 1 chip equivale a 1/1023 milisegundos.

Y ahora llega la propiedad mágica crucial de las secuencias PRN. La correlación entre dos secuencias PRN de dos satélites cualesquiera es casi nula. La autocorrelación de cada secuencia PRN consigo misma desplazada es también casi nula, salvo para el desplazamiento nulo.

La correlación entre secuencias cualesquiera de ceros y unos, de la misma longitud, es algo muy sencillo. Poniendolas en paralelo se van comparando los elementos que ocupan la misma posición. Si sus valores coinciden, sumo 1. Si difieren, resto 1. Lo que me dé la suma al final es la correlación, si bien puedo normalizar diviendo entre la longitud de las secuencias. Es como hacer una operación de lógica booleana, el OR exclusivo o XOR.

Otra forma de cálcular el mismo resultado, empleando los niveles del "chip", +1 para los ceros y -1 para los unos, es la multiplicación. Cuando hay coincidencia tanto 1x1 como (-1)x(-1) dan 1. Cuando hay discrepancia (-1)x1 = 1x(-1)= -1. Agregando los resultados de tales multiplicaciones llegamos a lo mismo que antes, aunque así tenemos una definición más general, aplicable a cualquier señal con amplitudes tanto negativas como positivas, lo que nos será útil más adelante.

Si comparamos una secuencia con una copia de sí misma desplazada un número variable de pasos, se habla de autocorrelación. Si la copia no se desplaza (o se "desplaza" 0 pasos), siempre se suma 1 por cada dígito o chip, nunca se resta. Para una secuencia genuinamente aleatoria de 1023 dígitos la autocorrelación (sin normalizar) será de 1023. Por otro lado si la comparo con una copia exacta pero desplazada a derecha o izquierda, por la misma aleatoriedad de la situación de ceros y unos, casi tantas veces tendré coincidencias, sumando 1, como discrepancias, restando 1, y la suma total estará muy cerca de 0.

Este tipo de correlación es "circular". Al comparar una secuencia con otra desplazada a la derecha 10 pasos (cada paso corresponde a un dígito o a un chip), entendemos que el valor de la posición 1 pasa a la posición 11, el de la 2 a la 12, ..., el de la 1013 a la 1023, el de la 1014 a la 1, ..., y el de la 1023 a la 10. Si el desplazamiento es de 10 pasos a la izquierda, la posición 1 pasa a la 1014, ..., la 10 a la 1023, la 11 a la 1, ..., y la 1023 a la 1013.

La siguiente imagen refleja los valores de la autocorrelación para las PRN 05, 19 y 25, cada una con una copia desplazada un número de pasos (eje horizontal) entre 0 y 50. El aspecto para los pasos del 51 al 1022 es parecido al de los pasos 1 a 50. Un desplazamiento de 1023 pasos lleva a la configuración inicial.



Por cierto, esos pasos ¿son hacia la derecha o hacia la izquierda? En el caso de la autocorrelación da igual, por la simetría. Si comparo A con una copia C, de A desplazada a la izquierda 10 pasos, es igual que comparar C con una copia suya desplazada a la derecha 10 pasos, que resulta ser A.

Hay que fijarse que para desplazamiento cero siempre se obtiene la máxima correlación de 1023, faltaría más. Las funciones de autocorrelación en cada uno de los tres casos no coinciden exactamente, pero todas (esas tres y las del resto de satélites) comparten algo en común: los únicos valores de autocorrelación (sin normalizar) que se dan para cualquier paso entre 1 y 1022 son los tres siguientes: 63, -1 y -65. La clave es que esos valores son mucho más pequeños que el máximo, 1023.

Además la correlación entre la PRN de un satélite y la de cualquier otro satélite distinto, para cualquier paso (incluido el cero) también tiene siempre uno de esos tres valores: 63, -1, -65.

Las PRN de los satélites GPS para el código C/A son un tipo de secuencias matemáticas estudiadas por el Dr. Robert Gold, por lo que se conocen como códigos de Gold. Las propiedades enunciadas arriba se pueden probar rigurosamente.

El caso es que cada satélite emite una especie de huella dactilar. Y los receptores GPS tienen un "detector de huellas", un correlacionador. Es un circuito electrónico que genera la secuencia PRN de chips de un satélite dado, o su versión desplazada cierto número de pasos, y durante un milisegundo calcula la correlación de dicha secuencia con la señal que le llega.

Imaginemos que el correlacionador va a ver si recibe algo del satélite 19 en la señal que, captada por la antena del receptor, es acondicionada previamente por otros componentes del receptor. Para ello calculará la correlación entre la PRN19 generada localmente en el receptor y la señal. Por cada valor de la función de correlación para un cierto número de pasos de desplazamiento, hay que generar durante un milisegundo la secuencia PRN19 desplazada ese número de pasos, multiplicarla chip a chip por la señal recibida e ir sumando para tener al final el valor buscado. Este cálculo hay que repetirlo para los pasos del cero al 1022, hasta obtener los 1023 puntos de la función correlación.

Imaginemos también (problema 2.2(f) de Misra&Enge) que la señal que le llega es una suma de una señal completamente aleatoria, con distribución normal, media cero y desviación estándar de 4 unidades, junto con tres señales de los satélites 05, 19 y 25, de amplitud 1 (inferior pues a la del ruido aleatorio).

Como la distancia de cada satélite hasta el receptor será diferente, el origen de la secuencia de cada satélite (emitido cada milisegundo del reloj del satélite), llegará al receptor con un retraso variable. Supongamos que respecto a la marca de milisegundo del reloj del receptor, cuando comienza este a generar su PRN particular, el origen de la PRN05 recibida del satélite va retrasada 75 chips, el origen de la PRN19 recibida va retrasada 350 chips, y el origen de la PRN25 recibida va retrasada 905 chips. (Aunque es equivalente decir, por ejemplo, que el origen de la PRN05 desde su satélite va adelantada 1023-75 = 948 chips).

Bueno, pues la mezcla de las señales indicadas tendrá el siguiente aspecto, que dificilmente da idea de satélite alguno.




Teniendo que calcular 1023 valores, empleando un milisegundo por cada uno, en poco más de un segundo el correlacionador obtiene ... ¡algo así!:



Tatacháaaaaaan. Un maravilloso pico que destaca primoroso ¡justo en el paso 350! ¡Magia!

Ese pico por un lado descubre el código PRN del satélite 19 entre el aparente ruido informe, y por otro no menos importante, mide la pseudo distancia entre el satélite 19 y el receptor, salvo por un número entero de milisegundos.

Con la información del correlacionador, el receptor no sabe si la señal que le llega salió del satélite 19 hace 66 milisegundos y pico, o hace 67 milisegundos y pico, o hace 68 milisegundos y pico, o... pero sí sabe que el "pico" es de 350/1023 milisegundos. Cada chip dura 1/1023 milisegundos, tiempo en que la luz recorre "sólo" unos 293 metros. Los circuitos electrónicos pueden afinar y precisar el paso fraccionario en que se da el máximo de correlación, llegando a las milésimas, por ejemplo 350,276 en vez de 350. En una milésima de chip la luz recorre casi 30 centímetros, y ese es el orden de error en la estimación de la pseudo distancia basada en el código C/A (sin olvidar que es pseudo, no la distancia real, por el tema de la falta de sincronía entre el reloj del receptor y del satélite, y que hay otras fuentes de errores).

Vamos a ver esto del "pico". La función de correlación de la figura se obtiene al comparar la señal recibida con una copia de la secuencia PRN19 desplazada n pasos a la derecha. En términos temporales "desplazar a la derecha" significa "emitir con retraso": cada valor en la copia se emite n chips más tarde que en el original. El máximo encontrado en el paso 350 se da al utilizar una copia de PRN19 retrasada 350 chips. Eso indica que hay una coincidencia con el PRN19 original del satélite cuando esperamos 350 chips tras la marca de milisegundo del receptor. Esos 350 chips son el intervalo de espera entre una marca de milisegundo del reloj del receptor y la siguiente marca de milisegundo del reloj del satélite.

Tras captar el máximo de correlación y adquirir la señal del satélite 19, consideremos el momento en que por el reloj del receptor llega la marca de milisegundo X. El receptor sabe que llegará una marca de milisegundo del reloj del satélite dentro de 350/1023 milisegundos, es decir, cuando sea el tiempo de recepción X+(350/1023)ms. Pero no sabe qué marca de milisegundo es la que llega, cuándo fué emitida esa marca por el satélite. Necesita aún que alguién le "sople" que esa marca del satélite es la del milisegundo Y. Entonces el receptor ya podrá estimar la pseudodistancia como c por el intervalo temporal entre emisión y recepción: 0,001.c.(X+(350/1023)-Y).

¿Quién le sopla al receptor el tiempo exacto de emisión, según el reloj del satélite, de la marca de milisegundo que le llega? Pues la propia señal del satélite, que además del código C/A con su PRN, incluye bits de datos con todo lo necesario. Aunque como se transmiten a un ritmo muy lento de 50 bits por segundo, hace falta un poco de paciencia para poder seguir la señal (tracking). Pero antes de poder seguirla, y de poder acceder a esos bits de datos, hay que adquirir (adquisition) la señal, que es donde interviene el correlacionador.

Pero nunca es tan facil. Hay tantos detalles... Hemos asumido que la antena capta la señal, algún sistema electrónico la "acondiciona", y ya puede el correlacionador calcular punto a punto la correlación. Ya es bastante trabajo que esa función tenga 1023 puntos, y que haya que probar las PRN de todos los satelites, porque si el receptor de entrada no sabe dónde está, no sabe qué satélites pueden estar a la vista. Es que además la señal que el satélite emite en una frecuencia muy determinada se modifica por el efecto Doppler, debido al movimiento del satélite y el posible del receptor. Este tiene que buscar en varios intervalos de frecuencias, sin saber exactamente en cúal estará la señal de cada satélite, y por cada intervalo elegido hay que repetir los cálculos del correlacionador.

Es como si éste fuese un labriego que busca un anillo perdido en un surco, recorriendo con cuidado el mismo a ver si descubre su brillo, y recibe una aclaración de su amada, y propietaria del anillo: "... pero no estoy segura de que fuese en ese surco, salté por todo el terreno, y perdí no uno sino todos mis anillos, 4 ... o 10, no estoy segura".

Menos mal que estos cacharricos electrónicos son muy rápidos, y no tienen venas.

La magia del GPS

Parece que Arthur C. Clarke, además de escribir libros de Ciencia Ficción y de idear el satélite de comunicaciones geoestacionario, dijo eso de que

Cualquier tecnología suficientemente avanzada es indistinguible de la magia.

A los efectos prácticos de cualquier usuario de navegador parlanchín que nos manda girar a izquierda o derecha, que ese cacharro sepa dónde estamos y por dónde ir es asimilable a la magia.

¿Qué "magia" usa el navegador o el móvil o cualquier dispositivo con GPS para saber dónde está?

Aunque llamamos a menudo "el gepeese" a los cacharritos que nos indican dónde estamos, realmente GPS hace referencia al Sistema de Posicionamiento Global, propiedad del ejército USA. Los soviéticos dieron la réplica con su GLONASS, y la Unión Europea anda tras una versión propia y civil, Galileo, al que citaba en noticias del espacio. También China está poniendo en marcha su propio sistema Compass.

Todos esos sistemas se engloban bajo unas siglas, GNSS, Sistemas de Navegación Global por Satélite. Permiten determinar (para ayudar en la Navegación) la posición y velocidad de un receptor, en cualquier parte de la Tierra (por eso lo de Global), con la ayuda de una constelación de Satélites. Para conseguirlo se hace uso no de una, sino de muchas, muchas magias acumuladas por los hechiceros científicos a lo largo de los siglos.

La magia más esencial y poderosa, que sustenta al resto de magias, es la Matemática. Por ejemplo, la geometría de un espacio euclídeo tridimensional nos permite resolver el problema siguiente: Dados tres puntos distintos S₁, S₂ y S₃ no colineales, determinar el punto R tal que su distancia a S₁ es D₁, a S₂ es D₂ y a S₃ es D₃. Basta considerar las esferas centradas en cada punto S_i, de radios dados por las distancias D_i, y en su intersección común, de existir, se encontrará R. Por la simetría de las 3 esferas consideradas respecto al plano determinado por sus 3 centros, si hay una solución a un lado del plano, habrá otra "reflejada", y si la solución resulta estar en ese plano, será única.

Nuestros modernos cacharros electrónicos no saben hacer construcciones geométricas, pero se les da de maravilla operar con números. Un gran paso en la historia de las Matemáticas fue el desarrollo de la geometría analítica. Las construcciones geométricas se convierten en ecuaciones algebraicas, y la solución concreta de encontrar las 3 coordenadas cartesianas de R dadas las 9 de S₁, S₂ y S₃, y las tres distancias D₁, D₂, D₃, implica resolver 3 ecuaciones (de segundo grado) con 3 incógnitas.

(R_x-S_ix)²+(R_y-S_iy)²+(R_z-S_iz)² = D_i²   (i=1,2,3)

Ya podemos usar los rápidos calculistas electrónicos, cada vez más chicos y poderosos gracias a las magias cuánticas y lógicas, para resolver las ecuaciones.

Pero antes de resolver las ecuaciones y despejar las incógnitas, hay que obtener los datos, ¡y esa es la parte más difícil! Esos 3 puntos S₁, S₂ y S₃ no colineales van a ser 3 satélites, no 3 puntos fijos. La razón para usar una constelación de satélites esta en la G de GPS, que allá donde vayamos sobre la Tierra y sus inmediaciones aéreas, encontremos los puntos de referencia necesarios, sobre nuestras cabezas. Resulta que con unos 24 a 30 satélites bien repartidos hacemos el apaño.

Mediante las magias correspondientes, largas de enumerar, se construyen y lanzan los satélites, y se les "sigue" estrechamente para saber no sólo dónde están sino por dónde van a estar, decírselo a los propios satélites, y que estos a su vez se lo digan a todo receptor que quiera y pueda entenderlos. Un satélite de un sistema GNSS en esencia va emitiendo continuamente una información como "soy satelifulanito de tal, ahora que son las tantas estoy justo en tales coordenadas".

Bueno, si soy el punto R del problema geométrico y a la vez un Receptor que entiende las señales del satélite S₁, del S₂ y del S₃, podría en un instante dado tener los datos de las 9 coordenadas de los 3 puntos que ocupaban los satélites cuando salieron los 3 mensajes que me llegan a la vez. Solo me hace falta obtener la distancia a cada uno de esos 3 puntos.

Aquí la geometría recibe una ayuda de la física, que le dice "esa información desde el satélite llega al receptor por una onda electromagnética que viaja a la velocidad de la luz, c=299.792.458 m/s". Como el espacio recorrido es igual a la velocidad por el intervalo de tiempo empleado, y en el mensaje desde el satélite me dan "la hora" en que estaba en "esas coordenadas", miro el momento en que me ha llegado el mensaje (más tarde que cuando salió del satélite), hago la diferencia entre tiempo de salida y tiempo de llegada, multiplico por c, y ya tengo (repitiendo el proceso para cada satélite) D₁, D₂, D₃, y a resolver las ecuaciones.

Llegados a este punto nos damos de bruces con un problemilla. Dar con la distancia entre S₁ y R implica restar tiempos, medidos por diferentes relojes. S₁ y R cada uno tiene su reloj. Si ambos estuviesen marcando exactamente la hora en perfecta sincronía, nada que objetar. Los relojes que van en los satélites son de lo mejorcito, relojes atómicos de Cesio o de Rubidio, varios en cada satélite por redundancia. Pero el reloj de un receptor, sin ánimo de ofender, no les llega a las suelas de los zapatos.

¿Entonces? Imaginemos un caso aproximadamente real. El satélite S₁ está sobre la cabeza de nuestro receptor R a 20.100km de altura, y pasa por allí a las 12h 0minutos 0segundos 0 microsegundos mandando su mensaje, que le llega a R cuando por el reloj de R han pasado 66 milisegundos de las 12h en punto. Redondeando c=300.000km/s, la distancia estimada por R será 0,066.c = 19.800km. Por llevar el reloj adelantado 1 milisegundo (marca 66 cuando son 67), R comete un error en su estimación de la distancia a S₁ de 300km. Incluso un adelanto tan pequeño como 10 microsegundos en el reloj de R implicaría un error en D₁ de 3 km.

La conclusión es que R no puede estimar directamente la distancia verdadera a un satélite, D, pero sí una pseudodistancia pseudoD, que encierra una discrepancia c. deltaR respecto a la verdadera por culpa del adelanto (deltaR > 0) o atraso (deltaR < 0) de su reloj.

La nueva incógnita deltaR es desconocida pero la misma para cualquier distancia a diferentes satélites. Una nueva incógnita requiere, para resolver el sistema de ecuaciones, una nueva ecuación, que se obtiene fácilmente si además de los satélites S₁, S₂ y S₃ interviene un cuarto satélite S₄ favorablemente situado. El nuevo sistema de ecuaciones a resolver es

(R_x-S_ix)²+(R_y-S_iy)²+(R_z-S_iz)² = (c . deltaR)² + pseudoD_i²   (i=1,2,3,4)

En el ejemplo de antes, uno de los datos obtenidos por R sería que pseudoD₁=19.800km. Al resolver el sistema con las otras 3 pseudodistancias y las posiciones de los 4 satélites, debería obtenerse la posición de R y que deltaR es aproximadamente 1 milisegundo de adelanto sobre el tiempo verdadero.

En esencia así se obtiene la posición de R. Aunque nunca es tan fácil. En la vida real todo es fuente de error. Las posiciones de los satélites no se conocen con total exactitud. Las ondas electromagnéticas que llevan información del satélite al receptor pasan por la ionosfera y la estratosfera, donde su velocidad varía según las condiciones ambientales. Las características de esa señal se ven influidas además por el movimiento tanto del satélite como del receptor, con el efecto Doppler cambiando frecuencias. También hay que contar con los efectos relativistas. Y con el ruido y demoras en los circuitos electrónicos. Incluso los avanzados relojes atómicos de los satélites tienen sus errores.

Ante tanta "fuente de error", la Física y la Matemática contraatacan modelizando esas fuentes y sus "errores". En general se sigue una estrategia como la anterior cuando se ampliaban las ecuaciones de 3 a 4, con una nueva incógnita. Cada fuente de variabilidad adecuadamente modelizada conlleva una serie de nuevos parámetros con valores desconocidos pero que se pueden estimar incorporando nuevas ecuaciones. Cuanto más datos y más ecuaciones, más se mejoran las estimaciones de los parámetros de los modelos, y más se reducen los errores.

La resolución de ecuaciones y estimación de parámetros es un asunto denso y un tanto aburrido, no toda la magia es igual de atrayente. Pero hay números de magia especialmente sorprendentes, como el que hace cada Satelifulanito para decir quién es y, el más difícil todavía, con el que un Receptor avispado le escucha e identifica incluso en medio de un ensordecedor ruido. Es ... la magia del PRN.

noticias del espacio

El año 2010 ha empezado con bastantes novedades relacionadas con la exploración espacial.

Por el lado americano la propuesta de Obama para el presupuesto de la NASA ha cambiado radicalmente los planes de la era Bush de volver a la Luna en 2020, cancelando el programa Constellation, con sus lanzadores Ares I y V y el vehículo Orion.

Cancelar un proyecto con sobrecoste, retrasado y con unas miras prácticas escasas, salvo el gustirrinín que debió tener pasajeramente George al anunciarlo en su día, parece una gran noticia. Pero no es que se le haya ocurrido enteramente a Obama. Parece determinante el "informe Augustine", obra de un comite para estudiar las planes de exploración espacial humana por parte de U.S.A., presidido por Norman R. Augustine. Los miembros del comite tienen todos un curriculum excelente y parece que saben algo del asunto. En las observaciones finales, en el punto 9.2 sobre "cuadrar los recursos con las metas" tenemos

Quizá la mayor contribución al riesgo en el programa espacial, tanto humano como financiero, es el pretender alcanzar logros extraordinariamente difíciles mediante recursos inconsistentes con las demandas de tales tareas.
...
En el programa Constellation, el coste estimado del desarrollo del lanzador Ares I creció conforme la NASA determinó que el plan original de usar los cohetes principales de la Lanzadera Espacial en la fase superior del Ares I sería demasiado costoso, en parte por tener que añadir un autoarranque. Pero el motor sustituto tenía menos empuje y menor economía de carburante, con lo que los cohetes sólidos de la primera fase tenían que modificarse para proporcionar mayor empuje. Esto a su vez contribuía a un fenómeno vibratorio, cuya corrección aún está por demostrar plenamente. Esta es la naturaleza de los programas de desarrollo complejo, con presupuestos que más probablemente se recortan en vez de aumentar.
...
Los presupuestos menguantes y unas reservas inadecuadas, no sólo en dolares sino también en tiempo y tecnología, son una fórmula para un casi seguro fracaso en el vuelo espacial tripulado. Si no hay disponibles recursos acordes con las metas establecidas, deben adoptarse nuevas metas.

Sobre la decisión de Obama me ha llamado la atención la opinión de David Mindell, autor de un libro cuya lectura me dió para dos incursiones lunáticas, Luna 1 y Luna 2. En su libro "Digital Apollo" se reflejaba la tensión entre el lado humano, con los astronautas en primer plano, y el tecnológico, con la automatización y el control. La primera visita a la Luna tenía un trasfondo de guerra fría y mediático, pero con lo que han dado de sí las misiones robóticas en Marte, ahora hay mucho a ganar ahorrando el factor humano en la ecuación de explorar el entorno más allá de la baja órbita.

Y en la baja órbita hay mucho que aprender sobre la vida humana en ingravidez en largos periodos, como los que requeriría una estancia lunar permanente o un viaje a Marte. La gran beneficiada de los nuevos planes va a ser la Estación Espacial Internacional, cuya vida se prolongará hasta 2020. De acuerdo al administrador de la NASA y ex-piloto de transbordador Charlie Bolden

Empezaremos a usar la Estación Espacial Internacional (ISS) como el laboratorio nacional que se imaginó. Haremos uso pleno de su increible potencial, mejorando nuestro uso de sus capacidades de investigación y desarrollo a bordo.
...
Hay mucho por conocer antes de aventurarse con seguridad lejos de la baja órbita terrestre de forma prolongada.
...
También instalaremos instrumentos de observación en la ISS para ampliar el conocimiento de nuestro planeta y usar esta plataforma para probar futuras tecnologías de exploración.

Precisamente ayer despegaba la misión STS-130 de la lanzadera espacial, la quinta última. Y en ella se plasma la colaboración internacional, en este caso de USA y la Unión Europea, o la NASA y la ESA, pues el Endeavour lleva en sus bodegas el Nodo 3, Tranquility, y la espectacular Cupola que se aprecia en el logo de la misión. Estas dos piezas se han construido en Europa principalmente por Thales Alenia Italia ( de ahí el italiano Cupola, supongo que por Cúpula), y como se indica en la nota de prensa esto se enmarca en el acuerdo de intercambio de la ESA con la NASA por lanzar ésta el Columbus. Bueno, la industria europea puede presumir. Es interesante el anterior documento de la misión STS-130 pues hay fotografías de los citados módulos. En una foto se ve la estructura de la Cupola en octubre de 2002, tras ser forjada en una sola pieza en Francia. El acuerdo de desarrollo del nodo 3 fué en 1997. Se ve que esto requiere mucho tiempo, ¡como para no extender al máximo el tiempo de vida de la ISS!.

Y hablando de espacio y Europa, puede correr el cava y el champan, y tañer las campanas, y festejar cual carnaval por las calles que ... ¡Galileo vive!

No es que el magnífico científico italiano se levante de entre los muertos, pero es casi tan increible, tras años de demoras, más que nada de tipo politico-económico, el sistema europeo de navegación por satélite parece despegar definitivamente. A primeros de enero se anunció, y se han firmado el 27 de enero tres contratos, uno para la construcción de 14 satélites a OHB (segmento espacial), otro para el soporte del sistema industrial a Thales Alenia, y otro a Arianespace para servicios de lanzamiento. En representación de la Comisión Europea firma los contratos la ESA. Faltan por asignar otros contratos, pero al menos parece que se ha superado el parón de tantos años, y quizá lo que tenía que estar operativo en teoría este año (una teoría muy pasada de moda), pueda con suerte estar en marcha en 2015.

código Octave

Como epílogo a la ecuación del tiempo copiaré el código Octave usado para obtener las gráficas de las entradas de la ecuación del tiempo, duración del día solar y analema. Además ese código genera varias tablas con los datos numéricos calculados. La ruta empleada debería ajustarse según dónde tenga cada uno instalado Octave. Hay un fichero principal, EDTinvocar_todo.txt, que llama a los demás. En Octave la instrucción source permite cargar y ejecutar un fichero. Basta hacer el source del fichero principal, y las gráficas van saliendo con pausas entre ellas. Como mi experiencia con Octave es limitada, es muy posible que la potencia de Octave se pueda usar de maneras mucho mejores, pero en el tema del plotting quizá sirvan los ejemplos, que en el manual de Octave tampoco es que abunden, aunque es "la" referencia a consultar siempre. El código puede usarse libremente bajo la única responsabilidad de cada uno (at your own risk), sin ninguna garantía de que sea útil ni de que sea correcto, ni de que no perjudique la salud física o mental. Disclaimer total.
Por cierto, para ver código fuente en el navegador, con las líneas largas ocupando una sola línea, es muy recomendable desactivar los estilos (suele ser en el menú Ver/Estilo).

EDTinvocar_todo.txt

#source("c:/octave/EDT/EDTinvocar_todo.txt");
source("c:/octave/EDT/EDTfunciones.txt");
source("c:/octave/EDT/EDTcalculo_principal.txt");
source("c:/octave/EDT/EDTgrafica_EDT_EEX_EIE.txt");pause(5);
source("c:/octave/EDT/EDTgrafica_lambda_alfa_gr.txt");pause(5);
source("c:/octave/EDT/EDTgrafica_M_nu_gr.txt");pause(5);
source("c:/octave/EDT/EDTgrafica_EDT_DiaSolar.txt");pause(5);subplot(1,1,1);
source("c:/octave/EDT/EDTgrafica_analema.txt");pause(2);

EDTfunciones.txt

#source("c:/octave/ecuTiempo.txt")
format long
global incliECLI=23.43928108*pi/180 #23°26'21".4119 (sin(incliECLI)=0.397776995) inclinación de la eclíptica
global anyotropico=365.242190402 # en días solares medios
global anyosidereo = 365.256363004 # en días solares medios, en la época J2000
global excentricidadTierra = 0.016729 

global factorEDTsegrad = (24*60*60)/(2*pi) # escala vertical dif. angular (rad) <-> tiempo (s)
global factorEDTminrad = (24*60)/(2*pi)    # escala vertical dif. angular (rad) <-> tiempo (min)
global omega_T = 2*pi / anyotropico # fecuencia angular del Sol Ecuatorial Medio, sistema equinoccial, rad/dia
global omega_S = 2*pi / anyosidereo  # fecuencia angular de la anomalía media M, sistema orbital, rad/dia
global omega_P = omega_S*(1-(anyotropico/anyosidereo)); #frecuencia angular de precesión de los equinoccios, rad/dia

#{ 
La anomalía Excentrica E se calcula en base a la anomalía Media M, 
y a la excentricidad e de la órbita, mediante la ecuación de Kepler
M = E - e sen(E)
#}
function ret = anomaliaExcentrica (anomaliaMedia, excentricidad)
 fsol = @(x) (x - excentricidad *sin(x) - anomaliaMedia);
 ret = fsolve(fsol,anomaliaMedia);
endfunction

#{ 
La anomalía verdadera, nu, se relaciona con la anomalia Excentrica E,
y la excentricidad e de la órbita, mediante 
tan(nu/2)=sqrt( tan(E/2) * (1+e)/(1-e)  )
Aqui se obtiene en función de la anomalía media M.
nu-M es periodica de periodo 2*pi. 
Si M0=M-n*2*pi tal que 0<=M0<2*pi entonces 
nu = nu0-M0+M = nu0 + n*2*pi
#}
function retv = anomaliaVerdadera (M, exc)
 retv = zeros(1,length(M));
 for i = [1:length(M)]
  n = sign(M(i))*abs(floor(M(i)/(2*pi)));
  M0=M(i)-n*2*pi; 
  E = anomaliaExcentrica(M0, exc);
  if (E == pi) ret0 = pi;
  elseif (E < pi)  ret0 = 2*atan(sqrt((1+exc)/(1-exc))*tan(E/2));
  elseif (E > pi)  E = 2*pi-E; ret0 = 2*(pi - atan(sqrt((1+exc)/(1-exc))*tan(E/2)));
  endif
  retv(i) = ret0 + n*2*pi;
 endfor
endfunction

#{
Anomalía media correspondiente a una anomalía verdadera dada. Función inversa de anomaliaVerdadera().
#}
function retv = MdeAnomaliaVerdadera (nu, exc)
 retv = zeros(1,length(nu));
 for i = [1:length(nu)]
  n = sign(nu(i))*abs(floor(nu(i)/(2*pi)));
  nu0=nu(i)-n*2*pi; 
  if (nu0 == pi) E = pi;
  elseif (nu0 < pi)  E = 2*atan(sqrt((1-exc)/(1+exc))*tan(nu0/2));
  elseif (nu0 > pi)  E = 2*(pi - atan(sqrt((1-exc)/(1+exc))*tan((2*pi-nu0)/2)));
  endif
  M = E - exc * sin(E);
  retv(i) = M + n*2*pi;
 endfor
endfunction

#{ 
Devuelve 'y' (alfa) Ascension Recta de un punto en la eclíptica con Longitud Eclíptica 'x' (lambda) 
y una inclinación 'incli' de la ecliptica.
Por tanto x medida desde el punto de Aries
La diferencia  x-y es una función periodica de periodo pi, impar respecto a 0, par respecto a pi/4. 
Si x = X+n*pi, con 0<=X<=pi
x-y(x) = X-y(X)  --> y(x)=y(X)+n*pi   ,  n = sign(x)*abs(floor(x/pi))
#}
function retv = ARdeLE(xv, incli)
 retv = zeros(1,length(xv));
 for i = [1:length(xv)]
  x = xv(i);
  n = sign(x)*abs(floor(x/pi));
  x0=x-n*pi; # x0 es X 
  ret = acos(cos(x0)/sqrt((cos(x0))^2+(sin(x0)*cos(incli))^2)) + n*pi;
  retv(i)=ret;
 endfor
endfunction

#{
DeClinación (delta) de un punto con Longitud ecliptica 'x' (lambda), para inclinación 'incli' de la eclíptica.
sin(delta)=sin(lambda)*sin(incli)
#}
function retv = DCdeLE(xv, incli)
 retv = zeros(1,length(xv));
 for i = [1:length(xv)]
  x = xv(i);
  ret = asin(sin(x)*sin(incli));
  retv(i)=ret;
 endfor
endfunction

EDTcalculo_principal.txt

#invocada desde EDTinvocarTodo.txt
#{
d: tiempo medido en dias solares medios, desde el 1 de enero de 2009 (d=0 a las 0h de ese día)
sufijo _e para valores en el equinoccio de marzo y sufijo _p para valores en el perihelio
nu anomalía verdadera del Sol, medida en sistema orbital, origen en el perihelio
L  Longitud eclíptica del Sol, medida en sistema equinoccial, origen en el punto de Aries
AR_Sol ascension recta del Sol, medida en sistema equinoccial, origen en el punto de Aries
M anomalía media, medida en sistema orbtial, origen en perihelio. Unida al discurrir uniforme del tiempo, y a ECLI.
lamb longitud eclíptica de ECLI (igual en valor a la ascensión recta de ECUA), sistema equinoccial, origen en el punto de Aries.
AR_lamb ascensión recta de ECLI, sistema equinoccial, origen en el punto de Aries
#}
d_ini=-10; # incluiremos unos pocos dias del año anterior 
d_fin=370; # y otros pocos del siguiente
dias=[d_ini:d_fin]; 

# En 2009 el equinoccio de primavera se producirá a las 11h 44m UT del 20 de marzo
d_e=31+28+19+((11+(44/60))/24);

# En http://aom.giss.nasa.gov/srorbpar.html obtenemos para la distancia angular desde el perihelio al nodo ascendente, en grados
AOM2009_gr = 283.05
# El angulo medido (en sistema orbital) desde el perihelio hasta el punto de Aries, en radianes, 
nu_e = (360 - AOM2009_gr)*pi/180;
# A partir de este angulo determinamos cuando se pasa por el perihelio, d_p.
# Este angulo coincide con nu_e, la anomalia verdadera en el equinoccio de marzo, cuando el Sol pasa por el punto de Aries
# La anomalia media en el equinoccio es 
M_e = MdeAnomaliaVerdadera(nu_e, excentricidadTierra);
# y M_e=(d_e-d_p)*omega_S de donde
d_p = d_e - (M_e/omega_S);
# anomalia media 
M = omega_S*(dias-d_p);

nu_e = anomaliaVerdadera(M_e, excentricidadTierra); 
nu = anomaliaVerdadera(M, excentricidadTierra); 
# La diferencia angular en d= d_e entre ECLI y Sol en el S.R. orbital es igual a(M_e - nu_e), 
# y coincide con la diferencia angular vista desde el S.R. equinoccial (lamb_e - L_e), siendo L_e=0
# Por tanto lamb_e es negativo, e igual a
lamb_e=M_e-nu_e 
lamb_p = lamb_e + omega_T*(d_p - d_e);

# es importante recordar que lamb se refiere a ECLI visto en el S. R. Equinoccial, moviendose al ritmo omega_T.
lamb = lamb_e + omega_T * (dias - d_e);

# Precisamente en el S.R. Ecuatorial la AR (ascensión recta) de ECUA es igual a lamb de ECLI
AR_ECUA = lamb;

# La longitud ecliptica del Sol, L, seria (nu - nu_e) de no haber precesión. Pero como el equinoccio se mueve 
# en sentido contrario al Sol, hay un pequeñisimo incremento de L adicional de modo que 
L = nu - nu_e + omega_P * (dias - d_e);
L_e = 0;
L_p = 0 - nu_e + omega_P * (d_p - d_e);

#ascensión recta del Sol
AR_Sol = ARdeLE(L, incliECLI);
AR_Sol_e = ARdeLE(L_e, incliECLI);
AR_Sol_p = ARdeLE(L_p, incliECLI);

#ECUACION DEL TIEMPO EDT
valorEDT = (AR_ECUA - AR_Sol) * factorEDTminrad;
# ascensión recta de ECLI
AR_ECLI = ARdeLE(lamb, incliECLI);
#Ecuación de la inclinación de la eclíptica
valorEIE = -(AR_ECLI - AR_ECUA) * factorEDTminrad;
#Ecuación del centro o de la excentricidad de la órbita
valorEEX = -(AR_Sol - AR_ECLI) * factorEDTminrad;

EDTgrafica_EDT_EEX_EIE.txt

#invocada desde EDTinvocar_todo.txt
plot(dias, valorEDT,"-k;EDT;","linewidth",3, dias,valorEIE,"-b;EIE;",dias,valorEEX,"-r;EEX;",d_e,0,"+b",d_p,0,"+r");
xlabel("días desde el inicio del año 2009");
ylabel("minutos");
title("Gráfica 1 - ecuación del tiempo EDT");
grid("minor");
set (gca (), "xlim", [d_ini-5, d_fin+5]);
set(gca(),"position",[0.15,0.15,0.75,0.75])
legend("show","location","northwest")
set(text(36,1,"equinoccio"),"color","blue")
set(text(0,2,"perihelio"),"color","red")

fid = fopen ("c:/Octave/EDT/tabla_EDT_EEX_EIE.txt", "w");
fdisp(fid,"indice,d,fecha,EDT,EEX,EIE");
for i = [1:length(dias)]
 fprintf(fid,"%d,%d,%s,%f,%f,%f\n",i,dias(i),datestr(datenum(2009,1,1)+dias(i),29),valorEDT(i),valorEEX(i),valorEIE(i));
endfor
fclose(fid);

EDTgrafica_lambda_alfa_gr.txt

#invocada desde EDTinvocar_todo.txt
lambd=[0:360];
alfa = ARdeLE(lambd*pi/180, incliECLI)*180/pi;

plot(lambd, lambd-alfa,"-k",0,0,"+k");
xlabel('{\fontsize{12} \lambda (grados)}');
ylabel('{\fontsize{12} \lambda - \alpha  (grados)}');
title("Gráfica 2");
text(60,2.75,'{\fontsize{12}oblicuidad de la eclíptica = 23,439º}');
text(0,-0.25,'{\fontsize{12}punto de Aries}');
grid("minor");
set (gca (), "xlim", [-5, 365]);
set(gca(),"position",[0.15,0.15,0.75,0.75])
#print -dpng "c:/octave/EDT/grafica_lambda_alfa_gr.png";

fid = fopen ("c:/Octave/EDT/tabla_lambda_alfa_gr.txt", "w");
fdisp(fid,"lambda_gr,alfa_rad,lambda_rad,alfa_gr");
for i = [1:length(lambd)]
 fprintf(fid,"%d,%f,%f,%f\n",lambd(i),alfa(i),lambd(i)*pi/180,alfa(i)*pi/180);
endfor
fclose(fid);

EDTgrafica_M_nu_gr.txt

#invocada desde EDTinvocar_todo.txt
M=[0:360];
nu = anomaliaVerdadera (M*pi/180, excentricidadTierra)*180/pi;

plot(M, M-nu,"-k",0,0,"+k");
xlabel('{\fontsize{12} M (grados)}');
ylabel('{\fontsize{12} M - \nu (grados)}');
title("Gráfica 3");
text(50,1.5,'{\fontsize{12}excentricidad = 0,016729}');
text(0,0.2,'{\fontsize{12}perihelio}');
grid("minor");
set (gca (), "xlim", [-5, 365]);
set(gca(),"position",[0.15,0.15,0.75,0.75])
#print -dpng "c:/octave/EDT/grafica_M_nu_gr.png";

fid = fopen ("c:/Octave/EDT/tabla_M_nu_gr.txt", "w");
fdisp(fid,"M_gr,nu_rad,M_rad,nu_gr");
for i = [1:length(M)]
 fprintf(fid,"%d,%f,%f,%f\n",M(i),nu(i),M(i)*pi/180,nu(i)*pi/180);
endfor
fclose(fid);

EDTgrafica_EDT_DiaSolar.txt

#invocada desde EDTinvocar_todo.txt
difer = zeros(1,length(dias));
for i = [2:length(difer)]
 difer(i)=valorEDT(i)-valorEDT(i-1);
endfor
difer(1)=difer(366);
difer=-difer*60;#diferencia dia a dia en segundos CAMBIADA DE SIGNO
diasEnMes=[31,28,31,30,31,30,31,31,30,31,30,31];
nombreMes=["Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dic","Ene"];
diaIniMes=[0,cumsum(diasEnMes)];

subplot(2,1,1);
plot(dias, valorEDT,diaIniMes,valorEDT(-d_ini+1+diaIniMes),'+');
set (gca (), "xlim", [d_ini-5, d_fin+5]);
ylabel("EDT (minutos)");
#set(text(50,18,"ecuación del tiempo"),"color","black");
for i=[0:12] set(text(diaIniMes(i+1)-4,valorEDT(-d_ini+1+diaIniMes(i+1))+5,nombreMes(3*i+1:3*i+3)),"color","black") ;endfor

grid("on");
subplot(2,1,2);
plot(dias,difer,diaIniMes,difer(-d_ini+1+diaIniMes),'+');
set (gca (), "xlim", [d_ini-5, d_fin+5]);

xlabel("días desde el inicio del año 2009");
ylabel(strcat('\Delta D_v (segundos)'));
set(text(52,26,"exceso del día solar verdadero sobre el día solar medio"),"color","black")
set(text(52,23,"valores positivos: días más largos, negativos: días más cortos"),"color","black")
for i=[0:12] set(text(diaIniMes(i+1)-4,difer(-d_ini+1+diaIniMes(i+1))-5,nombreMes(3*i+1:3*i+3)),"color","black") ;endfor
grid("on");

fid = fopen ("c:/Octave/EDT/tabla_EDT_DiaSolar.txt", "w");
fdisp(fid,"indice,d,fecha,EDT_min,difer_s");
for i = [1:length(dias)]
 fprintf(fid,"%d,%d,%s,%f,%f\n",i,dias(i),datestr(datenum(2009,1,1)+dias(i),29),valorEDT(i),difer(i));
endfor
fclose(fid);

EDTgrafica_analema.txt

#invocada desde EDTinvocar_todo.txt
#declinacion del Sol en grados
DC_Sol = DCdeLE(L, incliECLI)*180/pi;
DC_Sol_e = DCdeLE(L_e, incliECLI)*180/pi;
DC_Sol_p = DCdeLE(L_p, incliECLI)*180/pi;

# ángulo horario: diferencia angular entre la ascension recta del Sol y de ECUA, Horizontal, en grados
xH=(lamb-AR_Sol)*180/pi;
xH_e=(lamb_e-AR_Sol_e)*180/pi;
xH_p=(lamb_p-AR_Sol_p)*180/pi;

# plot(dias,DC_Sol);pause(2);#print -dpng "c:/octave/EDT/grafica_dias_DC_gr.png";
diasEnMes=[31,28,31,30,31,30,31,31,30,31,30,31];
nombreMes=["Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dic","Ene"];
diaIniMes=[0,cumsum(diasEnMes)];

#plot(xH,DC_Sol,'-k',"linewidth",1,xH_p,DC_Sol_p,'+r',xH_e,DC_Sol_e,'+b')
latN=0 #latitud Norte del punto de observacion
plot(xH,DC_Sol+latN,'-k',"linewidth",1,xH_p,DC_Sol_p+latN,'+r',xH_e,DC_Sol_e+latN,'+b',xH(-d_ini+1+diaIniMes),DC_Sol(-d_ini+1+diaIniMes),'+g')
offx=[-1,-3,-3,-3,+1,+1,-3,-3,+1,+1,+1,+1];
offy=[-1.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5];
for i=[0:11] set(text(xH(-d_ini+1+diaIniMes(i+1))+offx(i+1),DC_Sol(-d_ini+1+diaIniMes(i+1))+offy(i+1),nombreMes(3*i+1:3*i+3)),"color","black") ;endfor

xlabel("Ángulo entre Sol y meridiano (grados)");
ylabel("Declinación del Sol (grados)");
title("Analema en el Ecuador");
grid("on");
#grid("minor");
set (gca (), "xlim", [-15, 15]);
set (gca (), "ylim", [-30+latN, 30+latN]);
set (gca (), "dataaspectratio", [1, 2, 1]);
set(gca(),"position",[0.15,0.15,0.75,0.75])
set(text(-9,0.5+latN,"equinoccio"),"color","blue")

set(text(-7,-23+latN,"perihelio"),"color","red")
set(text(10,0.5+latN,"OESTE"),"color","black")
set(text(-14,0.5+latN,"ESTE"),"color","black")
set(text(-2,27+latN,"NORTE"),"color","black")
set(text(-1,-27+latN,"SUR"),"color","black")
#print -dpng "c:/octave/EDT/grafica_analema.png";
pause(5);
#-----------------------------------------------------------------------------------
# grafica de los analemas a diferentes horas vistos desde una latitud no ecuatorial
source("c:/octave/EDT/EDTfunciones_cambio_coordenadas.txt")
lat=40; #latitud de la observacion
for i=[0:23]
 [z,A]= zAhoriz(DC_Sol*pi/180,(xH+15*i)*pi/180,lat*pi/180);z=z*180/pi;A=A*180/pi;
 zz{i+1}=90-z;AA{i+1}=A;#zz distancia al horizonte en vez de al cenit, en grados
endfor
plot(AA{1},zz{1},"k",AA{2},zz{2},AA{3},zz{3},AA{4},zz{4},AA{5},zz{5},AA{6},zz{6},AA{7},zz{7},AA{8},zz{8},AA{9},zz{9},AA{10},zz{10},AA{11},zz{11},AA{12},zz{12},AA{13},zz{13},AA{14},zz{14},AA{15},zz{15},AA{16},zz{16},AA{17},zz{17},AA{18},zz{18},AA{19},zz{19},AA{20},zz{20},AA{21},zz{21},AA{22},zz{22},AA{23},zz{23},AA{24},zz{24});
title("analema cada hora a 40º latitud Norte");
xlabel("Azimut desde el Sur (+ Oeste, - Este)(grados)");
ylabel("Altura sobre el horizonte (grados)");
grid("on");
set (gca (), "xlim", [-180, 180]);
set (gca (), "ylim", [-90, 90]);
#print -dpng "c:/octave/EDT/grafica_analemas40N.png";
pause(5);
#-----------------------------------------------------------------------------------
lat=-40; #latitud de la observacion
for i=[0:23]
 [z,A]= zAhoriz(DC_Sol*pi/180,(xH+15*i)*pi/180,lat*pi/180);z=z*180/pi;A=A*180/pi;
 zz{i+1}=90-z;AA{i+1}=A;#zz distancia al horizonte en vez de al cenit, en grados
endfor
plot(AA{1},zz{1},"k",AA{2},zz{2},AA{3},zz{3},AA{4},zz{4},AA{5},zz{5},AA{6},zz{6},AA{7},zz{7},AA{8},zz{8},AA{9},zz{9},AA{10},zz{10},AA{11},zz{11},AA{12},zz{12},AA{13},zz{13},AA{14},zz{14},AA{15},zz{15},AA{16},zz{16},AA{17},zz{17},AA{18},zz{18},AA{19},zz{19},AA{20},zz{20},AA{21},zz{21},AA{22},zz{22},AA{23},zz{23},AA{24},zz{24});
title("analema cada hora a 40º latitud Sur");
xlabel("Azimut desde el Norte (+ Oeste, - Este)(grados)");
ylabel("Altura sobre el horizonte (grados)");
grid("on");
set (gca (), "xlim", [-180, 180]);
set (gca (), "ylim", [-90, 90]);
#print -dpng "c:/octave/EDT/grafica_analemas40S.png";
pause(5);
#-----------------------------------------------------------------------------------
lat=0; #latitud de la observacion
for i=[0:23]
 [z,A]= zAhoriz(DC_Sol*pi/180,(xH+15*i)*pi/180,lat*pi/180);z=z*180/pi;A=A*180/pi;
 zz{i+1}=90-z;AA{i+1}=A;#zz distancia al horizonte en vez de al cenit, en grados
endfor
plot(AA{1},zz{1},"k",AA{2},zz{2},AA{3},zz{3},AA{4},zz{4},AA{5},zz{5},AA{6},zz{6},AA{7},zz{7},AA{8},zz{8},AA{9},zz{9},AA{10},zz{10},AA{11},zz{11},AA{12},zz{12},AA{13},zz{13},AA{14},zz{14},AA{15},zz{15},AA{16},zz{16},AA{17},zz{17},AA{18},zz{18},AA{19},zz{19},AA{20},zz{20},AA{21},zz{21},AA{22},zz{22},AA{23},zz{23},AA{24},zz{24});
title("analema cada hora en el ecuador");
xlabel("Azimut desde el Sur (+ Oeste, - Este)(grados)");
ylabel("Altura sobre el horizonte (grados)");
grid("on");
set (gca (), "xlim", [-180, 180]);
set (gca (), "ylim", [-90, 90]);
#print -dpng "c:/octave/EDT/grafica_analemasEcuador.png";

fid = fopen ("c:/Octave/EDT/tabla_analema.txt", "w");
fdisp(fid,"indice,d,fecha,difAR_gr,DC_gr,AR_Sol_gr,AR_ECUA_gr");
for i = [1:length(dias)]
 fprintf(fid,"%d,%d,%s,%f,%f,%f,%f\n",i,dias(i),datestr(datenum(2009,1,1)+dias(i),29),xH(i),DC_Sol(i),AR_Sol(i)*180/pi,AR_ECUA(i)*180/pi);
endfor
fclose(fid);

EDTfunciones_cambio_coordenadas.txt

#{
paso de coordenadas ecuatoriales declinación delta, y ángulo horario t,
a coordenadas horizontales distancia cenital z y Azimut A, para una latitud lat
delta y t pueden ser vectores, pero de la misma dimensión. lat se considera escalar.
#}
function [z, A] = zAhoriz(delta, t, lat)
 z = zeros(1,length(delta))';
 A = zeros(1,length(delta))';
 slat=lat;#valor original con signo
 lat=abs(lat);#valor absoluto
 for i = [1:length(z)]
  if (slat<0) t(i)=t(i)+pi; #para latitudes negativas hay desfase pi en el angulo horario
  endif
  if (t(i)<0 || t(i)>2*pi) 
   t(i)=mod(t(i),2*pi); # reducimos SIEMPRE el ángulo horario t al intervalo [0,2*pi]
  endif
  z(i)=acos(sin(lat)*sin(delta(i))+cos(lat)*cos(delta(i))*cos(t(i)));
  
  if (abs(z(i)) < 0.0001) sA=0; # para z=0 el valor de A es indefinido, elegimos A=0.
  else sA=cos(delta(i))*sin(t(i))/sin(z(i)); # sA es sin(A)
       cA=(-cos(lat)*sin(delta(i))+sin(lat)*cos(delta(i))*cos(t(i)))/sin(z(i)); #cA es cos(A)
  endif
  # para obtener A(i) en [0,2*pi]
  # cuadrante coseno seno para angulo inverso usar
  #     1       +     +    acos o asin
  #     2       -     +    acos o pi-asin
  #     3       -     -    2*pi-acos o pi-asin
  #     4       +     -    2*pi-acos o 2*pi+asin
  # En octave asin toma valores en [-pi/2,pi/2], y acos toma valores en [0,pi], para el dominio [-1,1].
  if (sA>0) A(i)=acos(cA);
  elseif (sA<0) A(i)=2*pi-acos(cA);
  else A(i)=0;
  endif
  #Pero nos interesa A(i) en [-pi,pi]
  if (A(i)>pi) A(i)=A(i)-2*pi;
  endif
  if (slat<0) z(i)=pi-z(i); 
   # ahora se deshace el desfase pi en latitudes negativas. 
   # Queremos seguir teniendo las 12h del mediodia a 0º de acimut, que ahora se cuenta desde el Norte, no desde el Sur
   # y asi la grafica de analemas para latitudes sur queda bien
   if (A(i)<=0) A(i)=A(i)+pi;
   else A(i)=A(i)-pi;
   endif
  endif
 endfor
endfunction

analema

Aquí va otra secuela de la anterior entrada sobre la ecuación del tiempo, EDT. En ella se citaba el analema, la figura que puede construirse sobre el cielo con una cámara, un tripode, un reloj, y paciencia. Se trata de fotografiar exactamente la misma zona del cielo, todos los días del año (cada semana puede servir también), a la misma hora de nuestro buen reloj, que marca las horas "normales" de tiempo ecuatorial medio.

Como alternativa a la fotografía, vamos a calcular la posición del Sol a una misma hora de cada día, eligiendo una especialmente adecuada, las 12h de un buen reloj ajustado a la hora local de nuestro meridiano.

La posición del Sol a esa hora puede determinarse por dos coordenadas angulares. Una es lo que se denomina ángulo horario, que en nuestro caso coincide con la EDT expresada como ángulo en vez de en minutos de tiempo (recordando que 24h corresponden a 360º). Ese ángulo es la separación, en ascensión recta, entre el Sol verdadero y nuestro meridiano local, pues justo a las 12h de reloj en ese meridiano culmina el Sol Ecuatorial Medio, ECUA. La otra coordenada necesaria es la declinación del Sol.

Mostramos aquí de nuevo la Figura 2 de la citada entrada, copiando la explicación de la misma

• A es nuestro origen de medida de arcos en el sistema equinoccial, el punto de Aries.
• AB es un arco sobre la eclíptica, cuya magnitud λ da la longitud eclíptica del punto B.
• AC es un arco sobre el ecuador, cuya magnitud α da la ascensión recta del punto B, igual a la del punto C, que es la proyección de B sobre el ecuador, dado que el arco BC es parte de un meridiano celeste.
• Los lados AB y AC forman un ángulo ε igual a la inclinación (u oblicuidad) de la eclíptica.

Podemos añadir ahora que la declinación δ del punto B es el ángulo correspondiente al arco BC. En nuestro caso el punto B es el ocupado por el Sol, y la posición de este se fija inequívocamente con las dos coordenadas α, ascensión recta, y δ, declinación.

La trigonometría esférica nos proporcionaba, como se vió, una relación entre α y λ :



Pues bien, también nos ofrece esta otra relación entre δ y λ:



De modo que a los cálculos para obtener EDT, que nos daban α, les añadimos una función más para calcular δ en función de λ, que dependía a su vez del tiempo, y tenemos las dos coordenadas del Sol para cada día del año, a las 12h local, lo que nos lleva a la siguiente gráfica del Analema.



Esta típica figura de 8 desigual corresponde al analema visto desde un lugar en el ecuador terrestre. Nos podemos imaginar tumbados allí, en el suelo plano, boca arriba, con la cabeza hacia el Norte, los pies hacia el Sur, el brazo derecho extendido hacia el Oeste, el brazo izquierdo extendido hacia el Este, divididos imaginariamente en dos mitades por la línea meridiana que corre de Norte a Sur. Sobre nosotros el cenit corresponde justamente al origen de coordenadas (0,0) en la gráfica anterior. En el eje de abscisas se representa el ángulo horario del Sol, antes mencionado. Respecto a la Tierra, la esfera celeste va de Este a Oeste. Si a las 12h de buen reloj local el Sol está al Oeste del meridiano, es que ya ha culminado, su hora solar verdadera es mayor de las 12h, y EDT es positiva. Cuando el Sol se sitúa al Este es cuando la EDT (y la abscisa en la gráfica) es negativa. En cuanto a la declinación, mostrada en el eje de ordenadas, se toma como positiva hacia el Norte y como negativa hacia el Sur.

Si en vez de estar en el ecuador imaginamos que nos deslizamos sobre el meridiano hacia el Norte hasta una latitud de 40º, el origen (0,0) se nos habrá movido respecto a nuestro nuevo cenit justo 40º hacia el Sur. Si quisieramos poner las coordenadas referidas al nuevo cenit, deberíamos cambiar los valores del eje de ordenadas. El punto de cruce del 8 del analema está en el ecuador a unos 9º al norte del cenit. En la latitud de 40º Norte, ese mismo punto estará a 31º al Sur del cenit de esa latitud, o a 59º sobre el horizonte Sur. En la latitud de 40º Sur, ese punto estaría a 49º al norte del cenit, o 41º sobre el horizonte Norte. Esto nos podría ayudar a estimar la latitud a partir de una fotografía de un analema.

Pero en las fotografías el analema suele salir inclinado respecto al horizonte, no perpendicular como es el caso mostrado en la gráfica. La "culpa" es dos dos factores. Por un lado las fotos suelen hacerse en latitudes medias septentrionales. Además se han tomado por lo general a horas tempranas o tardías, no a las 12h hora local. Eso facilita que el analema esté cerca del horizonte y salga en la foto junto a algún rasgo del paisaje.

El analema puede considerarse una figura "adherida" a la esfera celeste. El Sol se desliza por ella a lo largo del año, pero además cada día la rotación terrestre hace dar vueltas al analema desde el Este al Oeste por el cielo visible, y del Oeste al Este bajo el horizonte.

¿Cómo se vería el analema en una latitud de 40º Norte a lo largo de un día? Nuestra amiga la trigonometría esférica nos permite pasar de unas coordenas ecuatoriales, como la declinación y ángulo horario, a las llamadas coordenadas horizontales, que nos situan un punto respecto al terreno con dos coordenadas angulares: la distancia o ángulo cenital y el azimut. En lugar de la distancia cenital puede usarse la altura sobre el horizonte.

En las dos gráficas siguientes se muestra el analema cada hora del día para una latitud 40ºN y luego para otra 40º Sur. En el eje de abscisas se usa un azimut particular en cada caso, midiendo desde el Sur en la latitud Norte, y desde el Norte en la latitud Sur, para que el origen 0 de azimut corresponda a la hora del mediodía medio, las 12h de reloj local. Desde ese origen, en ambas gráficas los valores positivos de azimut corresponden a ir al Oeste (y a horas despues de las 12h), y los valores negativos de azimut a ir hacia el Este (y a horas antes de las 12h). En el eje de ordenadas se ha preferido, en vez de la distancia cenital, representar la altura sobre el horizonte, siendo éste el origen 0, con valores positivos para el hemisferio visible, y valores negativos para el hemisferio oculto bajo el horizonte.





Aunque estas gráficas dan una idea del aspecto del analema a cada una de las 24 horas del día, son representaciones cartesianas de dos coordenadas esféricas, que habría que deformar para "pegarlas" sobre la esfera celeste o colocarlas sobre una fotografía. Las coordenadas esféricas son especialmente malas cerca de los puntos singulares. Para las coordenadas horizontales (distancia cenital o altura y azimut) los puntos singulares son el cenit y el nadir, a una altura sobre el horizonte de +90 y -90º respectivamente, pero con azimut indeterminado, pues en la inmediata proximidad de ambos puntos hay otros con cualquier valor de azimut que queramos. Este "mal comportamiento" se ve de forma espectacular al representar la gráfica del mismo tipo que las dos anteriores, pero para el ecuador.



En la primera gráfica de esta entrada el analema se representaba en unas coordenadas ecuatoriales cuyos puntos singulares son los polos, pero que se portan muy bien cerca del cenit. Sin embargo en esta última gráfica la figura del analema a las 12h es casi irreconocible, pues está en las proximidades del cenit, punto singular para las coordenadas horizontales. Sin embargo lo que se visualiza muy bien aquí es que la salida y la puesta del analema sobre el horizonte es perpendicular al mismo. En el ecuador el eje de giro terrestre es paralelo a la línea meridiana, y todos los objetos celestes ascienden rectos (perpendiculares) sobre el horizonte al Este, describen un semicirculo y se ponen por el horizonte al Oeste, formando también un ángulo recto con él. En esa latitud a la esfera celeste se le denomina esfera recta. Precisamente la ascensión recta se mide sobre esos circulos en que los astros dan vueltas, perpendiculares al horizonte ecuatorial.

Duración del día solar

Tras la entrada dedicada a la ecuación del tiempo pensaba emplear mi idem en otra cosa, pero hay alguna otra información útil extraible de mi querida EDT. Por ejemplo cómo varía la duración del día solar verdadero respecto al día solar medio.

Usando t_v para tiempo solar verdadero, t_m para el tiempo solar medio, la ecuación del tiempo es

EDT= t_v - t_m

Podemos poner esto mismo así

t_m = t_v - EDT

Consideremos dos culminaciones sucesivas sobre nuestro meridiano del Sol verdadero, de un mediodía al siguiente. En la primera culminación tenemos que

t_m1 = t_v1 - EDT₁    (1)

y en la siguiente, un día solar verdadero despúes

t_m2 = t_v2 - EDT₂    (2)

En ambos casos t_v1=t_v2=12h, es el mediodía, y la diferencia entre los tiempos de "reloj de 24h" de ambas culminaciones es justo la diferencia de los valores de EDT, cambiada de signo:

t_m2-t_m1 = - (EDT₂-EDT₁)    (3)


Planteado así, a (2) le hemos restado (1) y tomado t_v1=t_v2=12h con lo que al restar ese término desaparece y llegamos a (3). Esto es correcto tomando valores de tiempos entre 0 y 24h, que es lo que se ha supuesto.

Adoptando alternativamente un criterio de cómputo continuo del tiempo (sin vuelta a cero a las 24h), si en la primera culminación el tiempo solar verdadero era de 12h "verdaderas", en la siguiente (un día solar verdadero después) será de 12+24=36h "verdaderas". También el tiempo solar medio se incrementará casi 24h (de las "normales"). Lo que nos interesa en este caso no es t_m2-t_m1, sino t_m2-t_m1-24h. El resultado final será el mismo. Conviene tener presente que el valor de t_m y de EDT se expresa en unidades como días, horas o segundos "normales", mientras que el valor de t_m se expresa en unidades "verdaderas", que no usamos en la práctica precisamente por ser variables.

En definitiva, la diferencia de la duración del día solar verdadero D_v respecto al día solar medio D_m, está dada por la ecuación (3):

ΔD_v = D_v -D_m = - (EDT₂-EDT₁)

Por tanto de nuestra ecuación del tiempo sacamos cómo son de largos o cortos los días solares verdaderos respecto a las 24h de reloj. Podemos hacer la diferencia de la EDT en un día con respecto al día anterior, para cada día. Esto es muy aproximadamente la derivada de EDT respecto al tiempo, cuando este se mide en unidades de día solar medio, que es como se representa EDT habitualmente. Por tanto ΔD_v es la derivada de EDT respecto al tiempo, cambiada de signo. Con nuestro Octave podemos ver en la siguiente gráfica la EDT arriba y debajo ΔD_v.



En efecto se aprecía que el máximo de ΔD_v corresponde a la mayor pendiente negativa de EDT, el mínimo de ΔD_v corresponde a la mayor pendiente positiva de EDT, y que ΔD_v se anula cuando EDT alcanza sus máximos o mínimos.

Los días solares más cortos, con unos 21 segundos de menos, se dan hacia el 17 de septiembre, y desde entonces van creciendo rápidamente en duración hasta su máximo hacia el 23 de diciembre en que llegan a casi 30 segundos de más. En estas fechas se encuentra la máxima pendiente de EDT positiva y negativa respectivamente. La contribución mayor a la pendiente de EDT es de EIE (ver la figura de la entrada anterior), pero además en esta zona del año, sobre todo al final, EEX suma su variación a la de EIE de forma muy acusada, lo que propicia el máximo de ΔD_v aún mayor en valor absoluto que el mínimo.

Hay cuatro fechas en que ΔD_v se anula, y por tanto el día solar verdadero dura lo mismo que el día solar medio: el 11 de febrero, el 14 de mayo, el 26 de julio y el 3 de noviembre. Pero esto no implica ni mucho menos que coincidan los tiempos solar verdadero y solar medio, lo que sucede cuando EDT=0, circunstancia que se da el 15 de abril, el 13 de Junio, el 1 de septiembre y el 25 de diciembre.