El modelizado de órbitas es el proceso de crear modelos matemáticos para simular el movimiento de un cuerpo masivo a medida que se mueve en órbita alrededor de otro cuerpo masivo debido a la gravedad. Otras fuerzas como la atracción gravitatoria de terceros cuerpos, la resistencia del aire, el viento solar o el empuje de un motor cohete, se suelen incluir en el modelo como efectos de segundo orden. Modelar directamente una órbita puede llevar a rebasar los límites de precisión de cálculo de los ordenadores, debido a la necesidad de modelar perturbaciones pequeñas en órbitas muy grandes. Debido a esto, los métodos de cálculo de perturbaciones a menudo se usan para modelar órbitas con el fin de lograr una mayor precisión.

El estudio del movimiento orbital y el modelado matemático de las órbitas comenzó con los primeros intentos de predecir los movimientos planetarios en el cielo, aunque en la antigüedad las causas seguían siendo un misterio. Newton, en el momento en que formuló sus leyes del movimiento y de la gravitación, las aplicó al primer análisis de perturbaciones,^[1]​ reconociendo las complejas dificultades de su cálculo.^[1]​ Muchos de los grandes matemáticos desde entonces han prestado atención a los diversos problemas involucrados. Así, durante los siglos XVIII y XIX hubo una demanda de tablas precisas de la posición de la Luna y los planetas para fines de navegación en alta mar (véase la historia de la longitud).

Los movimientos complejos de las órbitas se pueden descomponer. El movimiento hipotético que el cuerpo sigue bajo el efecto gravitacional de otro cuerpo es típicamente una sección cónica, y puede modelarse fácilmente con métodos geométricos. A este caso se le denomina problema de los dos cuerpos o una órbita kepleriana imperturbable. Las diferencias entre una órbita de Kepler y el movimiento real del cuerpo son causadas por perturbaciones, consecuencia de fuerzas ajenas al efecto gravitatorio entre el cuerpo primario y el secundario. Estas fuerzas deben modelizarse para crear una simulación precisa de la órbita considerada. La mayoría de los enfoques de modelado de órbitas, abordan primero el problema de los dos cuerpos y luego agregan modelos de estas fuerzas perturbadoras y simulan estos modelos para lapsos de tiempo determinados. Las fuerzas perturbadoras pueden incluir la atracción gravitacional de otros cuerpos además de la primaria, el viento solar, la resistencia aerodinámica, los campos magnéticos y las fuerzas de propulsión.

Existen soluciones analíticas (expresiones matemáticas para predecir las posiciones y movimientos en cualquier momento futuro) para dos cuerpos simples y para el problema de los tres cuerpos; pero no se ha encontrado ninguna para el problema de los n cuerpos a excepción de ciertos casos especiales. Incluso el problema de dos cuerpos se vuelve irresoluble si uno de los cuerpos tiene una forma irregular.^[2]​

Debido a la dificultad de encontrar soluciones analíticas para la mayoría de los problemas de interés, generalmente se utilizan programas computarizados de modelado y simulación para analizar movimientos orbitales. Aplicaciones de software comercial como Satellite Tool Kit se han creado con el propósito específico de simular órbitas y trayectorias de naves espaciales.

En su forma más simple, se puede crear un modelo de órbita suponiendo que solo están involucrados dos cuerpos, ambos se comportan como masas puntuales esféricas, y que ninguna otra fuerza actúa sobre los cuerpos. Para este caso, el modelo se simplifica a una órbita de Kepler.

Las órbitas keplerianas siguen trayectorias con la forma de secciones cónicas. El modelo matemático de la órbita que da la distancia entre un cuerpo central y un cuerpo en órbita se puede expresar como:

${\displaystyle r(\nu )={\frac {a(1-e^{2})}{1+e\cos(\nu )}}}$

donde:

${\displaystyle r}$ es la distancia

${\displaystyle a}$ es el semieje mayor, que define el tamaño de la órbita

${\displaystyle e}$ es la excentricidad orbital, que define la forma de la órbita

${\displaystyle \nu }$ es la anomalía verdadera, que es el ángulo entre la posición actual del objeto en órbita y la ubicación en la órbita más cercana al cuerpo central (llamada ápside)

Alternativamente, la ecuación se puede expresar como:

${\displaystyle r(\nu )={\frac {p}{1+e\cos(\nu )}}}$

donde ${\displaystyle p}$ es el semiancho recto de la curva. Esta forma de la ecuación es particularmente útil cuando se trata de trayectorias parabólicas, para las que el semieje mayor es infinito.

