Orbit modeling

Cowell's method. Forces from all perturbing bodies (black and gray) are summed to form the total force on body i (red), and this is numerically integrated starting from the initial position (the epoch of osculation).

Cowell's method is a special perturbation method;[3] mathematically, for ${\displaystyle n}$ mutually interacting bodies, Newtonian forces on body ${\displaystyle i}$ from the other bodies ${\displaystyle j}$ are simply summed thus,

${\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}}}$

where

${\displaystyle \mathbf {\ddot {r}} _{i}}$ is the acceleration vector of body ${\displaystyle i}$

${\displaystyle G}$ is the gravitational constant

${\displaystyle m_{j}}$ is the mass of body ${\displaystyle j}$

${\displaystyle \mathbf {r} _{i}}$ and ${\displaystyle \mathbf {r} _{j}}$ are the position vectors of objects ${\displaystyle i}$ and ${\displaystyle j}$

${\displaystyle r_{ij}}$ is the distance from object ${\displaystyle i}$ to object ${\displaystyle j}$

with all vectors being referred to the barycenter of the system. This equation is resolved into components in ${\displaystyle x}$, ${\displaystyle y}$, ${\displaystyle z}$ and these are integrated numerically to form the new velocity and position vectors as the simulation moves forward in time. The advantage of Cowell's method is ease of application and programming. A disadvantage is that when perturbations become large in magnitude (as when an object makes a close approach to another) the errors of the method also become large.[4] Another disadvantage is that in systems with a dominant central body, such as the Sun, it is necessary to carry many significant digits in the arithmetic because of the large difference in the forces of the central body and the perturbing bodies.[5]

Encke's method. Greatly exaggerated here, the small difference δr (blue) between the osculating, unperturbed orbit (black) and the perturbed orbit (red), is numerically integrated starting from the initial position (the epoch of osculation).

Encke's method begins with the osculating orbit as a reference and integrates numerically to solve for the variation from the reference as a function of time.[6] Its advantages are that perturbations are generally small in magnitude, so the integration can proceed in larger steps (with resulting lesser errors), and the method is much less affected by extreme perturbations than Cowell's method. Its disadvantage is complexity; it cannot be used indefinitely without occasionally updating the osculating orbit and continuing from there, a process known as rectification.[4] [7]

Letting ${\displaystyle {\boldsymbol {\rho }}}$ be the radius vector of the osculating orbit, ${\displaystyle \mathbf {r} }$ the radius vector of the perturbed orbit, and ${\displaystyle \delta \mathbf {r} }$ the variation from the osculating orbit,

${\displaystyle \delta \mathbf {r} =\mathbf {r} -{\boldsymbol {\rho }},}$ and the equation of motion of ${\displaystyle \delta \mathbf {r} }$ is simply

(1)

${\displaystyle {\ddot {\delta \mathbf {r} }}=\mathbf {\ddot {r}} -{\boldsymbol {\ddot {\rho }}}.}$

(2)

${\displaystyle \mathbf {\ddot {r}} }$ and ${\displaystyle {\boldsymbol {\ddot {\rho }}}}$ are just the equations of motion of ${\displaystyle \mathbf {r} }$ and ${\displaystyle {\boldsymbol {\rho }}}$,

${\displaystyle \mathbf {\ddot {r}} =\mathbf {a} _{\text{per}}-{\mu \over r^{3}}\mathbf {r} }$ for the perturbed orbit and

(3)

${\displaystyle {\boldsymbol {\ddot {\rho }}}=-{\mu \over \rho ^{3}}{\boldsymbol {\rho }}}$ for the unperturbed orbit,

(4)

where ${\displaystyle \mu =G(M+m)}$ is the gravitational parameter with ${\displaystyle M}$ and ${\displaystyle m}$ the masses of the central body and the perturbed body, ${\displaystyle \mathbf {a} _{\text{per}}}$ is the perturbing acceleration, and ${\displaystyle r}$ and ${\displaystyle \rho }$ are the magnitudes of ${\displaystyle \mathbf {r} }$ and ${\displaystyle {\boldsymbol {\rho }}}$.

