Cálculo de la posición geocéntrica de los planetas

Un problema clásico dentro de la astronomía es calcular la posición geocéntrica de los planetas en cualquier instante.
Se recomienda leer del libro Astronomía de Posición de T.J. Vives el apartado:
3.10 Coordenadas planetarias heliocéntricas y geocéntricas (pag. 161 a 178 y también el apartado 3.9 del mismo libro.
El cálculo se puede subdividir en los siguientes apartados:

Cálculo de los elementos orbitales Longitud del Nodo ascendente(W), Longitud del Perihelio LP , Oblicuidad de la Eclíptica (e) actualizados a la época T:
La causa de esta variación de la órbita se debe a la acción sobre el planeta del resto de los planetas. Si solo existiesen el Sol y el planeta el caso se reduciría a un problema de dos cuerpos y la órbita elíptica sería fija. La acción de los demás planetas se toma como una perturbación que hace variar lentamente las características y orientación de esta órbita.
Deberemos partir de los valores W0, 0, correspondientes a un instante T0 caracterizado por una fecha juliana DJ(T0) y que viene especificado en el Anuario de San Fernando.
Aquí se ha tomado T0= 10 junio de 1980 de fecha juliana DJ(T0)=2444400.5. Si representamos por k y P las variaciones respectivas en estos dos elementos correspondientes a 100 días, tendremos para los elementos eclípticos actualizados:
W=W0+k[DJ(T)-DJ(T0)]
LP=LP0+P[DJ(T0)-DJ(T0)]
Finalmente el argumento de latitud actualizado w=LP-W=-W
El cálculo de la Oblicuidad de la Eclíptica se efectúa con la expresión de Newcomb:
e =23°27'8.26"-0".46845(t-1900)-0".00000059(t-1900)2
donde t expresa el año en años trópicos. Esta es la oblicuidad media y difiere de la verdadera en la Nutación .
Cálculo de las cantidades auxiliares P (Px,Py,Pz) y Q (Qx,Qy,Qz) que vienen dados en función de w,I,e, W actualizados para la época T.
En el libro del T.J. Vives "Astronomía de Posición" página 168 se explica el valor de estos parámetros.
Hay una parte del programa que se efectúa dos veces pues hay que calcular las posiciones heliocéntricas del planeta y de la Tierra. Primero calcula las Anomalías Media, Excéntrica (para lo que tiene que resolver la ley de Kepler), y Verdadera junto con el radio vector del planeta. La segunda vez la posición geocéntrica del Sol, cuyas componentes cartesianas difieren en un signo de las de de la Tierra. Para ello tiene que volver a calcular las Anomalías anteriores. En el programa Basic que aparece al final, se incluye la resolución de la Ec. de Kepler como un subprograma.
Resolución de la Ec. de Kepler.
4a) Dado que como época de referencia se toma T0 es lógico que uno de los datos sea la A. Media ese día M0.
4b) La A. Media en el instante T valdría M=M0+n[DJ(T)-DJ(T0)] siendo n la variación diaria en la Anomalía Media.
4c) Si las épocas T y T0 están bastante separadas, la A. Media no será del primer giro para convertirlo se usa la función parte entera E(x):
M=M-360*E (M/360) con M en grados, siendo E(x) denominada INT en Basic y floor en Java. En Java se ha definido la función PrimerGiro(x).
El siguiente paso se ha hecho diferente en Java y Basic
4d) en Java
Se calcula directamente la anomalía Verdadera V para la Tierra y el planeta basándose en el desarrollo en serie de la ecuación de centro en función de la excentricidad y hasta tercer grado. Este puede ser una fuente de errores si la excentricidad del planeta es grande. El radio vector se calcula a partir de la ecuación polar de la elipse.
Por ejemplo para la Tierra MT anomalía media y ET la excentricidad.
V2=MT+2*ET*sin(MT)+1.25*ET*ET*sin(2*MT);
V2=V2+(13*sin(3*MT)-3*sin(MT))*ET*ET*ET/12;
R2=(1-ET*ET)/(1+ET*cos(V2));

4d) en Basic
Calculada la A. Excéntrica E por resolución de la ec. de Kepler, mediante iteración hasta un error dado, el cálculo de la A. Verdadera se hace mediante la expresión:

En el programa se usan letras diferentes para la excentricidad e y la A. Excéntrica E y se usa la función ARCTAN. Esto genera un problema de ambigüedad que se resuelve postulando que V/2 y E/2 pertenecen al mismo cuadrante, ello conlleva a resolver la ambigüedad mediante la condición: Si V<0 sumar a V 360°.
La distancia R se calcula mediante R=a(1-e cosE)
Se calculan las componentes cartesianas heliocéntricas de la posición solar mediante:
X=R cos (w+V)
Y=R cos e sen (w+V)
Z=R sen e sen (w+V)
y del planeta mediante:
x=Px r cos V +Qx r sen V
y=Py r cos V +Qy r sen V
z=Pz r cos V +Qz r sen V

