Datos para "comprobar" la primera ley de Kepler

Quiero "comprobar" la primera ley de Kepler utilizando datos reales de Marte. De la ecuación de la elipse, derivé

1 r = a b 2 + a b 2 ϵ porque ( φ ) ,

dónde a es el semieje mayor, b es el semieje menor y ϵ es la excentricidad de la órbita elíptica. Estoy buscando el siguiente tipo de datos:

  1. La distancia de Marte al Sol r
  2. el ángulo φ entre Marte, el Sol y el eje principal de la órbita elíptica.

Entonces, quiero comprobar, si r y φ ajustar los valores medidos de a , b y ϵ . Si no hay tales datos disponibles (vista perpendicular en el plano orbital de Marte), ¿cómo puedo transformar los datos dados en otros sistemas de coordenadas a los que necesito? En un sitio web de la NASA ( https://omniweb.gsfc.nasa.gov/coho/helios/heli.html ) encontré datos en las coordenadas "Solar Ecliptic", "Heliographic" y "Heliographic Inertial", pero no sé que se acercan más a mi plan.

Actualizar:

Lo probé con las recomendaciones de uhoh. Lamentablemente fracasé.

Con el siguiente código de python, usando los datos de Horizons x, y, z almacenados en un archivo xlsx,

from __future__ import division
import numpy as np
from statsmodels.regression.linear_model import OLS
from statsmodels.tools import add_constant
from statsmodels.tools.eval_measures import aicc
import pandas as pd
import matplotlib.pyplot as plt
horizons = pd.read_excel("horizons2.xlsx")
horizons = np.array(horizons)
horizonsxyz=horizons[:,2:5]
horizonsxyz=np.array(horizonsxyz, dtype=np.float64)
hx=horizonsxyz[:,0]
hy=horizonsxyz[:,1]
hz=horizonsxyz[:,2]

horizonsr=np.sqrt(hx**2+hy**2+hz**2)
horizonsr=horizonsr*6.68459*(10**(-9))

phi=np.arctan2(hy, hx) * 180 / np.pi
phi2=np.mod(phi+360, 360)
phia=np.mod(phi-286, 360)
phiganz=add_constant(phia)
horizonsdurchr=1/horizonsr




horizons_regr=OLS(horizonsdurchr, phiganz).fit()
print(horizons_regr.params)
print(horizons_regr.summary())
y_pred_horizons=np.dot(phiganz, horizons_regr.params)
print(horizons_regr.params)

obtengo un valor de 7.1349 10 1 para a b 2 . Esto es malo, pero al menos en el orden correcto de magnitud. Sin embargo para a b 2 ϵ Obtengo un valor realmente malo de 2.89228 10 4 . Dividir los dos resultados produce una excentricidad estimada de 0.00044 que está muy lejos de la verdad 0.0934 .

También probé otro enfoque, utilizando los datos heliográficos mencionados anteriormente. Aquí, me acerco, pero solo si sumo 35 grados a los ángulos, lo cual no tiene sentido, ya que debo sumar 74 grados o restar 278 grados, para obtener el ángulo relativo al perihelio.

Respuestas (2)

¡Gran proyecto! y bienvenido a Stack Exchange. Publicaré una respuesta breve, pero creo que alguien puede agregar una respuesta más detallada, completa y perspicaz.

Creo que ese sitio web no se adapta bien, así que responderé basándome en que cambie a Horizons. Si te gusta Python, entonces es más divertido usar Skyfield .

Si desea aplicar una ecuación basada en un modelo de órbita de Kepler, deberá usar datos donde el Sol permanece en un lugar y Marte orbita alrededor de él. Eso sería heliocéntrico con el Sol en (0, 0, 0).

Que haya tres ceros plantea la cuestión del número de dimensiones; Las órbitas de Kepler adecuadas son una especie de 3D, es decir, tienen un plano orbital que se puede inclinar a un plano de referencia, pero las órbitas son planas. Dos problemas; su ecuación asume una órbita plana 2D debido a la forma φ se define. Idealmente, le gustaría obtener datos en el plano de la órbita de Marte y es posible que necesite transformar los datos de NASA/JPL Horizons en el plano orbital de Marte porque solo hay dos planos "oficiales" principales, ningún planeta real permanece perfectamente en un plano.

Entonces, lo que hagas depende de qué tan lejos en la madriguera del conejo de las órbitas simuladas estén los aviones que quieras ir.

Aproximación de orden cero

Ir a Horizontes

Use este tutorial y configúrelo para que coincida con lo siguiente:

Current Settings
Ephemeris Type:      VECTORS
Target Body:         Mars [499]
Coordinate Origin:   Sun (body center) [500@10]
Time Span:           Start=2020-10-04, Stop=2020-10-05, Step=1 d
Table Settings:      quantities code=2; output units=KM-S; CSV format=YES
Display/Output:      default (formatted HTML)
             -- OR --
Display/Output:      download/save (plain text file)

Configuración de JPL Horizontes

Aquí hay una línea de muestra para Marte para hoy usando el Sol como origen (he truncado algunos dígitos decimales). Ves de inmediato que Marte está a unos 201 millones de km del Sol, también está a unos 4 millones de km por debajo de la eclíptica J2000.0.

2459126.500, A.D. 2020-Oct-04 00:00:00.00,  2.036231544E+08,  5.355405115E+07, -3.872888712E+06...

Desde aquí se puede aproximar

r = X 2 + y 2 + z 2

y

φ = arcán 2 ( y , X ) 286.502°

Ya que está pasando por los cuatro cuadrantes, es mejor usar una computadora arctan2(y, x)o atan2(y, x)con dos argumentos, no arcán ( y / X ) que solo funciona en dos cuadrantes (es decir, 1/7 = -1/-7).

Aproximación de primer orden

Ves de inmediato que Marte está a unos 201 millones de km del Sol, también está a unos 4 millones de km por debajo de la eclíptica J2000.0.

Si desea corregir la inclinación de la órbita de Marte con respecto a la eclíptica, puede encontrar el mejor plano que se ajuste a los datos de un año marciano y crear su propia eclíptica de Marte.

Pero te recomiendo que primero hagas el orden cero y veas qué tan bien o mal funciona, luego puedes decidir si quieres inclinarte.

Eso es realmente útil. Sin embargo, no entiendo de dónde viene el grado 286.502.
Marte de @Joe_base Wikipedia da ese ángulo como el argumento del perihelio El X , y , z Los datos en Horizons están en un sistema de coordenadas estándar llamado J2000.0 (los detalles son bastante detallados, pero digamos que en un día determinado del año 2000, la eclíptica, el eje de la Tierra y el ecuador se "congelaron" y eso es, en términos generales, el sistema de coordenadas todos usan La dirección +x es como 0 grados, la dirección +y es como 90 grados.
Entonces, si obtienes un ángulo basado en X , y necesitas restar 286.502 para obtener el ángulo φ con respecto al perihelio.
@Joe_base está bien, ya veo. Por lo general, es malo cambiar una pregunta una vez que se publican las respuestas, pero este es un caso especial y básicamente estás refinando la misma pregunta. No tengo mucho tiempo en los próximos días para ver esto, así que agregaré una recompensa para llamar la atención y tal vez alguien más pueda intervenir. ¡Gracias!
ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf , tal vez necesite complementar con este documento
@AdrianR Esa es una gran idea, ¡genial! Los había visto hace años y los había olvidado, pero esa es una solución mucho más fácil. No tengo tiempo para ir más allá en esto, así que agregué la recompensa que vence el 20/10/2020 a las 17:00 UTC, además tiene un período de gracia de 24 horas. Si escribe una respuesta breve basada en eso, sería útil para el OP y competitivo para la recompensa, tal vez incluso agregue una segunda.
Solo para mantenerlos actualizados: finalmente descubrí que tenía que tener en cuenta la longitud del nodo ascendente. Por lo tanto, en el sistema de coordenadas J2000 (argumento del perihelio + longitud del nodo ascendente) debe deducirse para obtener los valores correctos. Doind que funcionó perfectamente bien. Al hacer eso, funcionan tanto las coordenadas heliográficas con las que comencé como la cuenta xyz de uhoh. No perseguí más los elementos keplerianos.

La primera ley de Kepler es que un planeta se mueve en una elipse con el sol en un foco. Tu ecuación es la de una elipse alrededor del foco, así que has probado la primera ley de Kepler. El φ es lo que los astrónomos llaman anomalía verdadera. Para poner tu ecuación en la forma usual, a / b 2 es 1 / pag entonces

r = pag 1 + ϵ C o s φ

Con esta ecuación, la elipse se puede trazar eligiendo muchos valores del ángulo φ y encontrar los valores r correspondientes, luego graficar.

La p es llamada parámetro por los astrónomos y semi-latus rectum por los matemáticos. Como puedes ver, cuando φ es de 90 grados, el valor de r es p. También, pag = a ( 1 ϵ 2 ) que se puede poner en la ecuación anterior como una forma alternativa de la ecuación.

La ley de Kepler no da información sobre dónde está el perihelio en el plano de la órbita.