When dealing with near-circular orbits I've seen the usage of the eccentricity vector defined as

$$(e_x,e_y) = (e\cos\omega, e\sin\omega)$$

I've seen this in Schaubs' "Analytical mechanics of space systems", where he says:

To avoid some numerical difficulties for near-circular orbits

Or in this paper: E/I-vector separation for safe switching of the GRACE formation

[...]the Keplerian elements $e$ and $\omega$ are commonly replaced by the eccentricity vector. It is free from singularities and well suited for the study of orbital perturbations

I think I understand where the singularity comes from: the argument of periapsis is calculated as $\cos\omega=\hat n\cdot \hat e$, where $\hat n$ is the line of nodes' versor and $\hat e =\frac{\bf{e}}{|\bf{e}|}$ is the versor of the eccentricity vector. What I don't understand is how this singularity is avoided for $e\rightarrow 0$ by defining $(e_x,e_y)$, which is still a function of $\omega$ thus still affected by the singularity.

Sorry for bad wording or formatting, I'm new on this platform, I can edit the question

  • 1
    $\begingroup$ In which calculation or equation does the singularity occur? If we see that equation and compare with the version using the eccentricity vector, we may be able to identify why the second version is singularity free. $\endgroup$
    – AJN
    Commented Feb 16, 2022 at 14:09
  • 1
    $\begingroup$ It should occour in computing $\omega$, since $\hat e = \bf{e}/|\bf{e}|$, but $e\approx 0$, so that fraction can give some numerical problems. The authors of these books/papers propose the usage of $(e_x,e_y)$ which should in their opinion be less troublesome. But this new vector still requires to compute $\omega$ EDIT: also, this new vectors also requires the eccentricity value $e$, which is affected by the same problem $\endgroup$
    – JuliusC
    Commented Feb 16, 2022 at 16:48
  • 1
    $\begingroup$ Neither this nor this answers your question directly, but it may be of some help and of interest to future readers who have landed here. $\endgroup$
    – uhoh
    Commented Feb 17, 2022 at 2:43
  • 1
    $\begingroup$ 1 "It should occur in computing $\omega$...". In a consistent equation set describing the orbital dynamics utilizing only the eccentric vector, calculation of $\omega$ might not even be required since it is not a state variable. 2 The linked pdf is not accessible without a subscription. So we can't know what equations they are using. Please post the equations relevant for orbital dynamics in the question itself using edit. 3 AFAIK, the linked pdf is actually quoting / citing a book and so the relevant equations may not be in the linked pdf; but in the cited book. $\endgroup$
    – AJN
    Commented Feb 17, 2022 at 5:10
  • $\begingroup$ 4 I could not find the book cited in the linked pdf in g__gle books. If you have access to the book, please edit in all the relevant equations from that book. $\endgroup$
    – AJN
    Commented Feb 17, 2022 at 5:20

1 Answer 1


The question is "how are singularities avoided for $e→0$ by defining $(e_x,e_y)$, which is still a function of $ω$, thus still affected by the singularity?".

(Worth noting first is that here $ω$ means argument of perigee, i.e. distance from node to perigee, and $e_x, e_y$ are rectangular components of the eccentricity vector. The formulation given in the question $$(e_x, e_y) = (e.\cos\omega, e.\sin\omega)$$ also presupposes that the xy plane is the orbit plane, also that the x-axis itself is in the direction of the node (vector) -- but these are not conditions often used in practice.)

While $e_x, e_y$ can be defined through a relation similar to that given in the question, they are not in practice calculated that way.

It would perhaps be too limited to leave the matter here, without saying how in practice such quantities have more suitably been calculated. For that, there are numerous ways, involving various different sets of orbital elements.

Relations such as given in the question can arise e.g. in problems of perturbed elliptic motion, where it is necessary to calculate relationships between many instantaneous positions and velocities of an orbiting body, and corresponding instantaneous values of some chosen set of elliptical orbital elements. It was found long ago that older sets of Keplerian elements including mean anomaly, argument of perigee and longitude of node could give rise to numerical difficulties where $e\rightarrow 0$ or $i\rightarrow 0$.

The problems, and their avoidance, can be made clearer when discussed not so much in terms of singularities, but rather of conditions under which some relations between orbital elements and position and velocity of the orbiting body become poorly determined. The 'offending' elements where the inclination and/or eccentricity are small are the locations of the perigee and node.