A partir de ambas se puede calcular las componentes cartesianas Geocéntricas del planeta mediante:

x'=X+x    y'=Y+y     z'=Z+z

Ahora se trata de pasar de esas coordenadas cartesianas geocéntricas a las polares (r,a,d) distancia geocéntrica, A. recta y Declinación Geocéntrica que determinan su posición en el cielo.

Relaciones directas	Relaciones inversas

En el cálculo de la declinación no hay ambigüedad. Se puede analizar que en el calculo de la A. Recta la ambigüedad se reduce al caso con x'<0 en cuyo caso hay que sumarle a a 180° o 12 horas.
El applet de Java y el programa (guardado como EFEMER1) calcula para todos los planetas del sistema solar AR y DEC y la distancia geocéntrica.
Se colocó un bucle de 1 a 8 para calcular las posiciones de los planetas sacando de dicho bucle el cálculo de V y R y de X2 Y2 Z2 y de la oblicuidad se hace una sola vez al principio del programa.
Se elaboró una tabla con la salida de los resultados.
El applet Java y el program Basic no contemplan las interacciones mutuas entre los planetas gigantes. La mas importante es la interacción mutua entre Jupiter y Saturno debido a la resonancia existente entre ambos que hace que aproximadamente cada 2 órbitas de Saturno Júpiter de 5 giros al Sol en unos 60 años. Ello hace que la longitud heliocéntrica de ambos planetas sufra unos adelantos o atrasos. Haberlos considerado obliga a hacer unprograma específico para cada planeta. Esta interacción si se ha considerado en el cálculo de la posición de los satélites de Júpiter.
El listado del programa en Basic es:


5 REM efemerides planetas          SAVE EFEMER1

8 DEF FNFR(X)=X-INT(X)

15 RESTORE:CLS:SCREEN 0:PI=3.1415928#:WIDTH 80

20 PRINT " COORDENADAS ECUATORIALES GEOCENTRICAS "

70 READ PT,VPT,AT,ET,NT,MT

80 DATA 282.60402,0.00471,1,0.016718,0.9856,155.9044

90 READ F1,D1,M1,A1

110 DATA 2415020.5,10,6,1980

120 PRINT  "INGRESE FECHA CALCULO EFEMERIDES:   ";

130 INPUT "DIA MES AÑO HORA(h.m) ";D2,M2,A2,H

150 HR=INT(H)+FNFR(H)*5/3:D2=D2+HR/24

190 DY=D1:A=A1:M=M1

220 GOSUB 1230

230 J1=J

240 DY=D2:A=A2:M=M2

270 GOSUB 1230

271 J2=J

272 PRINT :PRINT

273 U=(J2-F1)/365.2422

274 E=23.45229444#-(.46845*U+5.9E-07*U^2)/3600

275 E=E*PI/180

283 GOSUB 2000

285 PRINT TAB(1) "PLANETA";TAB(13) "R (UA)";TAB(32) 

    "ASCENSION RECTA";TAB(55)"DECLINACION"

290 PRINT "==========================================

    ===================================="

295 FOR XX=1 TO 8

300 READ P$,AP,EP,I,ND,VND,T0,LP,VLP,MP,NP

320 ND=ND+VND*((J2-J1)/100)

330 LP=LP+VLP*((J2-J1)/100)

340 W=LP-ND

430 ND=ND*PI/180

440 W=W*PI/180

460 I=I*PI/180

470 P1=COS (ND)*COS(W)-COS(I)*SIN(ND)*SIN(W)

480 P2=SIN(ND)*COS(W)*COS(E)+COS(I)*COS(ND)*SIN(W)*COS(E)

    -SIN(I)*SIN(W)*SIN(E)

500 P3=SIN(ND)*COS(W)*SIN(E)+COS(I)*COS(ND)*SIN(W)*SIN(E)

    +SIN(I)*SIN(W)*COS(E)

520 Q1=-COS (ND)*SIN(W)-COS(I)*SIN(ND)*COS(W)

530 Q2=-SIN(ND)*SIN(W)*COS(E)+COS(I)*COS(ND)*COS(W)*COS(E)

    -SIN(I)*COS(W)*SIN(E)

540 Q3=COS(ND)*COS(W)*SIN(E)*COS(I)-SIN(ND)*SIN(W)*SIN(E)

    +SIN(I)*COS(W)*COS(E)

620 REM CALCULO DE LA EC. KEPLER PARA PLANETA

