Plotting Ground Track of Elliptical Orbit with Inclination Angle

Question

I am trying to project the ground path of an Elliptical orbit with an arbitrary inclination angle of $\Delta i$ and a period of $T$. I have all the parameters of the orbit (semi-major axis length, semi-minor axis length, eccentricity, apogee and perigee velocity, energy, etc.).

The approach I was taking was plotting a position of the satellite in the elliptical orbit on an X-Y plane. I iterate through the time domain $t \in [0,\textrm{T}]$ where T is the period of the orbit with a specified time step of $\Delta t$. For each time, I calculate the mean anomaly and the magnitude of the distance between the satellite and the center body, $r$. When projecting the position vector of the satellite with respect to the center body on the semi-major axis, the length of that projected vector will be defined as $x$. After finding $x$, I found the vertical projection of that same position vector, thus forming the triangle below.

First of all, I need to find mean anomaly as a function of time. My code is currently not working, but in theory what I am doing is (from Curtis's Orbital Mechanics for Engineering Students, 1st Edition) I am solving numerically for a quantity labeled as $E$ using the equation for mean anomaly (eqn. between 3.12 and 3.13 for those who have the book), given as: $M_e = E - e\sin(E)$. From there, I solve numerically for true anomaly using eqn. 3.7a:

$M_e = 2\arctan{(\sqrt{\frac{1-e}{1+e}}\tan{\theta/2})} - \frac{e\sqrt{1-e^2}\sin\theta}{1+e\cos\theta}$.

After finding those distances, I found the spherical coordinates of the satellite with respect to the center of the center body. $r$ would be $r_p + x$, $\phi$ would simply be the complement of $\Delta i$. $\theta$ would be the true anomaly.

After that, I plot (x,y) using spherical-cartesian transformation, but I simply do not believe that it is that straightforward. For one, what I am essentially trying to do is projecting a 3D curve onto a sphere and then projecting that onto a 2D plane. If anyone has sources/reading material on how to actually do this math, that would be fantastic. I was referencing this chapter from a book on satellites: http://fgg-web.fgg.uni-lj.si/~/mkuhar/Pouk/SG/Seminar/Vrste_tirnic_um_Zemljinih_sat/Orbit_and_Ground_Track_of_a_Satellite-Capderou2005.pdf. Unfortunately I cannot reference other equations that the book references until I get it.

I have not gotten the code to run yet, so I do not know what my plots look like but from preliminary data, it does not seem to be going in the right way (e.g. True anomaly ($\theta$) should be going from 0 to $2\pi$ as time increases, but that is not what is outputting.

EDIT: Typos and to clarify, I do not want to use external packages such as STK. I want to be able to physically come up with these plots.

Answer 1

Let's first attach a right-handed coordinate system to the orbit, defined by the unit basis vectors $(\boldsymbol{e}_x, \boldsymbol{e}_y, \boldsymbol{e}_z)$, with $\boldsymbol{e}_x$ and $\boldsymbol{e}_y$ in the orbital plane and $\boldsymbol{e}_x$ in the direction of periapsis. First, compute the eccentric anomaly $E$ from the mean anomaly $M$: $$ M = \frac{2\pi t}{T} = E - \sin E, $$ assuming that the satellite is at periapsis at $t=0$. The position vector of the satellite is then given by $$ \boldsymbol{r} = a(\cos E -e)\,\boldsymbol{e}_x + b\sin E\,\boldsymbol{e}_y. $$ There's no need to compute the true anomaly.

Next we need to find the general transformation from the orbital frame to a given reference frame, defined by $(\boldsymbol{e}_x', \boldsymbol{e}_y', \boldsymbol{e}_z')$. This can be done using three Euler angles $\Omega$, $i$, and $\omega$:

• $\Omega$ is the longitude of the ascending node in the ($\boldsymbol{e}_x'$, $\boldsymbol{e}_y'$) plane, i.e. the angle measured from $\boldsymbol{e}_x'$ to the ascending node;
• $i$ is the inclination of the orbit with the ($\boldsymbol{e}_x'$, $\boldsymbol{e}_y'$) plane;
• $\omega$ is the argument of periapsis, i.e. the angle measured from the ascending node to the periapsis.

See also the wiki article on orbital elements. The transformation is then given by three rotations: $$ \begin{pmatrix} \boldsymbol{e}_x\\ \boldsymbol{e}_y\\ \boldsymbol{e}_z \end{pmatrix} \!=\! \begin{pmatrix} \cos\omega & \!\sin\omega & \!0 \\ -\sin\omega & \!\cos\omega & \!0 \\ 0 & \!0 & \!1 \end{pmatrix}\! \begin{pmatrix} 1 & \!0 & \!0 \\ 0 & \!\cos i & \!\sin i \\ 0 & \!-\sin i & \!\cos i \end{pmatrix}\! \begin{pmatrix} \cos\Omega & \!\sin\Omega & \!0 \\ -\sin\Omega & \!\cos\Omega & \!0 \\ 0 & \!0 & \!1 \end{pmatrix}\! \begin{pmatrix} \boldsymbol{e}_x'\\ \boldsymbol{e}_y'\\ \boldsymbol{e}_z' \end{pmatrix} $$ or explicitly, $$ \boldsymbol{e}_x = (\cos\Omega\cos\omega - \sin\Omega\sin\omega\cos i)\,\boldsymbol{e}_x' + (\sin\Omega\cos\omega + \cos\Omega\sin\omega\cos i)\,\boldsymbol{e}_y' + (\sin\omega\sin i)\,\boldsymbol{e}_z'\\[1em] \boldsymbol{e}_y = (-\cos\Omega\sin\omega - \sin\Omega\cos\omega\cos i)\,\boldsymbol{e}_x' + (-\sin\Omega\sin\omega + \cos\Omega\cos\omega\cos i)\,\boldsymbol{e}_y' + (\cos\omega\sin i)\,\boldsymbol{e}_z'\\[1em] \boldsymbol{e}_z = (\sin\Omega\sin i)\,\boldsymbol{e}_x' + (-\cos\Omega\sin i)\,\boldsymbol{e}_y' + \cos i\,\boldsymbol{e}_z'. $$ This enables you to express $\boldsymbol{r}$ in the reference frame coordinates. For the projected path, simply set the $z'$ component to zero.

Answer 2

How /is/ True Anomaly varying in your results? If it's varying from $-\pi$ to $\pi$, it might be okay, and merely a result of the computer restricting the ranges of the various trigonometric functions automatically.

You shouldn't need to solve numerically for True Anomaly ($\theta$) from Mean Anomaly (M). It is fairly straightforward to pull it from Eccentric Anomaly (E), and Orbital Eccentricity (e) once you have those values.

When I was working on a program to run Kepler's equations a few years ago, there were several things I did: First and foremost, for elliptical orbits, I would rectify the value I had for Mean Anomaly into the range $(-\pi,\pi]$ , which allowed for more code reuse when I did the version that dealt with hyperbolic orbits.

Given that assumption, the range on Eccentric Anomaly will also be $(-\pi,\pi]$, with the same signs as Mean Anomaly.

And that allows you to pull True Anomaly on the elliptical orbit with the following equation:

$$\theta = \pm arctan\left({\frac{\sqrt{1+e}}{\sqrt{1-e}}tan\frac{E}{2}}\right)$$

Pick the value of True Anomaly that has the same sign that the Eccentric Anomaly did. The result will, as with the other Anomalies, be in the range $(-\pi , \pi]$

If you truly need the value of True Anomaly in the range $(0,2\pi)$, you can subsequently add $2\pi$ to all negative values of $\theta$.