Thus, the relation given in the question would not be used in practice to calculate $e_x, e_y$, rectangular components of the eccentricity vector, partly because the location $ω$ of the perigee becomes poorly determined in terms of position and velocity when the magnitude $e$ of the eccentricity is close to zero. Also, the location of the node $Ω$ becomes poorly determined in terms of the orbiting body's position and velocity when the inclination is close to zero. Both those occurrences leave the relation given in the question with poor numerical definition.

For such reasons, other sets of orbital elements not involving the node, the argument of perigee/perihelion, and mean anomaly have often been used in place of the older orbital elements, so that the singularities or poor determinacies of the locations of node and perigee no longer arise.

The topic receives good explanation in A E Roy 'Orbital motion', (2005, 4th ed., ch.8, p.239 onwards), in terms that the "usual transformation to avoid this disadvantage" [i.e. the singularities] "consists of adopting the variables $h, k, p, q$ instead of the offending elements $e, ϖ, i$ and $Ω$" (where variable $ϖ$ is $(ω + Ω)$), and $$h = e.\sin ϖ$$ $$k = e.\cos ϖ$$ $$p = \sin i.\sin Ω$$ $$q = \sin i.\cos Ω .$$
Also, in place of mean anomaly $M$, mean longitude $\lambda = (M+ω+Ω)$ is used. So too are vectors $H$ and $e$, defined in terms of positions and velocities of the orbiting body in rectangular coordinates, where the instantaneous angular momentum or areal velocity, and eccentricity vector, are respectively

$H$ = ($r$ x $v$) ,

$e$ = ($v$ x $H$)/$\mu$ -$r$/r ,

with components and magnitudes respectively $H_x, H_y, H_z$ and $H_m$, and $e_x, e_y, e_z$ and $e_m$ .

As an example of principle to compare with the relation given in the question, the quantities $e_m.\cosω$ and $e_m.sinω$ would not be used, but instead the following related elements of the alternative set might be calculated when required from components of $e$ and $H$, in turn evaluated from the body's position and velocity:
$$k (= e_m.\cosϖ) = e_x - e_z.H_x/(H_m + H_z)$$ $$h (= e_m.\sinϖ) = e_y - e_z.H_y/(H_m + H_z)$$
(cf. Roy (2005) 8.30) (noting that if inclination i = 0 then e_z, H_x, and H_y are all zero as well), and these relations involve $ϖ$ i.e. $(ω + Ω)$ and do not depend on the individual parameters $ω, Ω,$ or $M$ which may become poorly-determined.

[Postscript/edit:] A visual form of explanation how the elements $ ϖ = (\omega + \Omega)$ and $\lambda = (M+\omega+\Omega)$ can remain well-determined even when the location of the node $\Omega$ and/or argument $\omega$ of perihelion become nearly indeterminate is presented by the diagram below.

Angles $ϖ$ and $\lambda$ are often portrayed as if they were (artificial) sums of angles in two different planes. In reality, they are also sums of angles in the same plane, i.e. the orbit plane, taken to define auxiliary coordinate axes, here $(F G W)$, inclined to the $(x y z)$ coordinate axes in main use (e.g. ecliptic for planetary orbits, equatorial for earth-satellites). The $F G$ plane is the orbit plane, and its intersection/node $N$ with the $x y$ plane is equally distant, by angle $\Omega = xN = FN$, from the initial points for reckoning longitudes in both planes, i.e. the $x$-vector of the $(x y z)$ coordinates and the $F$-vector of the $(F G W)$ coordinates. Argument of perigee $\omega$ is angle $NP$, and the mean position of the orbiting body is at $B$ with mean anomaly angle $PB$ = $M$.

(x y z) (F G W) coordinate diagram

It can be seen that if the inclination angle $i$ at node $N$ becomes very small, and the location of $N$ poorly determined, any error in the location of $N$ leaves the angle $FP$ practically unaffected, $XN$ and $FN$ remaining equal by construction. So the errors in angles $\omega$ and $\Omega$ are practically equal and opposite, and their sum, the distance $FP$ = $ϖ = (\omega + \Omega)$ practically unaffected. On similar considerations, $\lambda = (M + \omega + \Omega)$, angle $FB$, can remain well-determined even with poorly-determined locations of intermediate points $N$ and $P$.

The $(F G W)$ coordinate system is described for example in R M L Baker (1967) "Astrodynamics: Applications and Advanced Topics", especially Appendix C (without diagram), where $x- y- z-$ coordinates of unit vectors $W$ and $F$ are derived as