Substituting from equations (3) and (4) into equation (2),

${\displaystyle {\ddot {\delta \mathbf {r} }}=\mathbf {a} _{\text{per}}+\mu \left({{\boldsymbol {\rho }} \over \rho ^{3}}-{\mathbf {r} \over r^{3}}\right),}$

(5)

which, in theory, could be integrated twice to find ${\displaystyle \delta \mathbf {r} }$. Since the osculating orbit is easily calculated by two-body methods, ${\displaystyle {\boldsymbol {\rho }}}$ and ${\displaystyle \delta \mathbf {r} }$ are accounted for and ${\displaystyle \mathbf {r} }$ can be solved. In practice, the quantity in the brackets, ${\displaystyle {{\boldsymbol {\rho }} \over \rho ^{3}}-{\mathbf {r} \over r^{3}}}$, is the difference of two nearly equal vectors, and further manipulation is necessary to avoid the need for extra significant digits.[8][9]

In 1991 Victor R. Bond and Michael F. Fraietta created an efficient and highly accurate method for solving the two-body perturbed problem.[10] This method uses the linearized and regularized differential equations of motion derived by Hans Sperling and a perturbation theory based on these equations developed by C.A. Burdet in the year 1864. In 1973, Bond and Hanssen improved Burdet's set of differential equations by using the total energy of the perturbed system as a parameter instead of the two-body energy and by reducing the number of elements to 13. In 1989 Bond and Gottlieb embedded the Jacobian integral, which is a constant when the potential function is explicitly dependent upon time as well as position in the Newtonian equations. The Jacobian constant was used as an element to replace the total energy in a reformulation of the differential equations of motion. In this process, another element which is proportional to a component of the angular momentum is introduced. This brought the total number of elements back to 14. In 1991, Bond and Fraietta made further revisions by replacing the Laplace vector with another vector integral as well as another scalar integral which removed small secular terms which appeared in the differential equations for some of the elements.[11]

The Sperling–Burdet method is executed in a 5 step process as follows:[11]

Step 1: Initialization

Given an initial position, ${\displaystyle \mathbf {r} _{0}}$, an initial velocity, ${\displaystyle \mathbf {v} _{0}}$, and an initial time, ${\displaystyle t_{0}}$, the following variables are initialized:

${\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}}$

Perturbations due to perturbing masses, defined as ${\displaystyle V_{0}}$ and ${\displaystyle {\Bigg [}{\partial {V} \over {\partial {\mathbf {r} }}}{\Bigg ]}_{0}}$, are evaluated

Perturbations due to other accelerations, defined as ${\displaystyle \mathbf {P} _{0}}$, are evaluated

${\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}$

Step 2: Transform elements to coordinates

${\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}}$

where ${\displaystyle c_{0},c_{1},c_{2},c_{3}}$ are Stumpff functions

Step 3: Evaluate differential equations for the elements

${\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 ]}}$

Step 4: Integration

Here the differential equations are integrated over a period ${\displaystyle \Delta s}$ to obtain the element value at ${\displaystyle s+\Delta s}$

Step 5: Advance

Set ${\displaystyle s=s+\Delta s}$ and return to step 2 until simulation stopping conditions are met.

The Earth is not a perfect sphere nor is mass evenly distributed within the Earth. This results in the point-mass gravity model being inaccurate for orbits around the Earth, particularly Low Earth orbits. To account for variations in gravitational potential around the surface of the Earth, the gravitational field of the Earth is modeled with spherical harmonics[12] which are expressed through the equation:

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

where

${\displaystyle {\mu }}$ is the gravitational parameter defined as the product of G, the universal gravitational constant, and the mass of the primary body.

${\displaystyle \mathbf {\hat {r}} }$ is the unit vector defining the distance between the primary and secondary bodies, with ${\displaystyle {R}}$ being the magnitude of the distance.

${\displaystyle {\mathbf {f} }_{n,m}}$ represents the contribution to ${\displaystyle {\mathbf {f} }}$ of the spherical harmonic of degree n and order m, which is defined as:[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}}}$

where:

${\displaystyle R_{O}}$ is the mean equatorial radius of the primary body.