Un enfoque alternativo utiliza la Ley de gravitación universal de Isaac Newton como se define a continuación:

${\displaystyle F=G{\frac {m_{1}m_{2}}{r^{2}}}}$

donde:

${\displaystyle F}$ es la magnitud de la fuerza gravitacional entre las dos masas puntuales

${\displaystyle G}$ es la constante de gravitación universal

${\displaystyle m_{1}}$ es la masa del primer punto de material

${\displaystyle m_{2}}$ es la masa del segundo punto de material

${\displaystyle r}$ es la distancia entre las dos masas de puntuales

Haciendo una suposición adicional de que la masa del cuerpo primario es mucho mayor que la masa del cuerpo secundario, hace que la segunda ley de Newton, dé como resultado la siguiente ecuación diferencial

${\displaystyle {\ddot {\mathbf {r} }}={\frac {Gm_{1}}{r^{2}}}\mathbf {\hat {r}} }$

La solución de esta ecuación diferencial proporciona el movimiento de Kepler para una órbita. En la práctica, las órbitas de Kepler suelen ser útiles solo para aproximaciones de primer orden, casos especiales o como modelo base para una órbita con perturbaciones.

Los modelos de órbitas generalmente consideran la propagación de las perturbaciones en el tiempo y en el espacio utilizando métodos especiales. Para ello, se modeliza primero la órbita principal como una órbita Kepleriana y a continuación se agregan las perturbaciones al modelo para tener en cuenta las diversas alteraciones que afectan a la órbita.^[1]​ Se pueden aplicar perturbaciones especiales a cualquier problema en mecánica celeste, ya que no se limita a los casos en que las fuerzas son pequeñas.^[2]​ Los métodos de perturbación especiales son la base de las efemérides planetarias^[1]​ generadas por ordenadores muy potentes y de gran precisión (véase por ejemplo la entrada del Jet Propulsion Laboratory Development Ephemeris)

El método de Cowell es quizás el más simple de los métodos especiales de cálculo de perturbaciones:^[3]​ matemáticamente, para ${\displaystyle n}$ cuerpos que interactúan mutuamente, las fuerzas newtonianas ${\displaystyle i}$ que actúan sobre el cuerpo procedentes de los otros cuerpos ${\displaystyle j}$ simplemente se suman así,

${\displaystyle \mathbf {\ddot {r}} _{i}=\sum _{\underset {j\neq i}{j=1}}^{n}{Gm_{j}(\mathbf {r} _{j}-\mathbf {r} _{i}) \over r_{ij}^{3}}}$

donde

${\displaystyle \mathbf {\ddot {r}} _{i}}$ es el vector aceleración del cuerpo ${\displaystyle i}$

${\displaystyle G}$ es la constante de gravitación universal

${\displaystyle m_{j}}$ es la masa del cuerpo ${\displaystyle j}$

${\displaystyle \mathbf {r} _{i}}$ y ${\displaystyle \mathbf {r} _{j}}$ son las posiciones de los objetos ${\displaystyle i}$ y ${\displaystyle j}$

${\displaystyle r_{ij}}$ es la distancia del objeto ${\displaystyle i}$ al objeto ${\displaystyle j}$

con todos los vectores referidos al centro de masas del sistema. Esta ecuación se resuelve en componentes en ${\displaystyle x}$, ${\displaystyle y}$, ${\displaystyle z}$ y estos se integran numéricamente para formar los nuevos vectores de velocidad y posición a medida que la simulación avanza en el tiempo. La ventaja del método de Cowell es la facilidad de aplicación y programación. Una desventaja es que cuando las perturbaciones se vuelven grandes en magnitud (como cuando un objeto se acerca mucho a otro) los errores del método también se vuelven grandes.^[4]​ Otra desventaja es que en sistemas con un cuerpo central dominante, como el Sol, es necesario calcular muchas cifras significativas en las operaciones aritméticas debido a la gran diferencia entre la fuerza del cuerpo central y las de los cuerpos perturbadores.^[5]​

El método de Encke comienza con la órbita osculante como referencia y se integra numéricamente para resolver la variación de la referencia en función del tiempo.^[6]​ Sus ventajas son que las perturbaciones generalmente son de magnitud pequeña, por lo que la integración puede avanzar en pasos más amplios (con los consiguientes menores errores), y el método se ve mucho menos afectado por perturbaciones extremas que el método de Cowell. Su desventaja es la complejidad; no puede usarse indefinidamente sin actualizar ocasionalmente la órbita de osculación y continuar desde allí, en un proceso conocido como rectificación.^[4]​ ^[7]​

