Editing Orbital mechanics (section)

==Calculating trajectories==

===Kepler's equation===
One approach to calculating orbits (mainly used historically) is to use [[Kepler's equation]]:
:<math> M = E - \epsilon \cdot \sin E </math>.

where ''M'' is the [[mean anomaly]], ''E'' is the [[eccentric anomaly]], and <math> \epsilon </math> is the [[eccentricity (mathematics)|eccentricity]].

With Kepler's formula, finding the time-of-flight to reach an angle ([[true anomaly]]) of <math>\theta</math> from [[periapsis]] is broken into two steps:
# Compute the eccentric anomaly <math>E</math> from true anomaly <math>\theta</math>
# Compute the time-of-flight <math>t</math> from the eccentric anomaly <math>E</math>

Finding the eccentric anomaly at a given time ([[Kepler's Equation#Inverse problem|the inverse problem]]) is more difficult. Kepler's equation is [[transcendental function|transcendental]] in <math>E</math>, meaning it cannot be solved for <math>E</math> [[algebraic function|algebraically]]. Kepler's equation can be solved for <math>E</math> [[analytic function|analytically]] by inversion.

A solution of Kepler's equation, valid for all real values of <math> \textstyle \epsilon </math> is:

<math display="block">
E =
\begin{cases}

\displaystyle \sum_{n=1}^{\infty}
{\frac{M^{\frac{n}{3}}}{n!}} \lim_{\theta \to 0} \left(
\frac{\mathrm{d}^{\,n-1}}{\mathrm{d}\theta^{\,n-1}} \left[ \left(
\frac{\theta}{ \sqrt[3]{\theta - \sin(\theta)} } \right) ^n \right]
\right)
,  & \epsilon = 1  \\

\displaystyle \sum_{n=1}^{\infty}
{ \frac{ M^n }{ n! } }
\lim_{\theta \to 0} \left(
\frac{\mathrm{d}^{\,n-1}}{\mathrm{d}\theta^{\,n-1}} \left[ \left(
\frac{ \theta }{ \theta - \epsilon \cdot \sin(\theta)} \right) ^n \right]
\right)
, &  \epsilon \ne  1

\end{cases}
</math>

Evaluating this yields:

<math display="block">
E =
\begin{cases} \displaystyle
x + \frac{1}{60} x^3 + \frac{1}{1400}x^5 + \frac{1}{25200}x^7 + \frac{43}{17248000}x^9 + \frac{ 1213}{7207200000 }x^{11} +
\frac{151439}{12713500800000 }x^{13} \cdots \ | \ x = ( 6 M )^\frac{1}{3}
,  & \epsilon = 1  \\
\\
\displaystyle
\frac{1}{1-\epsilon} M
- \frac{\epsilon}{( 1-\epsilon)^4 } \frac{M^3}{3!}
+ \frac{(9 \epsilon^2 + \epsilon)}{(1-\epsilon)^7 } \frac{M^5}{5!}
- \frac{(225 \epsilon^3 + 54 \epsilon^2 + \epsilon ) }{(1-\epsilon)^{10} } \frac{M^7}{7!}
+ \frac{ (11025\epsilon^4 + 4131 \epsilon^3 + 243 \epsilon^2 + \epsilon ) }{(1-\epsilon)^{13} } \frac{M^9}{9!} \cdots

, &  \epsilon \ne  1

\end{cases}
</math>

<br/>Alternatively, Kepler's Equation can be solved numerically.  First one must guess a value of <math>E</math> and solve for time-of-flight; then adjust <math>E</math> as necessary to bring the computed time-of-flight closer to the desired value until the required precision is achieved.  Usually, [[Newton's method]] is used to achieve relatively fast convergence.

The main difficulty with this approach is that it can take prohibitively long to converge for the extreme elliptical orbits.  For near-parabolic orbits, eccentricity <math>\epsilon</math> is nearly 1, and substituting <math>e = 1</math> into the formula for mean anomaly, <math>E - \sin E</math>, we find ourselves subtracting two nearly-equal values, and accuracy suffers.  For near-circular orbits, it is hard to find the periapsis in the first place (and truly circular orbits have no periapsis at all).  Furthermore, the equation was derived on the assumption of an elliptical orbit, and so it does not hold for parabolic or hyperbolic orbits.  These difficulties are what led to the development of the [[universal variable formulation]], described below.

===Conic orbits===
For simple procedures, such as computing the [[delta-v]] for coplanar transfer ellipses, traditional approaches{{Clarify|date=February 2009}} are fairly effective. Others, such as time-of-flight are far more complicated, especially for near-circular and hyperbolic orbits.

===The patched conic approximation===
{{Main|Patched conic approximation}}
The [[Hohmann transfer orbit]] alone is a poor approximation for interplanetary trajectories because it neglects the planets' own gravity.  Planetary gravity dominates the behavior of the spacecraft in the vicinity of a planet and in most cases Hohmann severely overestimates delta-v, and produces highly inaccurate prescriptions for burn timings. A relatively simple way to get a [[orders of approximation|first-order approximation]] of delta-v is based on the 'Patched Conic Approximation' technique. One must choose the one dominant gravitating body in each region of space through which the trajectory will pass, and to model only that body's effects in that region.  For instance, on a trajectory from the Earth to Mars, one would begin by considering only the Earth's gravity until the trajectory reaches a distance where the Earth's gravity no longer dominates that of the Sun.  The spacecraft would be given [[escape velocity]] to send it on its way to interplanetary space.  Next, one would consider only the Sun's gravity until the trajectory reaches the neighborhood of Mars.  During this stage, the transfer orbit model is appropriate.  Finally, only Mars's gravity is considered during the final portion of the trajectory where Mars's gravity dominates the spacecraft's behavior.  The spacecraft would approach Mars on a hyperbolic orbit, and a final retrograde burn would slow the spacecraft enough to be captured by Mars. [[Friedrich Zander]] was one of the first to apply the patched-conics approach for astrodynamics purposes, when proposing the use of intermediary bodies' gravity for interplanetary travels, in what is known today as a [[gravity assist]].<ref>{{cite journal |last1=Negri |first1=Rodolfo Batista |last2=Prado |first2=Antônio Fernando Bertachini de Alme |title=A historical review of the theory of gravity-assists in the pre-spaceflight era |journal=Journal of the Brazilian Society of Mechanical Sciences and Engineering |date=August 2020 |volume=42 |issue=8 |pages=406 |doi=10.1007/s40430-020-02489-x|s2cid=220510617 |url=http://urlib.net/8JMKD3MGP3W34R/42T4NAH }}</ref>

The size of the "neighborhoods" (or [[sphere of influence (astrodynamics)|spheres of influence]]) vary with radius <math>r_{SOI}</math>:
:<math>r_{SOI} = a_p\left(\frac{m_p}{m_s}\right)^{2/5}</math>
where <math>a_p</math> is the [[semimajor axis]] of the planet's orbit relative to the [[Sun]]; <math>m_p</math> and <math>m_s</math> are the [[mass]]es of the planet and Sun, respectively.

This simplification is sufficient to compute rough estimates of fuel requirements, and rough time-of-flight estimates, but it is not generally accurate enough to guide a spacecraft to its destination.  For that, numerical methods are required.
<!--
====Derivation of formulas====
For a spacecraft to travel from one planet to another, the patched conic approximation breaks the trip into three phases:
#Escape from the origin planet
#Interplanetary transfer
#Capture by the destination planet

We begin by considering phase 2.  In this phase, the gravity of the planets is neglected, and only the Sun's gravity is considered.  Therefore, this phase can be treated as a standard Hohmann transfer orbit.

The [[specific orbital energy]] of an orbit is <math>-\mu/2a</math>, where <math>a</math> is the [[semimajor axis]] and <math>\mu</math> is the [[standard gravitational parameter]] of the gravitating body.  For phase 2, we are interested in the following orbits around the Sun:

{| border=1
|-
! Orbit              !! Major axis                 !! Total specific orbital energy
|-
| At origin          || <math> 2r_0 \,</math>      || <math>-\frac{\mu}{2r_0}</math>
|-
| Transfer           || <math> r_0 + r_1 \,</math> || <math>-\frac{\mu}{r_0+r_1}</math>
|-
| At destination     || <math> 2r_1 \,</math>      || <math>-\frac{\mu}{2r_1}</math>
|}

*<math>-\mu/2r_0</math> for the origin planet
*<math>-\mu/2r_1</math> for the destination planet
*<math>-\mu/(r_0+r_1)</math> for the Hohmann transfer orbit

We must first get from the origin planet's orbit to the Hohmann transfer orbit, which requires a delta-v sufficient to make up the difference in orbital energy between the two orbits involves a change in specific kinetic energy of

A useful concept is the speed of a circular orbit speed at a given distance from a given body, which we will refer to as <math>C_b(r)</math>, where <math>b</math> is the body and <math>r</math> is the distance.  (We may omit <math>b</math> and <math>r</math> for clarity.)  This speed is usually easy to determine, simplifies the formulas, and tells us a lot about the gravitational conditions of a location in space.  The [[specific orbital energy]] for a body in a circular orbit is <math>-\frac{1}{2C}</math>, and the [[specific kinetic energy]] is <math>\frac{1}{2C}</math>.
-->

===The universal variable formulation===
To address computational shortcomings of traditional approaches for solving the 2-body problem, the [[universal variable formulation]] was developed.  It works equally well for the circular, elliptical, parabolic, and hyperbolic cases, the differential equations converging well when integrated for any orbit.  It also generalizes well to problems incorporating perturbation theory.
<!-- TODO: Describe and possibly derive the universal variable formulation.  Talk about the f and g functions, and the X variable, and maybe outline an algorithm, maybe in wikicode.  Try not to just lift this from Fundamentals of Astrophysics text. -->

===Perturbations===
The universal variable formulation works well with the variation of parameters technique, except now, instead of the six Keplerian orbital elements, we use a different set of orbital elements: namely, the satellite's initial position and velocity vectors <math>x_0</math> and <math>v_0</math> at a given epoch <math>t = 0</math>.  In a two-body simulation, these elements are sufficient to compute the satellite's position and velocity at any time in the future, using the universal variable formulation.  Conversely, at any moment in the satellite's orbit, we can measure its position and velocity, and then use the universal variable approach to determine what its initial position and velocity ''would have been'' at the epoch.  In perfect two-body motion, these orbital elements would be invariant (just like the Keplerian elements would be).

However, perturbations cause the orbital elements to change over time. Hence, the position element is written as <math>x_0(t)</math> and the velocity element as <math>v_0(t)</math>, indicating that they vary with time.  The technique to compute the effect of perturbations becomes one of finding expressions, either exact or approximate, for the functions <math>x_0(t)</math> and <math>v_0(t)</math>.<!-- TODO: Explain it more -->

The following are some effects which make real orbits differ from the simple models based on a spherical Earth.  Most of them can be handled on short timescales (perhaps less than a few thousand orbits) by perturbation theory because they are small relative to the corresponding two-body effects.
*Equatorial bulges cause [[precession]] of the node and the perigee
*[[Spherical harmonics#Visualization of the spherical harmonics|Tesseral harmonic]]s<ref>{{MathWorld |title=Tesseral Harmonic |id=TesseralHarmonic |access-date=2019-10-07}}</ref> of the gravity field introduce additional perturbations
*Lunar and solar gravity perturbations alter the orbits
*Atmospheric drag reduces the semi-major axis unless make-up thrust is used

Over very long timescales (perhaps millions of orbits), even small perturbations can dominate, and the behavior can become [[Chaos theory|chaotic]]. On the other hand, the various perturbations can be orchestrated by clever astrodynamicists to assist with orbit maintenance tasks, such as [[Orbital station-keeping|station-keeping]], [[ground track]] maintenance or adjustment, or phasing of perigee to cover selected targets at low altitude.