Calculos Astronomicos - Capitulo 1

El texto que a continuación se presenta explica como se computa la posición en el cielo del Sol, la Luna, los planetas, cometas y asteroides.

Las fórmulas para estos cálculos pueden ser complicadas, pero si se toman en cuenta algunas simplificaciones es posible obtener unos cálculos sencillos con una exactitud de un minuto de arco = 1/60 grados. Para los usos corrientes del aficionado a la astronomía, son bastante buenos la exactitud aquí presentada.

Las simplificaciones realizadas aquí son las siguientes:

1.Nutación es ignorada.
2.Aberración planetaria es ignorada.
3.Las diferencia entre Tiempo Terrestre, Tiempo Efemérides y Tiempo Universal son ignorados.
4.La Precesión es calculada de una forma simplificada, es decir, agregándole simplemente un ángulo de corrección a la longitud de la eclíptica.
5.Términos de alto orden en las orbitas planetarias son ignorados. Esto nos dará un error de 2 minutos de arco (0º 2') en 1.000 años desde ahora.
6. Algunas perturbaciones planetarias son ignoradas. Sólo las mayores perturbaciones para la Luna, Júpiter, Saturno y Urano son incluidas. Si no importa mantener una exactitud elevada estas perturbaciones también pueden ser ignoradas. Índice

Elementos orbitales

Consisten en 6 cantidades las cuales definen completamente una orbita circular, elíptica, parabólica o hiperbólica. Tres de estas cantidades describen la forma, tamaño y posición del planeta en la orbita:

a: Distancia media

e: excentricidad

T: Tiempo en el perihelio

Los tres restantes elementos definen la orientación de la orbita en el espacio:

i: Inclinación de la orbita con relación a la eclíptica.

N: Longitud del nodo ascendente.

W: Es el ángulo desde el nodo ascendente al perihelio a lo largo de la orbita.

Existen otros elementos que se obtienen a partir de los valores anteriores que son usadas en los cálculos que son los siguientes:

AU: Unidad Astronómica es la distancia media de la Tierra al Sol.

Para describir la posición en la orbita, se usa tres ángulos Anomalía Media, Anomalía Verdadera y Anomalía Excéntrica. Ellas son cero cuando el planeta esta en el perihelio.

M: Anomalía Media, este ángulo se incrementa uniformemente hasta 360º por período orbital.

v: Anomalía Verdadera , es el actual ángulo entre el planeta y el perihelio como es visto desde el Sol. Este ángulo no se incrementa uniformemente, sino que cambia más rápidamente en el perihelio.

E: Anomalía Excéntrica, este es un ángulo auxiliar usado por la ecuación de Kepler.

w1 = N + w = Longitud del perihelio

L = M + w1 = Longitud media

q = a * (1 - e) = Distancia del perihelio

Q = a * (1 + e) = Distancia del afelio

P = a^1,5 = período orbital (en años)

T = época de M - (M / 360) / P = Tiempo del perihelio Índice


Escala del tiempo

La escala del tiempo está formulado para contar días. Las horas, minutos y segundos son expresados en fracción de días. El día cero comienza el 31 de diciembre de 1999 a las 12:00 am UT (Universal Time, es decir, Tiempo Universal).

d= 367*y - int( (7 * (y +(int( (m + 9) / 12) ) ) ) / 4) + int( 275*m / 9 ) + D - 730530 +UT / 24

Donde y = año (cuatro dígitos); m = mes; D = día y UT = en horas + decimales. int () es una función que sólo toma la parte entera de la división. Índice

Oblicuidad de la eclíptica

Otra cantidad que es necesaria es la oblicuidad de la eclíptica que es la desviación del eje de rotación de la Tierra y que actualmente es 23,4º que se encuentra decreciendo lentamente. La fórmula es:

ecl = 23,4393 - 0,0000003563*d

Índice
Corrección por Precesión

Simplificaremos la corrección por precesión agregando simplemente la siguiente fórmula a la longitud de la eclíptica:

3,82394E-5 * ( 365,2422 * ( época - 2000.0 ) - d )

Donde época es el año y fracción al cual vamos hacer los cálculos y d ya fue realmente calculado en la escala del tiempo. El resultado de esta fórmula es en grados y se agrega a la fórmula de la eclíptica.

Índice
Elementos Orbitales

Todos los ángulos son expresados en grados y es posible que al efectuar los cálculos den resultados mayores a 360º o negativos, por lo que es necesario aplicar la siguiente función ang = ang - int( ang / 360 ) * 360 a los elementos orbitales N, w y M.

Elementos Orbitales del Sol:

N = 0
i = 0
w = 282,9404 + 4,70935E-5 * d
a = 1 (AU)
e = 0,016709 - 1,151E-9 * d
M = 356,0470 + 0,9856002585 * d

Elementos Orbitales de la Luna:

N = 125,1228 - 0,0529538083 * d
i = 5,1454
w = 318,0634 + 0,1643573223 * d
a = 60,2666 (radio Terrestre)
e = 0,0549
M = 115,3654 + 13,0649929509 * d

Elementos Orbitales de Mercurio:

N = 48,3313 + 3,24587E-5 * d
i = 7,0047 + 5,00E-8 * d
w = 29,1241 + 1,01444E-5 * d
a = 0,387098 (AU)
e = 0,205635 + 5,59E-10 * d
M = 168,6562 + 4,0923344368 * d

Elementos Orbitales de Venus:

N = 76,6799 + 2,4659E-5 * d
i = 3,3946 + 2,75E-8 * d
w = 54,891 + 1,38374E-5 * d
a = 0,72333 (AU)
e = 0,006773 - 1,302E-9 * d
M = 48,0052 + 1,6021302244 * d

Elementos Orbitales de Marte:

N = 49,5574 + 2,11081E-5 * d
i = 1,8497 - 1,78E-8 * d
w = 286,5016 + 2,92961E-5 * d
a = 1,523688 (AU)
e = 0,093405 + 2,516E-9 * d
M = 18,6021 + 0,5240207766 * d

Elementos Orbitales de Júpiter:

N = 100,4542 + 2,76854E-5 * d
i = 1,303 - 1,557E-7 * d
w = 273,8777 + 1,64505E-5 * d
a = 5,20256 (AU)
e = 0,048498 + 4,469E-9 * d
M = 19,895 + 0,0830853001 * d

Elementos Orbitales de Saturno:

N = 113,6634 + 2,3898E-5 * d
i = 2,4886 - 1,081E-7 * d
w = 339,3939 + 2,97661E-5 * d
a = 9,55475 (AU)
e = 0,055546 - 9,499E-9 * d
M = 316,967 + 0,0334442282 * d

Elementos Orbitales de Urano:

N = 74,0005 + 1,3978E-5 * d
i = 0,7733 + 1,9E-8 * d
w = 96,6612 + 3,0565E-5 * d
a = 19,18171 - 1,55E-8 * d (AU)
e = 0,047318 + 7,45E-9 * d
M = 142,5905 + 0,011725806 * d

Elementos Orbitales de Neptuno:

N = 131,7806 + 3,0173E-5 * d
i = 1,77 - 2,55E-7 * d
w = 272,8461 - 6,027E-6 * d
a = 30,05826 + 3,313E-8 * d (AU)
e = 0,008606 + 2,15E-9 * d
M = 260,2471 + 0,005995147 * d

El cálculo de la posición del Sol es como la de cualquier otro planeta e inclusive la Luna, pero como el Sol siempre se mueve en el plano de la eclíptica, más simplificaciones pueden ser hechas.

Primero calculamos la anomalía excéntrica (E) la cual hay que hacer una iteración hasta que la diferencia sea de por lo menos 0,001º :

E = M + e * (180/pi) * sin(M) * (1 + e * cos(M))

E0 = E

E1 = E0 - (E0 - e * (180/pi) * sin(E0) - M) / (1 - e * cos(E0))

Luego se hace E1 igual E0 hasta que sus valores no difieran en 0,0001º. Todos los valores tiene que estar expresado en grados.

Ahora calcularemos la distancia (r) y la anomalía verdadera (v):

xv = a * (cos(E) - e)

yv = a * (raiz(1 - e^2) *sin(E))

v = atan2(xv, yv)

r = raiz(xv^2 + xv^2)

La longitud verdadera (lon) se calcula con la siguiente fórmula:

lon = v + w

Índice
Coordenadas heliocéntricas

Ahora hay que obtener la posición del cuerpo en el espacio 3-D. Estas son las coordenadas heliocéntricas, es decir, el origen del sistema se encuentra en el Sol. Las fórmulas son las siguientes:

xh = r * (cos(N) * cos(v + w ) - sin(N) * sin(v +w) * cos(i))

yh = r * (sin(N) * cos(v + w) + cos(N) * sin(v+w) * cos(i))

zh = r * (sin(v + w) * sin(i))

La longitud (lonh), la latitud (lath) y el radio (rh), es decir, pasar estas a coordenadas esféricas, se calcula con la fórmula siguiente:

lonh = atan2(xh, yh)

lath = atan2(raiz(xh^2 + yh^2), zh)

rh = raiz(xh^2 + yh^2 + zh^2)

Índice
Posición geocéntrica de los planetas

Ya conocemos las coordenadas con el Sol en el centro, ahora para el aficionado a la astronomía es interesante saber donde se encuentran los cuerpos desde la Tierra como centro de coordenadas. Entonces tenemos que trasladar el sistema de coordenadas heliocéntricas al sistema de coordenadas geocéntricas, es decir, con el origen en la Tierra. Para convertir estas coordenadas nosotros simplemente agregamos las coordenadas heliocéntricas del Sol y le sumamos las coordenadas heliocéntricas del planeta.

xg = xh(sol) + xh(planeta)

yg = xh(sol) + yh(planeta)

zg = zh(sol) + zh(planeta)

Índice
Coordenadas ecuatoriales

Ahora necesitamos pasar estas coordenadas geocéntricas a coordenadas ecuatoriales, esto es debido a que la Tierra se encuentra inclinada sobre su eje de rotación y hay conocer en que parte del cielo se encuentra el planeta. Para conocer esto se usa la Ascensión Recta (RA también denominado alfa), la Declinación (Decl también denominado delta) y el Radio (r), que definen la posición de un cuerpo en el cielo (esto no es más que transformar las coordenadas rectangulares geocéntricas a coordenadas esféricas geocéntricas).

xe = xg

ye = yg * cos(ecl) - zg * sin(ecl)

ze = yg * sin(ecl) + zg * cos(ecl)

Finalmente calculamos RA, Decl y r:

RA = atan2(xe, ye)

Decl = atan2(raiz(xe^2 + ye^2), ze)

r = raiz(xe^2 + ye^2 + ze^2)

La RA, Decl y r son usados mayormente en los Almanaques de Astronomía para definir la posición de un cuerpo en el cielo.

Existe perturbaciones de planetas grandes como Júpiter, Saturno y Urano. La Luna también sufre perturbaciones por parte de la Tierra, por lo que a continuación se presentan varias correcciones para mantener nuestra precisión. Si en caso de que al lector no estuviera interesado en perder hasta 1 º, se puede obviar la siguientes correcciones. Índice


Perturbaciones con Júpiter, Saturno y Urano

Sólo estos planetas tienen perturbaciones mayores a 0,001º por lo que se hace necesaria una corrección para estar con la exactitud esperada para los otros planetas.

Mj: Anomalía media de Júpiter
Ms: Anomalía media de Saturno
Mu: Anomalía media de Urano

Agregar estos resultados (todo en grados) a la longitud de Júpiter

-0,332 * sin(2 * Mj - 5 * Ms - 67,6)
-0,056 * sin(2 * Mj - 2 * Ms + 21)
+0,042 * sin(3 * Mj - 5 * Ms + 21)
-0,036 * sin(Mj - 2 * Ms)
+0,022 * cos(Mj - Ms)
+0,023 * sin(2 * Mj - 3 * Ms + 52)
-0,016 * sin(Mj - 5 * Ms - 69)

Agregar estos resultados (todo en grados) a la longitud de Saturno

+0,812 * sin(2 * Mj - 5 * Ms - 67,6)
-0,229 * cos(2 * Mj - 4 * Ms -2)
+0,119 * sin(Mj - 2 * Ms - 3)
+0,046 * sin(2 * Mj - 6 * Ms - 69)
+0,014 * sin(Mj - 3 * Ms +32)

También agregar estos resultados (todo en grados) a la latitud de Saturno

-0,02 * cos(2 * Mj - 4 * Ms - 2)
+0,018 * sin(2 * Mj - 6 * Ms - 49)

Agregar estos resultados (todo en grados) a la longitud de Urano

+0,04 * sin(Ms - 2 * Mu + 6)
+0,035 * sin(Ms - 3 * Mu + 33)
-0,015 * sin(Mj - Mu + 20)

Entonces a la lonh, lath y rh hay que agregar los valores antes descritos respectivamente y volver a transformar el sistema esférico a rectangular, es decir xh, yh, zh (perturbado).

xh = rh * cos(lonh) * cos(lath)

yh = rh * sin(lonh) * cos(lath)

zh = rh * sin(lath)

Índice
Cálculo de Plutón

No hay teoría analítica que sea construida para el planeta Plutón. Nuestra mayores representaciones del movimiento con la exactitud esperada para este planeta es una integración numérica. Luego de la integración se obtienen las fórmulas presentadas y son válidas desde el año 1800 hasta el 2100. Los cálculos de los ángulos S y P son los siguientes:

S = 50,03 + 0,033459652 * d

P = 238,95 + 0,003968789 * d

lon = 238,9508 + 0,00400703 * d
- 19,799 * sin(P) + 19,848 * cos(P)
+ 0,897 * sin(2*P) - 4,956 * cos(2*P)
+ 0,610 * sin(3*P) + 1,211 * cos(3*P)
- 0,341 * sin(4*P) - 0,190 * cos(4*P)
+ 0,128 * sin(5*P) - 0,034 * cos(5*P)
- 0,038 * sin(6*P) + 0,031 * cos(6*P)
+ 0,020 * sin(S-P) - 0,010 * cos(S-P)

lat = -3,9082
- 5,453 * sin(P) - 14,975 * cos(P)
+ 3,527 * sin(2*P) + 1,673 * cos(2*P)
- 1,051 * sin(3*P) + 0,328 * cos(3*P)
+ 0,179 * sin(4*P) - 0,292 * cos(4*P)
+ 0,019 * sin(5*P) + 0,100 * cos(5*P)
- 0,031 * sin(6*P) - 0,026 * cos(6*P)
+ 0,011 * cos(S-P)

r = 40,72
+ 6,68 * sin(P) + 6,90 * cos(P)
- 1,18 * sin(2*P) - 0,03 * cos(2*P)
+ 0,15 * sin(3*P) - 0,14 * cos(3*P)

Ahora que tenemos la lon, lat y r heliocéntricas, las convertimos a coordenadas geocéntricas de la misma forma que con los otros planetas. Índice


Perturbaciones con la Luna

Primero hay que realizar unos cálculos para las correcciones.



Ms, Ml


Anomalía media del Sol y la Luna

Nl


Longitud del nodo de la Luna

ws, wl


Argumento del perihelio del Sol y la Luna

Ls = Ms + ws


Longitud media del Sol

Ll = Ml + wl + Nl


Longitud media de la Luna

D = Ll + Ls


Elongación media de la Luna

F = Ll - Nl


Argumento de la latitud de la Luna

Agregar estos valores a la longitud de la Luna (en grados):

-1,274 * sin(Ml - 2*D)
+0,658 * sin(2*D)
-0,186 * sin(Ms)
-0,059 * sin(2*Ml - 2*D)
-0,057 * sin(Ml - 2*D + Ms)
+0,053 * sin(Ml + 2*D)
+0,046 * sin(2*D - Ms)
+0,041 * sin(Ml - Ms)
-0,035 * sin(D)
-0,031 * sin(Ml + Ms)
-0,015 * sin(2*F - 2*D)
+0,011 * sin(Ml - 4*D)

Agregar estos valores a la latitud de la Luna (en grados):

-0,173 * sin(F - 2*D)
-0,055 * sin(Ml - F - 2*D)
-0,046 * sin(Ml + F - 2*D)
+0,033 * sin(F + 2*D)
+0,017 * sin(2*Ml + F)

Agregar estos valores a la distancia de la Luna:

-0,58 * cos(Ml - 2*D)
-0,46 * cos(2*D)

Entonces a la lonh, lath y rh hay que agregar los valores antes descritos respectivamente y volver a transformar el sistema esférico a rectangular, es decir xh, yh, zh (perturbado).

xh = rh * cos(lonh) * cos(lath)

yh = rh * sin(lonh) * cos(lath)

zh = rh * sin(lath)
Índice
Posición Topocéntrica de la Luna

Al calcular las coordenadas geocéntricas de la Luna, es como si un observador imaginario estuviera observando desde el centro de la Tierra. El observador real está sobre la superficie de la Tierra. Esta posición puede diferir por más de un grado desde la posición geocéntrica. Para calcular la posición topocéntrica, nosotros debemos agregar una corrección a la posición geocéntrica.

Comenzamos por calcular el paralaje de la Luna, que es el tamaño aparente del radio ecuatorial de la Tierra, vista desde la Luna:

mpar = asin( 1/r )

Donde r es la distancia de la Luna a la Tierra en radios terrestre. Esto es simplemente aplicar una corrección en coordenadas horizontales (azimut y altitud) con nuestra exactitud esperada de 1 minuto, no se necesita corrección para el azimut. Sólo aplicaremos la corrección a la altitud por encima del horizonte.

alt_topoc = alt_geoc - mpar * cos(alt_geoc)

Algunas veces es necesario corregir para la posición topocéntrica directamente en las coordenadas ecuatoriales. Por ejemplo, si uno quiere dibujar un mapa de estrellas cuando la Luna pasa en frente de las Pleiades, como será visto desde un especifico punto de la Tierra. Entonces nosotros necesitamos conocer la RA, Decl, LST y la latitud de la Luna.

Nuestra latitud astronómica (lat) primero debe ser convertida a latitud geocéntrica (gclat) y la distancia desde el centro de la Tierra (rho) a distancia de radios terrestres. Podemos hacer una simplificación, diciendo que la Tierra es perfectamente esférica, entonces:

gclat = lat, rho = 1.0

Sin embargo, si queremos más exactitud entonces:

gclat = lat - 0,1924º * sin(2*lat)

rho = 0,99833 + 0,00167 * cos(2*lat)

Ahora nosotros calculamos el Angulo Horario (HA) de la Luna:

HA = LST - RA

LST = GMST0 + UT + LON/15

GMST0 = ( Ls + 180_deg ) / 15

Necesitamos un ángulo auxiliar g:

g = atan( tan(gclat) / cos(HA) )

Estamos listos para convertir la RA y Decl en valores topográficos (topRA, topDecl):

topRA = RA - mpar * rho * cos(gclat) * sin(HA) / cos(Decl)

topDecl = Decl - mpar * rho * sin(gclat) * sin(g - Decl) / sin(g)

Note que si la decl es exactamente 90º, cos(Decl) llega ser cero y nosotros tendríamos una división por cero cuando calculamos topRA, pero la fórmula también deja de servir cuando está muy cercano al polo. También gclat es precisamente cero, q llega a ser cero también y nosotros tendríamos una división por cero cuando calculamos topDecl. en este caso, remplace la fórmula de topDecl con:

topDecl = Decl - mpar * rho * sin(-Decl) * cos(HA)

El cálculo es valido para gclat igual a cero, este también puede ser usado para gclat extremadamente cercano a cero.

Esta corrección de la posición topocéntrica pueden ser aplicado al Sol y los planetas. Pero como ellos están mucho más lejos, las correcciones son muy pequeñas, menores que nuestra exactitud esperada de un minuto de arco.

Si tu quieres calcular las coordenadas topocéntricas de los planetas, se puede hacer de la misma forma que la Luna, con una sola excepción: el paralaje de los planetas (ppar) es calculada con esta fórmula:

ppar = (8,794º/3600)/ r

Donde r es la distancia del planeta desde la Tierra en unidades astronómicas. Índice
Cálculo de la posición de los Asteroides

Nosotros debemos de contar con los elementos orbitales del asteroide para la época a considerar su posición, donde N, i, w son validos para una época especifica (normalmente usamos la época 2000). Para nuestras simplificaciones en los cómputos, solamente cambios significativos ocurren en el elemento orbital N. Para convertir el N_época al N (actual) usaremos la siguiente fórmula que simplemente agrega una corrección para la Precesión.

N = N_época +0,013967 + (2000 - época) + 3,82394e-5 * d

Donde época es expresado como un año con fracciones (por ejemplo 1950,0 o 2000,0)

Normalmente, la Anomalía media (M) se quiere para otro día del que queremos para computar la posición del asteroide. Para el movimiento diario (n) lo obtenemos simplemente agregando n por la diferencia entre días a M. Si n no se tiene, pero el período se conoce entonces n = 360/P donde P está en días. Si no se conoce el período, entonces se puede calcular con la siguiente fórmula:

P = 365,2568984 * a^1,5 (días) = 1,00004024 * a^1,5 (años)

Todos los demás elementos orbitales se calculan de igual forma como los planetas.

Índice
Cálculo de la posición de los Cometas (orbita elíptica)

Para los cometas que tienen orbitas elípticas usualmente no se tiene M en el instante T (tiempo del perihelio) si se tiene, en el afelio M es cero. Para calcular M en cualquier otro momento, primero debemos calcular d de T que llamaremos dT. Entonces la Anomalía media se calcula como:

M = 360,0 * (d - dT) / P (en grados)

Donde P es en días. Entonces a puede ser calculado teniendo q y e: a = q / (1 - e ). Luego, todos los demás elementos orbitales se calculan de igual forma como los planetas. Índice


Cálculo de la posición de los Cometas (orbita parabólica)

Si el cometa tiene una orbita parabólica un método diferente tiene que ser usado. Cuando el período orbital del cometa es infinito, M es siempre cero. La excentricidad (e) siempre es igual a 1. Como a es infinito, nosotros debemos usar la distancia del perihelio (q) para calcular la orbita parabólica con la siguiente fórmula:

Siendo la constante gravitacional Gaussiana k = 0,0172020985

H = (d - dT) * (k / raíz(2)) / q^1,5

Luego se aplican las siguientes fórmulas:

h = 1,5 * H

g = raíz(1 + h^2)

s = (g + h)^1/3 - (g - h)^1/3

A continuación se obtienen la anomalía verdadera y la distancia heliocéntrica de las fórmulas:

v = 2 * atan(s)

r = q * (1 + s^2)

Entonces, todos los demás elementos orbitales se calculan de igual forma como los planetas. Índice


Cálculo de la posición de los Cometas (orbita cercano a parabólica)

Los más comunes casos de cometas descubiertos recientemente tiene orbitas que no son exactamente una parábola, pero se acercan mucho. Su excentricidad es ligeramente superior o inferior a 1. El algoritmo que se presenta, puede ser usado para excentricidades entre 0,98 y 1,02. Si es menor que 0,98, entonces aplicamos la de orbita elíptica. No hay cometa conocido hasta ahora que tenga una excentricidad mayor a 1,02.

De igual manera tenemos que calcular d, dT y q. Entonces procedemos a calcular lo siguiente:

a = 0,75 * (d - dT) * k * raíz((1 + e) / (q^3))

b = raíz(1 + a^2)

W = (b + a)^1/3 - (b - a)^1/3

c = 1 + 1/(W^2)

f = (1 - e) / (1 + e)

g = f / (c^2)

a1 = 2/3 + 2/5 * W^2

a2 = 7/5 + 33/55 * W^2 + 37/175 * W^4

a3 = 432/175 * W^2 + 956/1125 * W^2 + 84/1575 * W^4

w = W * (1 + g * c * (a1 + a2 *g + a3 * g^2))

v = 2 * atan(w)

r = q * (1 + w^2) / (1 + w^2 * f)

Nota: este algoritmo puede fallar por mucho, desde el perihelio; sin embargo, la exactitud es suficiente para todos los cometas cercanos a Plutón. Índice

Amanecer y ocultación de cuerpos celestes.

Efemérides de los Planetas.

*************************** EJEMPLO CÁLCULOS ASTRONÓMICOS********************************

Introducción

Para nuestro ejemplo numérico, vamos a utilizar la fecha 30 de mayo del año 2002 a las 11:45 am en la ciudad de Caracas, Venezuela con una longitud de -66º 55' y una latitud de 10º 30' existiendo una diferencia horaria con el meridiano de Greenwich de 4 horas. También usaremos como herramienta de cálculo una hoja de Excel y las funciones las expresaremos en Visual Basic.

Todos los cálculos presentados en esta página son comunes para todos los cómputos de las posiciones tanto para el Sol, los Planetas, Asteroides y Cometas.

Las fórmulas presentadas en esta página, son comunes para todos los cálculos.
Escala del tiempo

UT = tiempo Caracas (horas + decimales) - (-4)

Hora = 11:45 am = 11,75 horas UT = 11,75 + 4 = 15,75

d= 367*y - int((7 * (y +(int((m + 9) / 12)))) / 4) + int(275*m / 9) + D - 730530 +UT / 24

d = 367*(2002) - int((7 * (2002 + (int((5 + 9) /12)))) / 4) + int(275 * 5 / 9) + 30 - 730530 + 15,75/24

d = 734734 - int((7 * (2002 + (int(1,166667))) / 4) + int(152,777778) - 730500 + 0,65625

d = 4234,65625 - int((7 * (2002 + 1)) / 4) + 152

d = 4386,65625 - int(7 * 2003 / 4)

d = 4386,65625 - int(3505,25)

d = 4386,65625 - 3505

d = 881,65625 días

ecl = 23,4393 - 0,0000003563*d

ecl = 23,4393 - 0,0000003563 * 881,65625 = 23,43898587 º

Hay que usar la función ang(x) = x - int( x / 360 ) * 360 para los elementos N, w y M

Código en Visual Basic para la función ang(x):

Function ang(x)
ang = x - Int(x / 360) * 360
End Function

Ahora hay que calcular la anomalía excéntrica (E)

E = M + e * (180/pi) * sin(M) * (1 + e * cos(M))

E0 = E

E1 = E0 - (E0 - e * (180/pi) * sin(E0) - M) / (1 - e * cos(E0))

Código en Visual Basic para la función E(M; e)

Private Function ecc(M, emin)
Pi = 3.14159265358979
Emay = M + emin * 180 / Pi * Sin(M * Pi / 180) * (1 + emin * Cos(M * Pi / 180))
If emin > 0.04 Then
Emay1 = Emay - (Emay - emin * 180 / Pi * Sin(Emay * Pi / 180) - M) / (1 - emin * Cos(Emay * Pi / 180))
While Abs(Emay - Emay1) > 0.001
Emay = Emay1
Emay1 = Emay - (Emay - emin * 180 / Pi * Sin(Emay * Pi / 180) - M) / (1 - emin * Cos(Emay * Pi / 180))
Wend
ecc = Emay1
Else
ecc = Emay
End If
End Function

Las funciones trigonométricas de Excel seno(), cos(), tan() y en general todas las funciones trigonométricas tiene que estar en radianes, por lo que se hace necesario, pasarlas a radianes al tener que calcularlas, ya que todos los ángulos expresados en nuestras fórmulas están en grados. Entonces se multiplica por pi y se divide entre 180. Ejemplo, cos(x * pi /180), donde x está en grados.


Tiempo Sideral y Ángulo Horario. Azimut y Altitud

Tiempo Sideral (SIDTIME), Tiempo Sideral en el Meridiano de Greenwich (GMST) y el Tiempo Universal (UT), todos en horas + decimales, son necesarios calcularlos para obtener la altitud y azimut en nuestro sitio de observación. La Altitud y Azimut, cambian a medida que transcurre el tiempo y depende del sitio geográfico del observador.

SIDTIME = GMST + UT + Longitud / 15

GMST = ang(L + 180) / 15 = ang( ang(427,9895482) +180) / 15 = ang(67,98954819 + 180) / 15

GMST = ang(247,9895482) / 15 = 247,9895482 / 15 = 16,53263655

UT = Hora local - Diferencia Horaria = 11,75 - (-4) = 15,75

SIDTIME = 16,53263655 + 15,75 + (-66,9166666667) / 15 = 27,8215254389

Si SIDTIME es negativo, sumamos 24h. Si es mayor a 24h, entonces restamos 24h.

SIDTIME = 27,8459698833 - 24 = 3,8215254389

Ángulo Horario (HA) incrementa con el tiempo (a menos que se mueva más rápido que la rotación de la Tierra, como por ejemplo, satélites). SIDTIME y RA tienen que estar en la misma unidad, por tanto hay que pasar SIDTIME a grados, multiplicándolo por 15.

HA = SIDTIME - RA

Entonces pasamos estos datos a coordenadas rectangulares, donde el eje X apunta al ecuador celeste, el eje Y apunta hacia el horizonte oeste y el eje Z apunta al polo norte celeste.

x = cos(HA) * cos(Decl)

y = sin(HA) * cos(Decl)

z = sin(Decl)

Rotamos este sistema de coordenadas

xhor = x * sin(lat) - z * cos(lat)

yhor = y

zhor = x * cos(lat) + z * sin(lat)

Finalmente pasamos del sistema de coordenadas rectangulares al sistema esférico donde el radio puede ser la unidad.

Azimut = atan2(xhor; yhor) + 180º

Altitud = asin(zhor)

¿Por qué se agrega 180º al azimut? Para adaptarlo a la forma más común de representarlo, que es, 0º el norte, 90º el este, el sur 180º, el oeste 270º y volviendo al norte. Sino se agrega 180º comienza a contar desde el sur con 0º.