Haciendo que ${\displaystyle {\boldsymbol {\rho }}}$ sea el radio vector de la órbita osculante, ${\displaystyle \mathbf {r} }$ el vector del radio de la órbita perturbada y ${\displaystyle \delta \mathbf {r} }$ la variación de la órbita de osculación, entonces

${\displaystyle \mathbf {\ddot {r}} }$ y ${\displaystyle {\boldsymbol {\ddot {\rho }}}}$ son solo las ecuaciones de movimiento de ${\displaystyle \mathbf {r} }$ y ${\displaystyle {\boldsymbol {\rho }}}$,

donde ${\displaystyle \mu =G(M+m)}$ es el parámetro gravitacional estándar con ${\displaystyle M}$ y ${\displaystyle m}$ las masas del cuerpo central y del cuerpo perturbado, ${\displaystyle \mathbf {a} _{\text{per}}}$ es la aceleración perturbadora, y ${\displaystyle r}$ y ${\displaystyle \rho }$ son las magnitudes de ${\displaystyle \mathbf {r} }$ y ${\displaystyle {\boldsymbol {\rho }}}$.

Sustituyendo las ecuaciones (3) y (4) en la ecuación (2),

que, en teoría, podría integrarse dos veces para encontrar ${\displaystyle \delta \mathbf {r} }$. Dado que la órbita de osculación se calcula fácilmente mediante métodos de dos cuerpos, se toman en cuenta ${\displaystyle {\boldsymbol {\rho }}}$ y ${\displaystyle \delta \mathbf {r} }$ y se puede resolver ${\displaystyle \mathbf {r} }$. En la práctica, la cantidad en los corchetes, ${\displaystyle {{\boldsymbol {\rho }} \over \rho ^{3}}-{\mathbf {r} \over r^{3}}}$, es la diferencia de dos vectores casi iguales, y es necesaria un número mayor de cálculos para evitar la necesidad de utilizar numerosas cifras significativas adicionales.^[8]​^[9]​

En 1991 Victor R. Bond y Michael F. Fraietta crearon un método eficiente y altamente preciso para resolver el problema de dos cuerpos perturbado.^[10]​ Este método utiliza las ecuaciones diferenciales de movimiento linealizadas y regularizadas derivadas de Hans Sperling y una teoría de perturbación basada en estas ecuaciones desarrolladas por C.A. Burdet en el año 1864. En 1973, Bond y Hanssen mejoraron el conjunto de ecuaciones diferenciales de Burdet al usar la energía total del sistema perturbado como un parámetro en lugar de la energía de dos cuerpos y al reducir el número de elementos a 13. En 1989 Bond y Gottlieb incorporaron una integral jacobiana, que es una constante cuando la función potencial depende explícitamente del tiempo y de la posición en las ecuaciones newtonianas. La constante jacobiana se usó como un elemento para reemplazar la energía total en una reformulación de las ecuaciones diferenciales del movimiento. En este proceso, se introduce otro elemento que es proporcional a una componente del momento angular. Esto hizo que la cantidad total de elementos volviera a 14. En 1991, Bond y Fraietta realizaron nuevas revisiones al reemplazar el vector de Laplace con otra integral de vector así como otra integral escalar que eliminó los términos seculares pequeños que aparecían en las ecuaciones diferenciales para algunos de los elementos.^[11]​

El método de Sperling-Burdet se ejecuta en un proceso de 5 pasos de la siguiente manera:^[11]​

Paso 1: Inicialización

Dada una posición inicial, ${\displaystyle \mathbf {r} _{0}}$, una velocidad inicial, ${\displaystyle \mathbf {v} _{0}}$ y una hora inicial, ${\displaystyle t_{0}}$, se inicializan las siguientes variables:

${\displaystyle s=0}$

${\displaystyle r_{0}=(\mathbf {r} _{0}\cdot \mathbf {r} _{0})^{1/2}}$

${\displaystyle a=r_{0}}$

${\displaystyle b=\mathbf {r} _{0}\cdot \mathbf {v} _{0}}$

${\displaystyle \tau =t_{0}}$

${\displaystyle {\boldsymbol {\alpha }}=\mathbf {r} _{0}}$

${\displaystyle {\boldsymbol {\beta }}=a\mathbf {v} _{0}}$

