I would like to simulate a satellite's orbit iteratively using Cartesian coordinates. I am using the center of the earth for $(0, 0, 0).$ Given a heading, pitch, magnitude, and initial position, I want to compute the velocity vector.
Given that:
$$\begin{eqnarray}
\mbox{radius} &=& \sqrt{x^2 + y^2 + z^2}\\
\mbox{latitude} &=& \arccos(z / \mbox{radius})\\
\mbox{longitude} &=& \arctan(y / x)
\end{eqnarray}$$
For my polar coordinates, I assume the equator has latitude of $\pi/2,$ and the north pole is latitude $0.$
Right now I have unit velocity vector set with these equations, however they assume heading is $0$ and pitch is $0$ (I hope):
$$\begin{eqnarray}
V_x &=& -\cos(\mbox{latitude}) \cdot \cos(\mbox{longitude})\\
V_y &=& -\cos(\mbox{latitude}) \cdot \sin(\mbox{longitude})\\
V_z &=& \sin(\mbox{latitude})
\end{eqnarray}$$
Is this part correct? If so, what do I need to add to these equations to account for heading and pitch? My heading notation would be $0$ for pointing north, $\pi$ for pointing south, $\pi/2$ east, $3\pi/2$ for west. For pitch, my notation assumes $0$ means flat, $\pi/2$ straight up, $-\pi/2$ straight down (with respect to the normal plane of the surface of the sphere).
Best Answer
I find Cartesian coordinates let me easily solve a lot of navigation problems on a sphere. I use them even for motion constrained to the surface (so one of the coordinates is always redundant).
For notational convenience, let's write $r$ for radius, $\theta$ for latitude, and $\varphi$ for longitude. We want a new basis consisting of three orthogonal unit vectors, which we'll call $\hat r$, $\hat \theta$, and $\hat \varphi$, oriented according to the "up" direction at the point with spherical coordinates $(r, \theta, \varphi)$ and Cartesian coordinates $(x,y,z)$.
The unit vector "straight up", $\hat r$, is the same as the unit vector from the origin in the direction to your satellite's location, $$ \hat r = \begin{pmatrix} \sin\theta \, \cos\varphi \\ \sin\theta \, \sin\varphi \\ \cos\theta \end{pmatrix}. $$
The unit vector $\hat\varphi$ points "due east". The satellite is over the meridian at longitude $\varphi$; vector $\hat\varphi$ is perpendicular to the plane of that meridian and is equal to a unit vector from the center of the earth to the point on the equator at longitude $\varphi + \frac\pi2$, so $$ \hat \varphi = \begin{pmatrix} -\sin\varphi \\ \cos\varphi \\ 0 \end{pmatrix}. $$
The third vector, $\hat\theta$, is the "northward" vector. (When working in a coordinate system with $\theta = 0$ at the pole, I think it may be more usual to define $\hat \theta$ in the exact opposite direction, that is, pointing due "south", in the direction of increasing $\theta$, but since it's also more common to measure heading relative to due north than relative to due south, a northward vector is a bit more convenient here.) The northward vector $\hat\theta$ is the cross-product $\hat r \times \hat \varphi$ and can be immediately computed from the other two vectors using that formula. Another way to find the coordinates of $\hat\theta$ is to consider it as a vector an angle $\frac\pi2$ "northward" from $\hat r$. That is, it is parallel to a vector pointing from the center of the earth to a point at longitude $\varphi$ and latitude $\theta - \frac\pi2$. When $\theta \geq \frac\pi2$ it is easy to see where this point is, but we can also plot that point when $\theta < \frac\pi2$; in the latter case it is the same as the point at longitude $\varphi+\pi$, latitude $\frac\pi2 - \theta$. In any case, the unit vector's coordinates are $$ \hat \theta = \begin{pmatrix} -\cos\theta \, \cos\varphi \\ -\cos\theta \, \sin\varphi \\ \sin\theta \end{pmatrix} $$ (the same thing you wrote as a vector pointing due north, but in slightly more compact notation).
A much more succinct and elegant derivation of $\hat r$, $\hat\varphi$, and $\hat\theta$ is in this answer to a similar question. (I discovered this some time after initially writing up this answer, and decided to keep the long version here in the hope that it may help some people visualize the geometry.)
If you have the satellite's position in Cartesian coordinates, you can find $\hat r$, $\hat \theta$, and $\hat \varphi$ as follows. Using the fact that $$ \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} r \sin\theta \, \cos\varphi \\ r \sin\theta \, \sin\varphi \\ r \cos\theta \end{pmatrix}, $$ if $x^2 + y^2 \neq 0$ compute \begin{align} r &= \sqrt{x^2 + y^2 + z^2} \\ \hat r &= \begin{pmatrix} x/r \\ y/r \\ z/r \end{pmatrix} \\ \rho &= r \sin\theta = \sqrt{x^2 + y^2} \\ \hat \varphi &= \begin{pmatrix} -y/\rho \\ x/\rho \\ 0 \end{pmatrix} \\ \hat \theta &= \begin{pmatrix} -xz/(\rho r) \\ -yz/(\rho r) \\ \rho/r \end{pmatrix} \\ \end{align}
If $x^2 + y^2 = 0$ then the satellite is on the north-south axis and "heading" has no meaning, so heading-pitch-magnitude will not give us sufficient information to determine the direction of motion.
To convert a heading $\psi$ and pitch $\alpha$ (also known as a bearing or azimuth $\psi$ and elevation angle $\alpha$) to a vector, I would construct the vector using the three unit vectors $\hat \theta$, $\hat \varphi$, and $\hat r$ as components, much as a vector was derived from elevation and azimuth in this answer. That is, letting $v$ be your velocity vector, and assuming you know the magnitude $\|v\|$ of your velocity, $$ v = \|v\| ((\cos\psi\,\cos\alpha)\hat \theta + (\sin\psi\,\cos\alpha)\hat \varphi + (\sin\alpha)\hat r). $$