¿Cómo calcular los planetas y lunas más allá de la fuerza gravitatoria de Newton?

"Pregunta: ¿Cómo calcular los planetas y las lunas más allá de la fuerza gravitacional de Newton?"

Uhoh, tu comentario invitó a más fuentes sobre esto. (Felicitaciones, por cierto, por todo el trabajo y los interesantes resultados que diste en tu propia respuesta).

¿Has visto lo que hizo Steve Moshier (SL Moshier) a principios de la década de 1990?

Proporcionó una réplica completa del modelo físico para las efemérides numéricamente integradas DE200/LE200 del JPL (entonces estándar), (usadas como base de los datos del sistema solar Astronomical Almanac durante los años 1984-2002), incluido el código fuente completo en C junto con el integrador adecuado & c), lo que también permite la extensión del rango de tiempo de 250 años publicado entonces para DE200 por JPL. La integración de Moshier fue comparada de forma independiente con la integración del JPL en la parte común de 250 años del rango de tiempo por J Chapront en el Observatorio de París, quien encontró que para los cinco planetas exteriores "las discrepancias nunca van más allá de 4. 10 ^ -7 arco -segundo, que es superabundante", y en el peor de los casos (luna), las discrepancias en la longitud nunca superaron los 0,008 durante el intervalo de tiempo de 250 años de DE200.

Para completar el modelo físico para que se corresponda con el estándar de entonces, Moshier tuvo que buscar información/datos más allá de lo que se había publicado, y reconoció los datos adicionales del JPL y de otros lugares.

Hasta donde yo sé, esta es la única integración estándar de efemérides del sistema solar para la cual se ha puesto a disposición del público una implementación completa y viable, y esto parece convertirlo en un trabajo notable e incluso históricamente significativo.

Las referencias a la integración DE200 de Moshier (realizada como 'DE118' en el marco de referencia de 1950 y luego rotada) son:

(Resumen del trabajo en): Moshier, SL (1992), "Comparación de una efeméride lunar de 7000 años con la teoría analítica", Astronomy and Astrophysics 262, 613-616: en http://adsabs.harvard.edu/abs /1992A%26A...262..613M .

Los detalles de implementación importantes, con el código fuente completo (C), no están en el documento de 1992, pero aún están disponibles (hasta este escrito en marzo de 2018) en el sitio web del autor en http://www.moshier.net , especialmente en archivos

http://www.moshier.net/de118i.zip ;

http://www.moshier.net/de118i-1.zip ;

y http://www.moshier.net/de118i-2.zip ;

con comentarios en http://www.moshier.net/ssystem.html .

(Estos archivos datan de 1993 a 2004, las modificaciones posteriores no fueron para cambiar el modelo, sino para acomodar la sintaxis para más compiladores, agregar comentarios y permitir opciones adicionales, como la introducción de más cuerpos en la integración, etc.)

La "referencia de resumen principal" para el modelo de física fue:

Newhall, XX, EM Standish y JG Williams (1983), "DE102: efemérides numéricamente integradas de la Luna y los planetas que abarcan cuarenta y cuatro siglos", Astronomy and Astrophysics 125, 150-167, en http://adsabs.harvard .edu/abs/1983A%26A...125..150N .

La matriz de rotación para cambiar el marco de referencia 1950->2000 fue de Standish, EM (1982), "Orientation of the JPL Ephemerides, DE200/LE200, to the Dynamical Equinox of J2000", Astronomy and Astrophysics 114, 297-302, en http://adsabs.harvard.edu/abs/1982A%26A...114..297S .

La verificación independiente se menciona en

Chapront, J. (1995), "Representación de efemérides planetarias por análisis de frecuencia. Aplicación a los cinco planetas exteriores". Suplemento de astronomía y astrofísica, v.109, 181-192: en http://adsabs.harvard.edu/abs/1995A%26AS..109..181C .