Perturbaciones debidas a masas perturbadoras, definidas como ${\displaystyle V_{0}}$ y ${\displaystyle {\Bigg [}{\partial {V} \over {\partial {\mathbf {r} }}}{\Bigg ]}_{0}}$, son evaluadas

Perturbaciones debido a otras aceleraciones, definidas como ${\displaystyle \mathbf {P} _{0}}$, son evaluadas

${\displaystyle \alpha _{J}={\frac {2\mu }{r_{0}}}-\mathbf {v} _{0}\cdot \mathbf {v} _{0}-2V_{0}}$

${\displaystyle \gamma =\mu -\alpha _{J}a}$

${\displaystyle {\boldsymbol {\delta }}=-(\mathbf {v} _{0}\cdot \mathbf {v} _{0})\mathbf {r} _{0}+(\mathbf {r} _{0}\cdot \mathbf {v} _{0})\mathbf {v} _{0}+{\frac {\mu }{r_{0}}}\mathbf {r} _{0}-\alpha _{J}\mathbf {r} _{0}}$

${\displaystyle \sigma =0}$

Paso 2: Transformar elementos en coordenadas

${\displaystyle \mathbf {r} ={\boldsymbol {\alpha }}+{\boldsymbol {\beta }}sc_{1}+{\boldsymbol {\delta }}s^{2}c_{2}}$

${\displaystyle \mathbf {r'} ={\boldsymbol {\beta }}c_{0}+{\boldsymbol {\delta }}sc_{1}}$

${\displaystyle \mathbf {x} _{3}=\alpha _{J}({\boldsymbol {\alpha }}-\mathbf {r} )+{\boldsymbol {\delta }}}$

${\displaystyle \gamma =\mu -\alpha _{J}a}$

${\displaystyle r=a+bsc_{1}+\gamma s^{2}c_{2}}$

${\displaystyle \mathbf {v} =\mathbf {r'} /r}$

${\displaystyle r'=bc_{0}+\gamma sc_{1}}$

${\displaystyle t=\tau +as+bs^{2}c_{2}+\gamma s^{3}c_{3}}$

donde ${\displaystyle c_{0},c_{1},c_{2},c_{3}}$ son funciones de Stumpff

Paso 3: Evaluar ecuaciones diferenciales para los elementos

${\displaystyle \mathbf {F} =\mathbf {P} -{\partial {V} \over \partial {\mathbf {r} }}}$

${\displaystyle \mathbf {Q} =r^{2}\mathbf {F} +2\mathbf {r} (-V+\sigma )}$

${\displaystyle \alpha '_{J}=2(-\mathbf {r'} +r{\boldsymbol {\omega }}\times \mathbf {r} )\cdot \mathbf {P} }$

${\displaystyle \mu {\boldsymbol {\epsilon }}'=2(\mathbf {r'} \cdot \mathbf {F} )\mathbf {r} -(\mathbf {r} \cdot \mathbf {F} )\mathbf {r'} -(\mathbf {r} \cdot \mathbf {r'} )\mathbf {F} }$

