Sesgo del semieje mayor medio de Brouwer-Lyddane

Ampliaré mi comentario :

No estoy familiarizado con el software, pero si inicia el satélite en el mismo todo excepto el tiempo (diferencia de 12 horas), entonces el satélite comenzará con diferentes potenciales gravitacionales para los dos casos. Supongo que si lo ejecutó 100 veces espaciando sus épocas iniciales de manera uniforme durante 24 horas, vería que sus "sesgos de altitud" se distribuirían por encima y por debajo, dependiendo de si comenzó en un pico o valle de gravedad. No debe esperarse que la media sea igual al valor inicial instantáneo, está comenzando al azar.

Para mantenerlo simple, usaré solo el$J_2$componente multipolar del potencial gravitacional de la Tierra además del término monopolo$GM_E$. Escribí un breve script en python usando números y ecuaciones de Wikipedia (los enlaces se muestran dentro del script).

He construido una órbita de inclinación de 90 grados, a una altitud de unos 901 kilómetros y propagada durante 14 órbitas o alrededor de 1 día. (Elegí esta altitud para poder establecer la inclinación en 99 grados y confirmar que era sincronizada con el sol al propagar durante 91 días y observar cómo el nodo ascendente se movía del eje x al eje y, comprobando así que el gradiente de la$J_2$término potencial está bien.)

Calculé un eje semi-mayor y la velocidad orbital.$a_0$y$v_0$para una órbita circular en el campo monopolar$GM_E/r$y los usó para las condiciones iniciales.

La primera trama muestra$x(t), y(t), z(t)$y$r(t)-a_0$durante 24 horas (14 órbitas) comenzando en el ecuador con phase = 0. El cálculo se repitió por phase = 30, 60, 90grados, de forma que el último partía por encima del polo, también a la misma altura de 901 kilómetros.

En la segunda parcela, muestro las cuatro parcelas de$r(t)-a_0$que comienzan a la misma altitud, pero en diferentes lugares sobre la Tierra. Inmediatamente puede ver que el "sesgo" del que habla depende de la ubicación inicial y el potencial gravitatorio particular que existe allí.

Comenzando por el ecuador con phase = 0grados, donde uno se sienta más bajo en el pozo de gravedad de la Tierra, la velocidad es un poco baja para una órbita circular, por lo que el satélite está realmente en apoapsis, cayendo a una altitud ligeramente más baja medio período más tarde.

Sin embargo, comenzando sobre el polo con phase = 90grados, más alto en el pozo de gravedad de la Tierra, la velocidad es un poco alta para la órbita circular, por lo que ahora es periapsis, y la altitud aumenta medio período después.

resumen: Al elegir un campo gravitatorio más fácil de entender, el comportamiento del "sesgo" es aquí más sistemático. En su pregunta, está utilizando un campo de gravedad muy complejo y de aspecto irregular que gira constantemente. Al comenzar en el mismo lugar en el espacio, pero en diferentes épocas, en realidad estabas rotando la tierra, moviendo diferentes potenciales gravitacionales a la ubicación inicial de tu satélite .

Ahora predigo que si apaga el campo multipolar en JGM-2 y usa un potencial esféricamente simétrico, sus cálculos mostrarán una altitud plana, sin variaciones ni sesgos.

Pensamientos finales: Las órbitas realistas no son cónicas perfectas, por lo que ni ellas ni sus elementos keplerianos representan órbitas realistas . Son solo aproximaciones a la realidad, por lo que no son correctas aunque estén cerca.

Los elementos keplerianos se usaban cuando las personas escribían con plumas utilizando la luz de la quema de grasa animal (si no estaban ocupados siendo quemados en la hoguera). Son una bendición mixta en el siglo XXI cuando todo tiene muchos más dígitos.

def deriv(X, t):

    rr, vel = X.reshape(2, -1)

    acc0, acc2 = accs(rr)

    return np.hstack([vel, acc0 + acc2])


def accs(rr):

    x,   y,   z   = rr
    xsq, ysq, zsq = rr**2

    rsq = (rr**2).sum()
    rm3 = rsq**-1.5
    rm7 = rsq**-3.5

    acc0  = -GM_earth * rr * rm3

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

    acc2  = J2_earth * np.hstack((acc2x, acc2y, acc2z))

    return acc0, acc2


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

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

GM_earth = 3.986004418E+14 #  m^3/s^2  https://en.wikipedia.org/wiki/Standard_gravitational_parameter
J2_earth = 1.7555E+25      #  m^5/s^2  https://en.wikipedia.org/wiki/Geopotential_model

a0  = (901 + 6378) * 1E+03 # meters (roughly) for a n=14 sun-synchronous orbit

v0  = np.sqrt(GM_earth/a0)  # m/s  (vis-viva)
T   = twopi * a0 / v0       # sec

print T/60, ' minutes'

phases = [rads*d for d in [0, 30, 60, 90]]

X0s = []
for phase in phases:
    sp, cp = np.sin(phase), np.cos(phase)

    X0 = np.hstack(( [ a0*cp, 0, a0*sp],   # initial positions
                     [-v0*sp, 0, v0*cp] )) # initial velocitys
    X0s.append(X0)


n, days = 14, 1.  # orbits/day, days

time = np.linspace(0, n*days*T, int(n*days*100)+1)  

answers = []
for X0 in X0s:

    answer, info = ODEint(deriv, X0, time, full_output=True)

    answers.append(answer)

if 1 == 1:

    names  = 'x(t)', 'y(t)', 'z(t)', 'r(t)-a0'
    answer = answers[1]   # just choose one for the plot
    posn   = answer.T[:3]
    rdiff  = np.sqrt((posn**2).sum(axis=0)) - a0
    things = [thing for thing in posn] + [rdiff]

    plt.figure()
    for i, (thing, name) in enumerate(zip(things, names)):
        plt.subplot(4, 1, i+1)
        plt.plot(time/(24.*3600), thing)   # x-axis in days, not seconds
        plt.title(name, fontsize=16)
        plt.xlim(0, days+0.01)  # a little more than n days
    plt.show()

if 1 == 1:

    plt.figure()

    for answer in answers:

        posn   = answer.T[:3]
        rdif   = np.sqrt((posn**2).sum(axis=0)) - a0

        plt.plot(time/(24.*3600), rdif)   # x-axis in days, not seconds

    plt.xlim(0, days+0.1)  # a little more than n days
    plt.text(1.01, -6000, u' 0\u00B0', fontsize=15 )
    plt.text(1.01,  2000, u'30\u00B0', fontsize=15 )
    plt.text(1.01, 16000, u'60\u00B0', fontsize=15 )
    plt.text(1.01, 22000, u'90\u00B0', fontsize=15 )
    plt.title('r(t)-a0', fontsize=16)
    plt.show()

a continuación: una verificación rápida durante 91 días con una inclinación de 99 grados para confirmar que la órbita es más o menos sincronizada con el sol para asegurarme de que al menos no haya errores graves con la forma en que escribí las ecuaciones para el gradiente de la$J_2$componente, o los valores numéricos. La secuencia de comandos anterior no permite la inclinación, pero la incluiría mezclando el$v_y$y$v_z$valores cuando X0se calcula.