Agreguemos una aproximación para tener en cuenta algunos de los efectos de la Relatividad General (GR), al menos para los cuerpos que orbitan cerca del Sol masivo, y comencemos a observar$J_2$el término multipolar de orden más bajo para el potencial gravitacional de un cuerpo más allá del término monopolar$-GM/r$.

Newton:

La aceleración de un cuerpo en el campo gravitatorio de otro cuerpo de parámetro gravitatorio estándar $GM$puede ser escrito:

$$\mathbf{a_{Newton}} = -GM \frac{\mathbf{r}}{|r|^3},$$

donde$r$es el vector del cuerpo$M$al cuerpo cuya aceleración se está calculando. Recuerda que en la mecánica newtoniana la aceleración de cada cuerpo depende únicamente de la masa del otro cuerpo , aunque la fuerza depende de ambas masas, porque la primera masa se anula por$a=F/m$.

Relatividad general (aproximada):

Aunque no estoy familiarizado con GR, voy a recomendar una ecuación que parece funcionar bien y parece estar respaldada por varios enlaces. Es una corrección relativista aproximada a la gravedad newtoniana que se utiliza en simulaciones de mecánica orbital. Verá varios formularios en los siguientes enlaces, la mayoría con términos adicionales que no se muestran aquí:

• https://física.stackexchange.com/q/313146/83380
• ecuación 1 en https://www.lpi.usra.edu/books/CometsII/7009.pdf
• ecuación 27 en https://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf
• ecuación 4-26 en https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf
• también vea la discusión/consejo de David Hammen en esta respuesta .
• ecuación 3.11 en http://adsabs.harvard.edu/full/1994AJ....107.1885S (haga clic en "imprimir") vea esta respuesta
• ecuación A.8 en GF Rubilar1 & A. Eckart 2001 Periastron desplazamientos de órbitas estelares cerca del Centro Galáctico , excepto que el$\mathbf{r}/r^3$El newtoniano se incluye en lugar de tratarse por separado.
• ecuación 2 en M. Parsa et al. 2017 Investigando el Movimiento Relativista de las Estrellas Cerca del Agujero Negro Supermasivo en el Centro Galáctico , excepto que el$\mathbf{r}/r^3$El newtoniano se incluye en lugar de tratarse por separado.

Se debe agregar la siguiente aproximación al término newtoniano:

$$\mathbf{a_{GR}} = GM \frac{1}{c^2 |r|^3}\left(4 GM \frac{\mathbf{r}}{|r|} - (\mathbf{v} \cdot \mathbf{v}) \mathbf{r} + 4 (\mathbf{r} \cdot \mathbf{v}) \mathbf{v} \right),$$

Oblato ($J_2$solo):

Solo estoy usando las matemáticas del artículo de Wikipedia sobre el modelo geopotencial con una aproximación muy importante para recordar; Estoy suponiendo que el achatamiento está en el plano de la eclíptica, que el eje de rotación del cuerpo achatado está en el$\mathbf{\hat{z}}$dirección, perpendicular a la eclíptica. ¡No olvides que esto es una aproximación! Las ecuaciones vectoriales completas son más complicadas que esto, intentaré volver y actualizar esto una vez que esté seguro de que lo tengo correcto. Mientras tanto, aquí hay una aproximación:

$$\mathbf{r} = x \mathbf{\hat{x}} + y \mathbf{\hat{y}} + z \mathbf{\hat{z}}$$

$$a_x = J_2 \frac{x}{|r|^7} (6z^2 - 1.5(x^2+y^2))$$

$$a_y = J_2 \frac{y}{|r|^7} (6z^2 - 1.5(x^2+y^2))$$

$$a_z = J_2 \frac{z}{|r|^7} (3z^2 - 4.5(x^2+y^2))$$

Debe agregarse lo siguiente al término newtoniano:

$$\mathbf{a_{J2}} = a_x \mathbf{\hat{x}} + a_y \mathbf{\hat{y}} + a_z \mathbf{\hat{z}}$$

Fuerzas de marea:

Se complica aún más cuando se analizan términos que involucran distribuciones de masa no esféricas en ambos cuerpos al mismo tiempo, ya sea que sean estáticos o no. En este punto es probablemente necesario golpear los libros.

---

Aquí hay una prueba. He comparado datos descargados de Horizons de JPL . Para los planetas exteriores utilizo los datos de Horizons para el baricentro de cada planeta, lo que asegura que sea la velocidad del centro de masa del planeta más todas sus lunas. No he agregado la corrección a las masas del planeta, pero ese es un efecto mucho menor ya que solo afecta el movimiento de otros cuerpos distantes.

Para los datos de la Tierra, asegúrese de descargar el geocentro de la Tierra y la Luna por separado (no el baricentro Tierra-Luna). Para los planetas exteriores recuerda descargar los baricentros.

Me he integrado durante un año, y todo está dentro de aproximadamente un kilómetro de los datos de Horizons, excepto la Luna de la Tierra. Eso no es una sorpresa considerando todos los efectos de marea íntimos entre estos dos. Agregando la Tierra$J_2$al potencial sentido por la Luna solo ayuda parcialmente, realmente no es la forma correcta de hacerlo, especialmente considerando que el eje de la Tierra (y por lo tanto el achatamiento) está tan lejos de la eclíptica.

Así que, en general, esto es una ilustración de que cuanta más física pongas, ¡más cerca estarás de una simulación JPL realmente seria! Esta es la distancia absoluta entre las posiciones simuladas aquí y la salida de Horizons de 2017-01-01 00:00a 2018-01-01 00:00. A continuación está el script de Python que usé para generarlo.

---

Basado en la discusión de la rigidez en los comentarios a continuación, aquí hay un gráfico rápido del tamaño del paso en función del tiempo. Creo que la rigidez proviene del sistema Tierra-Luna, estas excursiones frecuentes se parecen un poco a las excursiones erróneas de la Tierra y la Luna. Creo que voy a tratar de reescalar el problema a unidades adimensionales, ahora mismo la tolerancia numérica se aplica a todas las velocidades y posiciones por igual, no es una buena idea.

def deriv_Newton_Only(X, t):
    x,  v  = X.reshape(2, -1)
    xs, vs = x.reshape(-1, 3), v.reshape(-1, 3)
    things = zip(bodies, xs, vs)

    accs, vels = [], []
    for a, xa, va in things:
        acc_a = np.zeros(3)
        for b, xb, vb in things:
            if b != a:
                r = xa - xb
                acc_a += -b.GM * r * ((r**2).sum())**-1.5
        accs.append(acc_a)
        vels.append(va)

    return np.hstack((np.hstack(vels), np.hstack(accs)))