${\displaystyle {\boldsymbol {\alpha }}'=-\mathbf {Q} sc_{1}-\mu {\boldsymbol {\epsilon }}'s^{2}c_{2}-\alpha '_{J}{\big [}{\boldsymbol {\alpha }}s^{2}c_{2}+2{\boldsymbol {\beta }}s^{3}{\bar {c}}_{3}+{\frac {1}{2}}{\boldsymbol {\delta }}s^{4}c_{2}^{2}{\big ]}}$

${\displaystyle {\boldsymbol {\beta }}'=\mathbf {Q} c_{0}+\mu {\boldsymbol {\epsilon }}'sc_{1}+\alpha '_{J}{\big [}{\boldsymbol {\alpha }}sc_{1}+{\boldsymbol {\beta }}s^{2}{\bar {c}}_{2}-{\boldsymbol {\delta }}s^{3}(2{\bar {c}}_{3}-c_{1}c_{2}){\big ]}}$

${\displaystyle {\boldsymbol {\delta }}'=\mathbf {Q} \alpha _{J}sc_{1}-\mu {\boldsymbol {\epsilon }}'c_{0}+\alpha '_{J}{\big [}-{\boldsymbol {\alpha }}c_{0}+2\alpha _{J}{\boldsymbol {\beta }}s^{3}{\bar {c}}_{3}+{\frac {1}{2}}{\boldsymbol {\delta }}\alpha _{J}s^{4}c_{2}^{2}{\big ]}}$

${\displaystyle \sigma '=r{\boldsymbol {\omega }}\cdot \mathbf {r} \times \mathbf {F} }$

${\displaystyle a'=-{\frac {1}{r}}\mathbf {r} \cdot \mathbf {Q} sc_{1}-\alpha _{J}'{\big [}as^{2}c_{2}+2bs^{3}{\bar {c}}_{3}+{\frac {1}{2}}\gamma s^{4}c_{2}^{2}{\big ]}}$

${\displaystyle b'={\frac {1}{r}}\mathbf {r} \cdot \mathbf {Q} c_{0}+\alpha _{J}'{\big [}asc_{1}+bs^{2}{\bar {c}}_{2}-\gamma s^{3}(2{\bar {c}}_{3}-c_{1}c_{2}){\big ]}}$

${\displaystyle \gamma '=-{\frac {1}{r}}\mathbf {r} \cdot \mathbf {Q} \alpha _{J}sc_{1}+\alpha _{J}'{\big [}-ac_{0}+2b\alpha _{J}s^{3}{\bar {c}}_{3}+{\frac {1}{2}}\gamma \alpha _{J}s^{4}c_{2}^{2}{\big ]}}$

${\displaystyle \tau '={\frac {1}{r}}\mathbf {r} \cdot \mathbf {Q} s^{2}c_{2}+\alpha _{J}'{\big [}as^{3}c_{3}+{\frac {1}{2}}bs^{4}c_{2}^{2}-2\gamma s^{5}(c_{5}-4{\bar {c}}_{5}){\big ]}}$

Paso 4: Integración

Aquí las ecuaciones diferenciales se integran durante un período ${\displaystyle \Delta s}$ para obtener el valor del elemento en ${\displaystyle s+\Delta s}$

Paso 5: Avance

Configurando ${\displaystyle s=s+\Delta s}$ y regresando al paso 2 hasta que se cumplan las condiciones de detención de la simulación.

Las fuerzas perturbadoras hacen que las órbitas se aparten de una órbita Kepleriana perfecta. Los modelos para cada una de estas fuerzas se crean y ejecutan durante la simulación orbital para que se puedan determinar sus efectos en la órbita.

La Tierra no es una esfera perfecta ni su masa está distribuida uniformemente en su interior. Esto hace que el modelo de gravedad de punto-masa sea inexacto para las órbitas alrededor de la Tierra, particularmente en orbita baja terrestre. Para tener en cuenta las variaciones en el potencial gravitatorio alrededor de la superficie de la Tierra, su campo gravitatorio se modela con armónicos esféricos.^[12]​ que se expresan a través de la ecuación:

${\displaystyle {\mathbf {f} }=-{\frac {\mu }{R^{2}}}\mathbf {\hat {r}} +\sum _{n=2}^{\infty }\sum _{m=0}^{n}{\mathbf {f} }_{n,m}}$

donde

${\displaystyle {\mu }}$ es el parámetro gravitatorio definido como el producto de G, la constante de gravitación universal, y la masa del cuerpo primario.

${\displaystyle \mathbf {\hat {r}} }$ es el vector unitario que define la distancia entre los cuerpos primario y secundario, siendo ${\displaystyle {R}}$ la magnitud de la distancia.

${\displaystyle {\mathbf {f} }_{n,m}}$ representa la contribución a ${\displaystyle {\mathbf {f} }}$ del armónico esférico de grado n y orden m, que se define como:^[12]​

${\displaystyle {\begin{aligned}\mathbf {f} _{n,m}&={\frac {\mu R_{O}^{2}}{R^{n+m+1}}}\left({\frac {C_{n,m}{\mathcal {C}}_{m}+S_{n,m}{\mathcal {S}}_{m}}{R}}(A_{n,m+1}\mathbf {\hat {e}} _{3}-\left(s_{\lambda }A_{n,m+1}+(n+m+1)A_{n,m}\right)\mathbf {\hat {r}} \right)\\[10pt]&{}\quad {}+mA_{n,m}((C_{n,m}{\mathcal {C}}_{m-1}+S_{n,m}{\mathcal {S}}_{m-1})\mathbf {\hat {e}} _{1}+(S_{n,m}{\mathcal {C}}_{m-1}-C_{n,m}{\mathcal {S}}_{m-1})\mathbf {\hat {e}} _{2}))\end{aligned}}}$

donde:

${\displaystyle R_{O}}$ es el radio ecuatorial medio del cuerpo primario.

${\displaystyle R}$ es la magnitud del vector de posición desde el centro del cuerpo primario hasta el centro del cuerpo secundario.

${\displaystyle C_{n,m}}$ y ${\displaystyle S_{n,m}}$ son coeficientes gravitacionales de grado n y orden m. Estos se encuentran típicamente a través de medidas de gravimetría.

Los vectores unitarios ${\displaystyle \mathbf {\hat {e}} _{1},\mathbf {\hat {e}} _{2},\mathbf {\hat {e}} _{3}}$ definen un sistema de coordenadas fijado en el cuerpo primario. Para la Tierra, ${\displaystyle \mathbf {\hat {e}} _{1}}$ se encuentra en el plano ecuatorial paralelo a una línea que se cruza con el centro geométrico de la Tierra y los puntos del meridiano de Greenwich, ${\displaystyle \mathbf {\hat {e}} _{3}}$ en la dirección del eje polar norte y ${\displaystyle \mathbf {\hat {e}} _{2}=\mathbf {\hat {e}} _{3}\times \mathbf {\hat {e}} _{1}}$

${\displaystyle A_{n,m}}$ se conocen como polinomios de Legendre derivados de grado n y orden m. Se resuelven a través de una relación de recurrencia: ${\displaystyle A_{n,m}(u)={\frac {1}{n-m}}((2n-1)uA_{n-1,m}(u)-(n+m-1)A_{n-2,m}(u))}$

${\displaystyle s_{\lambda }}$ es el seno de la latitud geográfica del cuerpo secundario, que es ${\displaystyle \mathbf {\hat {r}} \cdot \mathbf {\hat {e}} _{3}}$.

${\displaystyle {\mathcal {C}}_{m},{\mathcal {S}}_{m}}$ se define con la siguiente relación de recurrencia y condiciones iniciales: ${\displaystyle {\mathcal {C}}_{m}={\mathcal {C}}_{1}{\mathcal {C}}_{m-1}-{\mathcal {S}}_{1}{\mathcal {S}}_{m-1},{\mathcal {S}}_{m}={\mathcal {S}}_{1}{\mathcal {C}}_{m-1}+{\mathcal {C}}_{1}{\mathcal {S}}_{m-1},{\mathcal {S}}_{0}=0,{\mathcal {S}}_{1}=\mathbf {R} \cdot \mathbf {\hat {e}} _{2},{\mathcal {C}}_{0}=1,{\mathcal {C}}_{1}=\mathbf {R} \cdot \mathbf {\hat {e}} _{1}}$

Cuando se modelan perturbaciones de una órbita alrededor de un cuerpo primario, solo se debe incluir la suma de los términos ${\displaystyle {\mathbf {f} }_{n,m}}$ en la perturbación, ya que el modelo de gravedad punto-masa se cuenta en el término ${\displaystyle -{\frac {\mu }{R^{2}}}\mathbf {\hat {r}} }$.

Las fuerzas gravitacionales de terceros cuerpos pueden causar perturbaciones en una órbita. Por ejemplo, el Sol y la Luna causan perturbaciones en las órbitas alrededor de la Tierra.^[13]​ Estas fuerzas se modelan de la misma manera que la gravedad se modela para el cuerpo primario mediante una simulación de n cuerpos. Normalmente, solo se utiliza un modelo de gravedad de punto esférico para modelar los efectos de estos terceros cuerpos.^[14]​ Algunos casos especiales de perturbaciones de un tercer cuerpo tienen soluciones analíticas aproximadas. Por ejemplo, las perturbaciones para la ascensión recta del nodo ascendente y el argumento del perigeo para una órbita terrestre circular son:^[13]​

${\displaystyle {\dot {\Omega }}_{\mathrm {LUNA} }=-0.00338(\cos(i))/n}$

${\displaystyle {\dot {\omega }}_{\mathrm {LUNA} }=-0.00169(4-5\sin ^{2}(i))/n}$

donde:

${\displaystyle {\dot {\Omega }}}$ es el cambio a la ascensión recta del nodo ascendente en grados por día.

${\displaystyle {\dot {\omega }}}$ es el cambio al argumento del perigeo en grados por día.

${\displaystyle i}$ es la inclinación orbital.

${\displaystyle n}$ es el número de revoluciones orbitales por día.

La presión de la radiación solar causa perturbaciones a las órbitas. La magnitud de la aceleración que imparte a una nave espacial en órbita terrestre se modela utilizando la siguiente ecuación:^[13]​

${\displaystyle a_{R}\approx -4.5\times 10^{-6}(1+r)A/m}$

donde:

${\displaystyle a_{R}}$ es la magnitud de la aceleración en metros por segundo cuadrado.

${\displaystyle A}$ es el área de la sección transversal expuesta al Sol en metros cuadrados.

${\displaystyle m}$ es la masa de la nave espacial en kilogramos.

${\displaystyle r}$ es el factor de reflexión que depende de las propiedades del material. ${\displaystyle r=0}$ para absorción, ${\displaystyle r=1}$ para reflexión especular y ${\displaystyle r\approx 0.4}$ para reflexión difusa.

Para órbitas alrededor de la Tierra, la presión de radiación solar se convierte en una fuerza más fuerte que el arrastre por encima de 800 km de altitud.^[13]​

Hay muchos tipos diferentes de propulsión de naves espaciales. Los motores cohete son uno de los más utilizados. La fuerza de un motor cohete está modelada por la ecuación:^[15]​

${\displaystyle F_{n}={\dot {m}}\;v_{\text{e}}={\dot {m}}\;v_{\text{e-act}}+A_{\text{e}}(p_{\text{e}}-p_{\text{amb}})}$

donde:
${\displaystyle {\dot {m}}}$	=  Masa del flujo de gases de escape
${\displaystyle v_{\text{e}}}$	=  velocidad efectiva del escape
${\displaystyle v_{\text{e-act}}}$	=  velocidad del chorro real en el plano de salida de la boquilla
${\displaystyle A_{\text{e}}}$	=  área de flujo en el plano de salida de la boquilla (o el plano donde el chorro sale de la boquilla si se separa el flujo)
${\displaystyle p_{\text{e}}}$	=  presión estática en el plano de salida de la boquilla
${\displaystyle p_{\text{amb}}}$	=  presión ambiental (o atmosférica)

Otro posible método es una vela solar, que usa la presión de radiación para obtener la fuerza de propulsión deseada.^[16]​ El modelo de perturbación debido al viento solar se puede utilizar como un sistema de fuerza propulsora de una vela solar.

La fuerza primaria no gravitatoria que actúa sobre los satélites en la órbita baja de la Tierra es la resistencia atmosférica.^[13]​ El arrastre actuará en oposición a la dirección de la velocidad y eliminará la energía de la órbita. La fuerza debida al arrastre se modela mediante la siguiente ecuación:

${\displaystyle F_{D}\,=\,{\tfrac {1}{2}}\,\rho \,v^{2}\,C_{d}\,A,}$

donde

${\displaystyle \mathbf {F} _{D}}$ es la fuerza de arrastre,

${\displaystyle \mathbf {} \rho }$ es la densidad del fluido,^[17]​

${\displaystyle \mathbf {} v}$ es la velocidad del objeto respecto al fluido,

${\displaystyle \mathbf {} C_{d}}$ es el coeficiente de arrastre (un parámetro adimensional, por ejemplo, de 2 a 4 para la mayoría de los satélites^[13]​)

${\displaystyle \mathbf {} A}$ es el área de la sección de referencia.

Las órbitas con una altitud por debajo de los 120 km generalmente tienen un arrastre tan alto que se descomponen demasiado rápido para dar a un satélite la vida suficiente para cumplir cualquier misión práctica. Por otro lado, las órbitas con una altitud por encima de 600 km tienen una resistencia relativamente pequeña, por lo que la órbita se degrada lo suficientemente despacio como para no tener un impacto real en el satélite durante su vida útil.^[13]​ La densidad del aire puede variar significativamente en la termosfera, por la que circulan algunos satélites. La variación se debe principalmente a la actividad solar y, por lo tanto, este factor puede influir en gran medida en la fuerza de arrastre sobre una nave espacial y complicar la simulación de la órbita a largo plazo.^[13]​

Los campos magnéticos pueden jugar un papel importante como fuente de perturbación orbital, como se vio en la misión Long Duration Exposure Facility.^[12]​ Al igual que la gravedad, el campo magnético de la Tierra se puede expresar a través de armónicos esféricos, como se muestra a continuación:^[12]​

${\displaystyle {\mathbf {B} }=\sum _{n=1}^{\infty }\sum _{m=0}^{n}{\mathbf {B} }_{n,m}}$

donde

${\displaystyle {\mathbf {B} }}$ es el vector del campo magnético en un punto sobre la superficie de la Tierra.

${\displaystyle {\mathbf {B} }_{n,m}}$ representa la contribución a ${\displaystyle {\mathbf {B} }}$ del armónico esférico de grado n y orden m, definido como:^[12]​

${\displaystyle {\begin{aligned}\mathbf {B} _{n,m}={}&{\frac {K_{n,m}a^{n+2}}{R^{n+m+1}}}\left[{\frac {g_{n,m}{\mathcal {C}}_{m}+h_{n,m}{\mathcal {S}}_{m}}{R}}((s_{\lambda }A_{n,m+1}+(n+m+1)A_{n,m})\mathbf {\hat {r}} )-A_{n,m+1}\mathbf {\hat {e}} _{3}\right]\\[10pt]&{}-mA_{n,m}((g_{n,m}{\mathcal {C}}_{m-1}+h_{n,m}{\mathcal {S}}_{m-1})\mathbf {\hat {e}} _{1}+(h_{n,m}{\mathcal {C}}_{m-1}-g_{n,m}{\mathcal {S}}_{m-1})\mathbf {\hat {e}} _{2}))\end{aligned}}}$

donde:

${\displaystyle a}$ es el radio ecuatorial medio del cuerpo primario.

${\displaystyle R}$ es la magnitud del vector de posición desde el centro del cuerpo primario hasta el centro del cuerpo secundario.

${\displaystyle \mathbf {\hat {r}} }$ es un vector unitario en la dirección del cuerpo secundario con su origen en el centro del cuerpo primario.

${\displaystyle g_{n,m}}$ y ${\displaystyle h_{n,m}}$ son coeficientes de Gauss de grado n y orden m, que se conocen mediante medidas del campo magnético.

Los vectores unitarios ${\displaystyle \mathbf {\hat {e}} _{1},\mathbf {\hat {e}} _{2},\mathbf {\hat {e}} _{3}}$ definen un sistema de coordenadas fijado en el cuerpo primario. Para la Tierra, ${\displaystyle \mathbf {\hat {e}} _{1}}$ se encuentra en el plano ecuatorial paralelo a una línea que se cruza con el centro geométrico de la Tierra y los puntos del meridiano de Greenwich, ${\displaystyle \mathbf {\hat {e}} _{3}}$ en la dirección del eje polar norte y ${\displaystyle \mathbf {\hat {e}} _{2}=\mathbf {\hat {e}} _{3}\times \mathbf {\hat {e}} _{1}}$

${\displaystyle A_{n,m}}$ se conocen como polinomios de Legendre derivados de grado n y orden m. Se resuelven a través de la relación de recurrencia: ${\displaystyle A_{n,m}(u)={\frac {1}{n-m}}((2n-1)uA_{n-1,m}(u)-(n+m-1)A_{n-2,m}(u))}$

${\displaystyle K_{n,m}}$ se define como: 1 si m = 0, ${\displaystyle {\big [}{\frac {n-m}{n+m}}{\big ]}^{0.5}K_{n-1,m}}$ para ${\displaystyle n\geq (m+1)}$ y ${\displaystyle m=[1\ldots \infty ]}$, y ${\displaystyle [(n+m)(n-m+1)]^{-0.5}K_{n,m-1}}$ para ${\displaystyle n\geq m}$ y ${\displaystyle m=[2\ldots \infty ]}$

${\displaystyle s_{\lambda }}$ es el seno de la latitud geográfica del cuerpo secundario, que es ${\displaystyle \mathbf {\hat {r}} \cdot \mathbf {\hat {e}} _{3}}$.

${\displaystyle {\mathcal {C}}_{m},{\mathcal {S}}_{m}}$ se define con la siguiente relación de recurrencia y condiciones iniciales: ${\displaystyle {\mathcal {C}}_{m}={\mathcal {C}}_{1}{\mathcal {C}}_{m-1}-{\mathcal {S}}_{1}{\mathcal {S}}_{m-1},{\mathcal {S}}_{m}={\mathcal {S}}_{1}{\mathcal {C}}_{m-1}+{\mathcal {C}}_{1}{\mathcal {S}}_{m-1},{\mathcal {S}}_{0}=0,{\mathcal {S}}_{1}=\mathbf {R} \cdot \mathbf {\hat {e}} _{2},{\mathcal {C}}_{0}=1,{\mathcal {C}}_{1}=\mathbf {R} \cdot \mathbf {\hat {e}} _{1}}$

• Perturbación (astronomía)
• Órbita osculante
• Resonancia orbital
• Problema de los n cuerpos
• Problema de los dos cuerpos
• Esfera de influencia (astrodinámica)

• [1] Mapas de gravedad de la Tierra