$ W_x = + \sin \Omega . \sin i $ ( = $H_x/Hm$),
$ W_y = - \cos \Omega . \sin i $ ( = $H_y/Hm$),
$ W_z = \cos i $ ( = $H_z/Hm$);

$ F_x = 1 - W_x^2 / (1 + W_z) $ ,
$ F_y = - W_x W_y / (1 + W_z) $ ,
$ F_z = - W_x $ ;

and $G$ = ($W$ x $F$).

The eccentricity vector $e$ can also be expressed by elements $k$ and $h$, seen to be its components resolved in the directions of the orthogonal $F-$ and $G-$ axes:

$(k, h) = (e_F, e_G) = (e_m.cosϖ, e_m.sinϖ)$ .

This form of relation is clearly analogous to the formulation given in the question.

It has in addition the desirable property that the quantities involved do not become poorly-determined at low inclinations and/or eccentricities, as they can be computed as shown above from the elements of vectors $H$ and $e$ derived from the orbital positions and velocities. [end edit:]

A number of related alternative sets of orbital elements have also been proposed and used for similar reasons -- see the same chapter in Roy (2005) already cited; and a previous related set proposed in Broucke & Cefola (1972), 'Equinoctial orbital elements', Cel.mech. 5:303-310, which adopts $$p = \tan i/2.\sin Ω$$ $$q = \tan i/2.\cos Ω$$ and a yet further variant set which has been used and which adopts
$$p = \sin i/2.\sin Ω$$ $$q = \sin i/2.\cos Ω .$$

[Postscript/edit:] The last-mentioned variant of elements $p$ and $q$, based on the factor $\sin i/2$, has been effectively used in celestial-mechanics studies at the Observatoire de Paris/Bureau des Longitudes (see "Introduction aux éphémérides astronomiques: supplément explicatif à la Connaissance des Temps", 1997, esp. ch.8). This set of elements has an advantage in generality for inclinations that may exceed 90° (in which case the ($p$,$q$) elements based on $\sin i$ can result in ambiguity for the value of $i$), or may reach 180° (where the elements based on the factor $\tan i/2$ encounter a singularity. [end edit:]

Brouwer & Clemence (1961) 'Methods of Celestial mechanics', ch.XI on 'Variation of Arbitrary Constants' (a/k/a variation of parameters) (sec.7 at p.287 on) has further discussion of situations where alternative forms of the orbital elements are desirable in cases of small eccentricity and/or small inclination, and how they are used.

  • $\begingroup$ Thank you for the detailed answer and for all the references, $k$ and $h$ look well behaved for near circular orbits. But aren't the equinoctial orbital elements affected by a similar problem? For (near) equatorial orbits $\Omega$ isn't well defined. Let's say I'm coding a program that computes maneuvers and trajectories using the equinoctial O.E., will I get into some numerical problems when dealing with a plane change that equatorializes my orbit? $\endgroup$
    – JuliusC
    Commented Feb 17, 2022 at 22:23
  • 1
    $\begingroup$ @JuliusC - Yes, true that in near-equatorial orbits Ω is poorly-defined, but ω+Ω=ϖ is better-defined than Ω provided ω and Ω are evaluated from the same r- and v- data & intermediates. The Ω node vector has a v small z-component that hardly affects ϖ in practice, so ϖ is as well-determined as the size of e allows. Then if e is also v small, λ=(M+ω+Ω) is as well-determined as may be, under the same proviso. The k, h, q, p elements help to keep these indeterminacies under control. I'll return on the 'equinoctial' elements, the q/p variants of them mentioned last are more often used. $\endgroup$
    – terry-s
    Commented Feb 17, 2022 at 23:02
  • 1
    $\begingroup$ @JuliusC - While I don't have details of your problem, I suggest you get hold of an edition of Roy's 'Orbital motion' and use as needed his relations between perturbation variables, rectangular posn & vel coordinates and conic-section elements. If you want to calculate perturbations of the elements & progress the motion & orbit by variation of parameters, the fmlae fit well with his preceding derivations of perturbation eqns (it can be tricky with some parameter sets to correctly implement the needed additional equation for the integration of the mean motion.) Good luck. $\endgroup$
    – terry-s
    Commented Feb 17, 2022 at 23:41
  • 1
    $\begingroup$ @JuliusC -- Thanks for your response: I just added some material to give visual explanation of the point we already discussed in our comments. $\endgroup$
    – terry-s
    Commented Feb 19, 2022 at 14:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.