${\displaystyle R}$ is the magnitude of the position vector from the center of the primary body to the center of the secondary body.

${\displaystyle C_{n,m}}$ and ${\displaystyle S_{n,m}}$ are gravitational coefficients of degree n and order m. These are typically found through gravimetry measurements.

The unit vectors ${\displaystyle \mathbf {\hat {e}} _{1},\mathbf {\hat {e}} _{2},\mathbf {\hat {e}} _{3}}$ define a coordinate system fixed on the primary body. For the Earth, ${\displaystyle \mathbf {\hat {e}} _{1}}$ lies in the equatorial plane parallel to a line intersecting Earth's geometric center and the Greenwich meridian,${\displaystyle \mathbf {\hat {e}} _{3}}$ points in the direction of the North polar axis, and ${\displaystyle \mathbf {\hat {e}} _{2}=\mathbf {\hat {e}} _{3}\times \mathbf {\hat {e}} _{1}}$

${\displaystyle A_{n,m}}$ is referred to as a derived Legendre polynomial of degree n and order m. They are solved through the recurrence relation: ${\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 }}$ is sine of the geographic latitude of the secondary body, which is ${\displaystyle \mathbf {\hat {r}} \cdot \mathbf {\hat {e}} _{3}}$.

${\displaystyle {\mathcal {C}}_{m},{\mathcal {S}}_{m}}$ are defined with the following recurrence relation and initial conditions:${\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}}$

When modeling perturbations of an orbit around a primary body only the sum of the ${\displaystyle {\mathbf {f} }_{n,m}}$ terms need to be included in the perturbation since the point-mass gravity model is accounted for in the ${\displaystyle -{\frac {\mu }{R^{2}}}\mathbf {\hat {r}} }$ term

Gravitational forces from third bodies can cause perturbations to an orbit. For example, the Sun and Moon cause perturbations to Orbits around the Earth.[13] These forces are modeled in the same way that gravity is modeled for the primary body by means of direct gravitational N-body simulations. Typically, only a spherical point-mass gravity model is used for modeling effects from these third bodies.[14] Some special cases of third-body perturbations have approximate analytic solutions. For example, perturbations for the right ascension of the ascending node and argument of perigee for a circular Earth orbit are:[13]

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

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

where:

${\displaystyle {\dot {\Omega }}}$ is the change to the right ascension of the ascending node in degrees per day.

${\displaystyle {\dot {\omega }}}$ is the change to the argument of perigee in degrees per day.

${\displaystyle i}$ is the orbital inclination.

${\displaystyle n}$ is the number of orbital revolutions per day.

Solar radiation pressure causes perturbations to orbits. The magnitude of acceleration it imparts to a spacecraft in Earth orbit is modeled using the equation below:[13]

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

where:

${\displaystyle a_{R}}$ is the magnitude of acceleration in meters per second-squared.

${\displaystyle A}$ is the cross-sectional area exposed to the Sun in meters-squared.

${\displaystyle m}$ is the spacecraft mass in kilograms.

${\displaystyle r}$ is the reflection factor which depends on material properties. ${\displaystyle r=0}$ for absorption, ${\displaystyle r=1}$ for specular reflection, and ${\displaystyle r\approx 0.4}$ for diffuse reflection.

For orbits around the Earth, solar radiation pressure becomes a stronger force than drag above 800 km altitude.[13]

There are many different types of spacecraft propulsion. Rocket engines are one of the most widely used. The force of a rocket engine is modeled by the equation:[15]

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

where:
${\displaystyle {\dot {m}}}$	=  exhaust gas mass flow
${\displaystyle v_{\text{e}}}$	=  effective exhaust velocity
${\displaystyle v_{\text{e-act}}}$	=  actual jet velocity at nozzle exit plane
${\displaystyle A_{\text{e}}}$	=  flow area at nozzle exit plane (or the plane where the jet leaves the nozzle if separated flow)
${\displaystyle p_{\text{e}}}$	=  static pressure at nozzle exit plane
${\displaystyle p_{\text{amb}}}$	=  ambient (or atmospheric) pressure