630 M=MP+NP*(J1-T0):N=NP:E1=EP:A=AP

670 GOSUB 1480

680 V1=V:R1=R

820 V1=V1*PI/180

840 X1=P1*R1*COS(V1)+Q1*R1*SIN(V1)

850 Y1=P2*R1*COS(V1)+Q2*R1*SIN(V1)

860 Z1=P3*R1*COS(V1)+Q3*R1*SIN(V1)

970 X=X1+X2

980 Y=Y1+Y2

990 Z=Z1+Z2

1020 REM calculo de AR y DEC

1030 A1=ATN(Y/X)*12/PI

1040 IF X>=0 GOTO 1060

1050 A1=A1+12

1060 IF A1>=0 GOTO 1080

1070 A1=A1+24

1080 G=INT (A1)

1090 M=INT((A1-INT(A1))*60)

1100 S=((A1-INT(A1))*60-M)*60

1110 DC=ATN(Z/SQR(X*X+Y*Y))*180/PI

1111 D=ABS(DC)

1120 G1=SGN(DC)*INT (D)

1130 M1=INT((D-INT(D))*60)

1140 S1=((D-INT(D))*60-M1)*60

1150 R=SQR(X*X+Y*Y+Z*Z)

1160 PRINT TAB(1) P$;TAB(13) R;TAB(30) G;"H ";M;"M ";INT(S);

     "S ";TAB(52) G1;"G";M1;"M";INT (S1);"S"

1180 NEXT XX

1190 LOCATE 22,4:INPUT "OTRO CALCULO ";A$

1195 IF A$="s" OR A$="S" THEN 15

1200 END

1220 REM dia juliano

1230 MN=M:YR=A

1240 IF MN=1 OR MN=2 THEN 1260

1250 GOTO 1270

1260 YR=YR-1:MN=MN+12

1270 L=INT(YR/100)

1280 P=2-L+INT(L/4)

1290 Q=INT(365.25*YR)

1300 C= INT(30.6001*(MN+1))

1310 W=P+Q+C+DY

1320 J=W+1720994.5#

1330 RETURN

1470 REM ECUACION DE KEPLER

1480 M=M+N*(J2-J1)

1490 M=M-INT(M/360)*360

1500 K=M

1510 P=M+E1*180/PI*SIN(PI*K/180)

1520 IF ABS(P-K)>=.00001 GOTO 1550

1530 K=P

1540 GOTO 1510

1550 V=180/PI*ATN(SQR((1+E1)/(1-E1))*TAN(PI*P/360))*2

1560 IF V>0 THEN 1580

1570 V=V+360

1580 R=A*(1-E1*COS (PI*P/180))

1590 RETURN

1600 DATA MERCURIO,.3870986,.20563,7.0043579,48.094173,.00325

1610 DATA 2444238.5,77.14442128,.004265,154.15309,4.092339

1620 DATA VENUS,.7233316,.0067826,3.394435,76.4997524,.00247

1630 DATA 2444238.5,131.2895792,.003847,224.44394,1.602130

1660 DATA MARTE,1.5236883,.0933865,1.8498011,49.4032001,.00211

1670 DATA 2444238.5,335.690816,.005038,150.61701,.524033

1680 DATA JUPITER,5.202561,.0484658,1.3041819,

     100.2520175,.00276

1690 DATA 2444238.5,14.0095493,.004412,132.95682,.0830912

1700 DATA SATURNO,9.554747,.0556155,2.4893741,

     113.4888341,.00239

1710 DATA 2444238.5,92.665397,.005371,72.65684,.0334596

1720 DATA URANO,19.21814,.0463232,.7729895,73.8768642,.00137

1730 DATA 2444238.5,172.7363288,.004063,55.33453,0.0117308

1740 DATA NEPTUNO,30.10957,.0090021,1.7716016,

     131.5660649,.00301

1750 DATA 2444238.5,47.8672148,.003899,212.49069,.0059819

1760 DATA PLUTON,39.78459,.25387,17.137,109.941,.00305

1770 DATA 2444240.5,222.972,.00083,346.467,.0039746

2000 REM calculo coordenadas cartesianas tierra

2010 REM CALCULO DE LA EC. KEPLER PARA TIERRA

2020 M=MT:N=NT:E1=ET:A=AT

2030 GOSUB 1480

2040 V2=V:R2=R

2050 V2=V2*PI/180

2060 PT=PT+VPT*(J2-J1)/100

2070 PT=PT*PI/180

2080 X2=R2*COS(PT+V2)

2090 Y2=R2*COS(E)*SIN(PT+V2)

2100 Z2=R2*SIN(E)*SIN(PT+V2)

2110 RETURN