def deriv_sunGRJ2EarthJ2(X, t):
    x,  v  = X.reshape(2, -1)
    xs, vs = x.reshape(-1, 3), v.reshape(-1, 3)
    things = zip(bodies, xs, vs)

    accs, vels = [], []
    for a, xa, va in things:
        acc_a = np.zeros(3)
        for b, xb, vb in things:
            if b != a:
                r = xa - xb
                acc_a += -b.GM * r * ((r**2).sum())**-1.5

        if a.do_SunGR and not a.name == 'Sun':

            a.flag0 = True

            r   = xa - xs[0]
            v   = va - vs[0]
            rsq = (r**2).sum()
            rm3 = rsq**-1.5
            rm1 = rsq**-0.5

            # https://physics.stackexchange.com/q/313146/83380
            # Eq.    1 in https://www.lpi.usra.edu/books/CometsII/7009.pdf
            # Eq.   27 in https://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf
            # Eq. 4-26 in https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf
            # Eq. 3.11 in http://adsabs.harvard.edu/full/1994AJ....107.1885S
            
            term_0 = Sun.GM / (clight**2) * rm3
            term_1 = 4.*Sun.GM * r * rm1
            term_2 =   -np.dot(v, v) * r
            term_3 = 4.*np.dot(r, v) * v

            accGR = term_0*(term_1 + term_2 + term_3)
            acc_a += accGR
            
        if a.do_SunJ2 and not a.name == 'Sun':

            a.flag1 = True

            r = xa - xs[0] # position relative to Sun
            x,   y,   z   = r
            xsq, ysq, zsq = r**2

            rsq = (r**2).sum()
            rm7 = rsq**-3.5

            # https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere
            accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))

            accJ2  = J2s * np.hstack((accJ2x, accJ2y, accJ2z))
            acc_a += accJ2
            
        if a.do_EarthJ2 and not a.name == 'Earth':

            a.flag2 = True

            r = xa - xs[3] # position relative to Earth
            
            x,   y,   z   = r
            xsq, ysq, zsq = r**2

            rsq = (r**2).sum()
            rm7 = rsq**-3.5

            # https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere
            accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
            accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))

            accJ2  = J2e * np.hstack((accJ2x, accJ2y, accJ2z))
            acc_a += accJ2
            
        accs.append(acc_a)
        vels.append(va)

    return np.hstack((np.hstack(vels), np.hstack(accs)))

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

names = ['Sun', 'Mercury', 'Venus',
         'Earth', 'Moon', 'Mars',
         'Ceres', 'Pallas', 'Vesta',
         'Jupiter', 'Saturn', 'Uranus',
         'Neptune']

GMsDE430 = [1.32712440040944E+20, 2.203178E+13,  3.24858592E+14,
        3.98600435436E+14,    4.902800066E+12,  4.2828375214E+13,
        6.28093938E+10,       1.3923011E+10,    1.7288009E+10, 
        1.267127648E+17,      3.79405852E+16,   5.7945486E+15,
        6.83652719958E+15 ] # https://ipnpr.jpl.nasa.gov/progress_report/42-178/178C.pdf

# for masses also see ftp://ssd.jpl.nasa.gov/pub/xfr/gm_Horizons.pck

# https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/
# https://naif.jpl.nasa.gov/pub/naif/JUNO/kernels/spk/de436s.bsp.lbl
# https://astronomy.stackexchange.com/questions/13488/where-can-i-find-visualize-planets-stars-moons-etc-positions
# https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/jup310.cmt
# ftp://ssd.jpl.nasa.gov/pub/xfr/gm_Horizons.pck

GMs = GMsDE430

clight = 2.9979E+08 # m/s

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]

# J2 values
J2_sun = 2.110608853272684E-07  # unitless
R_sun  = 6.96E+08 # meters
J2s    = J2_sun * (GMs[0] * R_sun**2)   # is there a minus sign?

J2_earth = 1.08262545E-03  # unitless
R_earth  = 6378136.3 # meters
J2e      = J2_earth * (GMs[3] * R_earth**2)   # is there a minus sign?

JDs, positions, velocities, linez = [], [], [], []

use_outer_barycenters = True

for name in names:

    fname = name + ' horizons_results.txt'

    if use_outer_barycenters:
        if name in ['Jupiter', 'Saturn', 'Uranus', 'Neptune']:
            fname = name + ' barycenter horizons_results.txt'

    with open(fname, 'r') as infile:

        lines = infile.read().splitlines()

        iSOE = [i for i, line in enumerate(lines) if "$$SOE" in line][0]
        iEOE = [i for i, line in enumerate(lines) if "$$EOE" in line][0]

        # print name, iSOE, iEOE, lines[iSOE], lines[iEOE]

        lines = lines[iSOE+1:iEOE]

        lines = [line.split(',') for line in lines]
        JD  = np.array([float(line[0]) for line in lines])
        pos = np.array([[float(item) for item in line[2:5]] for line in lines])
        vel = np.array([[float(item) for item in line[5:8]] for line in lines])

        JDs.append(JD)
        positions.append(pos * 1000.)   # km   to m
        velocities.append(vel * 1000.)  # km/s to m/s
        linez.append(lines)