Another possible method is a solar sail. Solar sails use radiation pressure in a way to achieve a desired propulsive force.[16] The perturbation model due to the solar wind can be used as a model of propulsive force from a solar sail.

The primary non-gravitational force acting on satellites in low Earth orbit is atmospheric drag.[13] Drag will act in opposition to the direction of velocity and remove energy from an orbit. The force due to drag is modeled by the following equation:

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

where

${\displaystyle \mathbf {F} _{D}}$ is the force of drag,

${\displaystyle \mathbf {} \rho }$ is the density of the fluid,[17]

${\displaystyle \mathbf {} v}$ is the velocity of the object relative to the fluid,

${\displaystyle \mathbf {} C_{d}}$ is the drag coefficient (a dimensionless parameter, e.g. 2 to 4 for most satellites[13])

${\displaystyle \mathbf {} A}$ is the reference area.

Orbits with an altitude below 120 km generally have such high drag that the orbits decay too rapidly to give a satellite a sufficient lifetime to accomplish any practical mission. On the other hand, orbits with an altitude above 600 km have relatively small drag so that the orbit decays slow enough that it has no real impact on the satellite over its useful life.[13] Density of air can vary significantly in the thermosphere where most low Earth orbiting satellites reside. The variation is primarily due to solar activity, and thus solar activity can greatly influence the force of drag on a spacecraft and complicate long-term orbit simulation.[13]

Magnetic fields can play a significant role as a source of orbit perturbation as was seen in the Long Duration Exposure Facility.[12] Like gravity, the magnetic field of the Earth can be expressed through spherical harmonics as shown below:[12]

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

where

${\displaystyle {\mathbf {B} }}$ is the magnetic field vector at a point above the Earth's surface.

${\displaystyle {\mathbf {B} }_{n,m}}$ represents the contribution to ${\displaystyle {\mathbf {B} }}$ of the spherical harmonic of degree n and order m, defined as:[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}}}$

where:

${\displaystyle a}$ is the mean equatorial radius of the primary body.

${\displaystyle R}$ is the magnitude of the position vector from the center of the primary body to the center of the secondary body.

${\displaystyle \mathbf {\hat {r}} }$ is a unit vector in the direction of the secondary body with its origin at the center of the primary body.

${\displaystyle g_{n,m}}$ and ${\displaystyle h_{n,m}}$ are Gauss coefficients of degree n and order m. These are typically found through magnetic field measurements.

The unit vectors ${\displaystyle \mathbf {\hat {e}} _{1},\mathbf {\hat {e}} _{2},\mathbf {\hat {e}} _{3}}$ define a coordinate system fixed on the primary body. For the Earth, ${\displaystyle \mathbf {\hat {e}} _{1}}$ lies in the equatorial plane parallel to a line intersecting Earth's geometric center and the Greenwich meridian,${\displaystyle \mathbf {\hat {e}} _{3}}$ points in the direction of the North polar axis, and ${\displaystyle \mathbf {\hat {e}} _{2}=\mathbf {\hat {e}} _{3}\times \mathbf {\hat {e}} _{1}}$

${\displaystyle A_{n,m}}$ is referred to as a derived Legendre polynomial of degree n and order m. They are solved through the recurrence relation: ${\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}}$ is defined as: 1 if m = 0, ${\displaystyle {\big [}{\frac {n-m}{n+m}}{\big ]}^{0.5}K_{n-1,m}}$ for ${\displaystyle n\geq (m+1)}$ and ${\displaystyle m=[1\ldots \infty ]}$ , and ${\displaystyle [(n+m)(n-m+1)]^{-0.5}K_{n,m-1}}$ for ${\displaystyle n\geq m}$ and ${\displaystyle m=[2\ldots \infty ]}$

${\displaystyle s_{\lambda }}$ is sine of the geographic latitude of the secondary body, which is ${\displaystyle \mathbf {\hat {r}} \cdot \mathbf {\hat {e}} _{3}}$.

${\displaystyle {\mathcal {C}}_{m},{\mathcal {S}}_{m}}$ are defined with the following recurrence relation and initial conditions:${\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}}$