JD = JDs[0] # assume they are identical

class Body(object):
    def __init__(self, name):
        self.name = name

bodies = []
for name, GM, pos, vel in zip(names, GMs, positions, velocities):
    
    body = Body(name)
    bodies.append(body)
    
    body.GM = GM
    
    body.daily_positions  = pos
    body.daily_velocities = vel
    
    body.initial_position = pos[0]
    body.initial_velocity = vel[0]

x0s = np.hstack([b.initial_position for b in bodies])
v0s = np.hstack([b.initial_velocity for b in bodies])

X0  = np.hstack((x0s, v0s))

ndays   = 365
nperday = 144

time = np.arange(0, ndays*24*3600+1, 24*3600./nperday)
days = time[::nperday]/(24*3600.)

for body in bodies:
    body.do_SunGR   = False
    body.do_SunJ2   = False
    body.do_EarthJ2 = False
    body.flag0      = False
    body.flag1      = False
    body.flag2      = False

Sun, Mercury, Venus, Earth, Moon, Mars = bodies[:6]
Ceres, Pallas, Vesta = bodies[6:9]
Jupiter, Saturn, Uranus, Neptune = bodies[9:]

Mercury.do_SunGR = True
Venus.do_SunGR   = True
Earth.do_SunGR   = True
Moon.do_SunGR    = True
Mars.do_SunGR    = True
Ceres.do_SunGR   = True
Pallas.do_SunGR  = True
Vesta.do_SunGR   = True

Mercury.do_SunJ2 = True

Moon.do_EarthJ2  = True

rtol = 1E-12   # this is important!!!

answer, info = ODEint(deriv_sunGRJ2EarthJ2, X0, time,
                      rtol = rtol, full_output=True)

print answer.shape

nbodies = len(bodies)

xs, vs = answer.T.reshape(2, nbodies, 3, -1)

for body, x, v in zip(bodies, xs, vs):
    body.x = x
    body.v = v
    body.x_daily = body.x[:, ::nperday]
    body.v_daily = body.v[:, ::nperday]
    body.difference = np.sqrt(((body.x_daily - body.daily_positions.T)**2).sum(axis=0))

if True:
    for body in bodies[:6]:
        print body.name, body.flag0, body.flag1, body.flag2

if True:
    plt.figure()
    for i, body in enumerate(bodies[:12]):  # Sorry Neptune!!!
        plt.subplot(4, 3, i+1)
        plt.plot(days, 0.001*body.difference)
        plt.title(body.name, fontsize=14)
        plt.xlim(0, 365)
    plt.suptitle("calc vs JPL Horizons (km vs days)", fontsize=16)
    plt.show()

Solo quiero agregar que el término de corrección relativista mencionado en la respuesta de uhoh, que es la "expansión posnewtoniana" en el nivel "1PN", aproxima los efectos relativistas al introducir un repulsivo$1/r^3$término. La expresión es utilizada por el JPL, por lo que debe usarla si desea acercarse lo más posible a sus efemérides. A pesar de que obtiene el "desplazamiento anómalo del perihelio" correctamente, obtiene efectos muy extraños de "rebote" en el límite de campo fuerte, por lo que la expresión probablemente funcione principalmente en los campos débiles de nuestro sistema solar. Realicé algunas simulaciones a continuación, el círculo verde es el radio de Schwarzschild y el círculo rojo es el radio de la "órbita circular estable más interna", ubicada a una distancia radial de tres radios de Schwarzschild. El "rebote" que se ve se debe obviamente al término repulsivo del cubo inverso. Con más momentos angulares iniciales, las órbitas se vuelven menos extrañas .